ESTIMAÇÃO ROBUSTA EM PROCESSOS PERIÓDICOS AUTO …na estimativa dos parâmetros de modelos Auto-Regressivos Integrados de Médias Móveis (ARIMA). Por exemplo, Ledolter (1989) mostrou

ESTIMAÇÃO ROBUSTA EM PROCESSOS PERIÓDICOSAUTO-REGRESSIVOS NA PRESENÇA DE OUTLIERS

ADITIVOS

ALESSANDRO JOSÉ QUEIROZ SARNAGLIAE

VALDÉRIO ANSELMO REISEN

Resumo. Este artigo propõe uma metodologia de estimação dos parâ-metros do modelo Periódico Auto-Regressivo (PAR), a qual é robustana presença de observações atípicas ou outliers. O método de esti-mação é uma variação das equações de Yule-Walker periódicas (McLeod(1994)). Em particular o artigo deriva uma modi�cação da função deautocovariância robusta proposta por Ma & Genton (2000) para esti-mar a função de autocovariância periódica do processo PAR. A estima-tiva dessa função é substituída nas equações de Yule-Walker periódicasfornecendo a estimativa dos parâmetros do modelo PAR. Resultados deMonte Carlo mostraram que, em geral, o estimador proposto para os pa-râmetros auto-regressivos periódicos é robusto na presença de outliersaditivos.

Keywords. Outliers aditivos, modelo PAR, periodicidade, ro-bustês.

1. Introdução

Frequentemente, séries temporais econômicas, hidrológicas, climatológicase de poluição atmosférica contêm observações atípicas ou outliers, geralmenteprovocadas por eventos externos. Diversos autores têm investigado seu efeitona estimativa dos parâmetros de modelos Auto-Regressivos Integrados deMédias Móveis (ARIMA). Por exemplo, Ledolter (1989) mostrou que os in-tervalos de previsão são consideravelmente sensíveis à presença de outliersaditivos, diferente das previsões pontuais (a menos que o(s) outlier(s) este-jam próximos do �nal da série); Chang, Tiao & Chen (1988) e Chen & Liu(1993) mostraram que as estimativas dos parâmetros do modelo ARMA sãomais viesados quando os dados contêm outliers; Deutsch, Richards & Swain(1990) e Chan (1992, 1995) derivaram o vício das autocorrelações amostraisdevido à presença de outliers.

Por outro lado, um fenômeno presente em séries temporais dessas mes-mas áreas é a correlação periódica. Ultimamente, vários modelos capazes decaptar a correlação periódica têm sido propostos na literatura. Um modelobastante utilizado é o modelo Periódico Auto-Regressivo (PAR) (McLeod(1994)) que é uma generalização do modelo Auto-Regressivo (AR) (Box &

1

2 SARNAGLIA, A.J.Q. E REISEN, V.A.

Jenkins (1976)). O processo PAR é caracterizado pela variação dos parâ-metros no tempo. Lund & Basawa (2000) mostraram que o modelo PARpode ser representado na forma vetorial através do modelo Vetorial Auto-Regressivo (VAR). A condição de estacionariadade para o processo VAR édescrita detalhadamente por Lütkepohl (2005).

Franses & Paap (1999) utilizaram o processo PAR para modelar váriasséries de consumo do Reino Unido. McLeod (1993) faz comparações entrealgumas séries de �uxo de rios. Outras aplicações do modelo PAR são ap-resentadas em Bloom�eld, Hurd & Lund (1994), Lund, Hurd, Bloom�eld &Smith (1995) e Gardner & Franks (1975). Para recentes discussões teóricassobre modelos periódicos pode-se citar Lund & Basawa (1997) e Lund &Basawa (2000).

Hurd & Gerr (1991) e Vecchia & Ballerini (1991) apresentam procedi-mentos para detecção de correlação periódica em séries temporais. Existemvárias formas de estimar os parâmetros auto-regressivos do modelo PAR.McLeod (1994) sugere substituir as autocorrelações amostrais periódicas nasequações de Yule-Walker periódicas e então resolver o sistema de equaçãoencontrando assim as estimativas dos parâmetros do modelo.

Este artigo propõe um estimador robusto para os parâmetros auto-regressivosdo modelo PAR, na presença de outliers aditivos. O procedimento de esti-mação é baseado no estimador robusto da função de autocovariância propostopor Ma & Genton (2000), do qual deriva-se um estimador robusto da funçãode autocovariância periódica. Esse estimador é substituído nas equaçõesde Yule-Walker periódicas para obter estimadivas robustas dos coe�cientesauto-regressivos do modelo. Nós mostramos que o nosso estimador é robustoatravés de simulações de Monte Carlo sob vários cenários.

O artigo é estruturado da seguinte forma: Na Seção 2 nós apresentamosa função de autocovariância e de autocorrelação periódica, os processos pe-riodicamente correlacionados, o processo periódico com outlier aditivo e omodelo PAR; Na Seção 3 são apresentadas as metodologias usuais de iden-ti�cação e estimação do modelo PAR; Na Seção 4 é apresentado o impactoque as funções de autocovariância periódica, tanto teórica quanto amostral,sofrem na presença de outliers aditivos; Na Seção 5 é proposto um estimadorrobusto para a função de autocovariância periódica e para os parâmetrosauto-regressivos do modelo PAR; Na Seção 6 são apresentados resultadosempíricos de ensaios de Monte Carlo; Por �m, na Seção 7 são apresentadasalgumas considerações �nais.

2. Processos periódicos

Seja {Yt}, t ∈ Z, um processo estocástico com características sazonais deperíodo s. O índice t pode ser escrito, por meio de divisão inteira, comot = t(r,m) = (r − 1)s + m onde m = 1, . . . , s e r = 1, 2, . . .. Por exemplo,no caso de dados mensais, m seria o mês e r o ano.

ESTIMAÇÃO ROBUSTA EM PROCESSOS PERIÓDICOS AUTO-REGRESSIVOS 3

As Funções de Autocovariância (ACV) e de Autocorrelação (ACF) sãoferramentas de grande importância em análise de séries temporais. A seguirserão de�nidas as extensões da ACV e da ACF usuais para o processo{Yt(r,m)}, que são chamadas de Função de Autocovariância Periódica (PeACV)e de Função de Autocorrelação Periódica (PeACF).

DEFINIÇÃO 2.1. Seja {Yt(r,m)}, um processo tal que V (Yt(r,m)) < ∞,∀ t(r,m) ∈ Z. A PeACV e a PeACF de {Yt(r,m)} são de�nidas, respectiva-mente, como sendo

γ(r,m)Y (h) = Cov

[Yt(r,m), Yt(r,m)−h

]= E[Yt(r,m)Yt(r,m)−h]− E[Yt(r,m)]E[Yt(r,m)−h]

= E[Yt(r,m)Yt(r,m)−h]− µr,mµr,m−he

ρ(r,m)Y (h) =

γ(r,m)Y (h)√

γ(r,m)Y (0)γ

(r,m−h)Y (0)

,

onde E[Yt(r,m)] = µr,m.

Note que a de�nição 2.1 é bastante ampla, entretanto, existem proces-sos estocásticos com propriedades que requerem características especiais daPeACV, tais como os processos periodicamente estacionários, de�nidos aseguir.

DEFINIÇÃO 2.2. Seja {Yt(r,m)} um processo estocástico com PeACV de-

notada por γ(r,m)Y (h). {Yt(r,m)} é dito ser periodicamente estacionário se

µr,m = µm e γ(r,m)Y (h) = γ

(m)Y (h) (1)

existirem e dependerem somente da defasagem (ou "lag") h e do período m,isto é, não dependerem de r.

Note que, de acordo com a de�nição 2.2, se {Yt(r,m)} é um processo pe-

riodicamente estacionário então γ(m)Y (h) e ρ(m)

Y (h) são funções periódicas de

período s, ou seja, γ(m)Y (h) = γ

(m+ks)Y (h) e ρ(m)

Y (h) = ρ(m+ks)Y (h), ∀ k ∈ Z.

Se s = 1 então a condição de estacionariedade periódica equivale a esta-cionariedade usual de processos homogêneos (Tiao & Grupe (1980)). Noque segue, neste artigo, serão considerados apenas os processos periodica-mente estacionários.

A in�uência de outliers na estimação e inferência de parâmetros, em qual-quer análise de dados, já é bastante discutida na literatura, por exemplo, verMartin & Yohai (1985) e Bustos & Yohai (1986). Existem vários tipos deoutliers que provocam efeitos diferentes no processo, entretanto, em geral, ostrês tipos a seguir são os mais considerados (Denby & Martin (1979)): out-liers de inovação (IO), que têm efeito em todas observações subsequentes;outliers aditivos (AO) e outliers de reposição (RO), que têm efeito apenas


na observação afetada. Os AO são conhecidos por provocarem efeitos muitomais prejudiciais do que os IO e possuirem o mesmo efeito que os RO (Ma &Genton (2000)). Por esse motivo, neste artigo serão considerados apenas osoutliers aditivos. A seguir será apresentada a de�nição do processo periódicona presença de outliers aditivos.

DEFINIÇÃO 2.3. Seja {Yt(r,m)} um processo periodicamente estacionário

e Vt(r,m)i , i = 1, . . . , l, l variáveis aleatórias (v.a.) independentes entre si e

independentes de {Yt(r,m)}, com P (Vt(r,m)i = −1) = P (V

t(r,m)i = 1) = pi/2

e P (Vt(r,m)i = 0) = 1 − pi. {Zt(r,m)} é de�nido como sendo um processo

periódico com outlier aditivo se é solução da equação

Zt(r,m) = Yt(r,m) +l∑

i=1

(ωiV

t(r,m)i

)(2)

onde l é a quantidade de outliers, e o parâmetro desconhecido ωi é a magni-tude do i-ésimo outlier.

Note que, de acordo com a de�nição 2.3, a presença de outliers em Zt(r,m)

é aleatória. Essa suposição é bastante plausível já que, em geral, a ocorrênciade observações atípicas em séries temporais não é determinística.

Vários modelos capazes de captar e modelar séries provenientes de pro-cessos periodicamente estacionários têm sido estudados e aplicados à dadosreais de diversas áreas da ciência, por exemplo ver Bloom�eld et al. (1994),Lund et al. (1995) e Gardner & Franks (1975). O processo PAR, de�nidoa seguir, é uma generalização do processo Auto-Regressivo (AR) (Box &Jenkins (1976) e Brockwell & Davis (1991)) e é amplamente utilizado paramodelar séries temporais periódicas.

DEFINIÇÃO 2.4. O processo periodicamente estacionário {Yt(r,m)} é ditoser um processo PAR se e somente se {Yt(r,m)} é qualquer solução da equaçãode diferenças dada por

Yt(r,m) = µm +

pm∑i=1

{φ(m)i

[Yt(r,m)−i − µm−i

]}+ εt(r,m) (3)

onde εt(r,m) ∼ iid(0, σ2(m)), pm é a ordem do polinômio auto-regressivo no

período m, e φ(m)j são os parâmetros auto-regressivos do processo.

Note que, equivalentemente à propriedade de γ(m)Y (h) e ρ(m)

Y (h) (de�nição

2.2), os parâmetros φ(m)j na Equação 3 são funções periódicas de período s,

isto é, φ(m)j = φ

(m+ks)j , j = 1, . . . , pm e k ∈ Z.

Segundo Lund & Basawa (2000), o processo PAR apresentado na de�nição2.4 pode ser representado como um processo AR vetorial (VAR) da forma(

Φ0 − Φ1B − · · · − ΦPBP) (−→

Y r −−→µ)

=−→Υ t


onde, B é o operador de retardo (BkYt = Yt−k),−→Y r =

⟨Yt(r,1), . . . , Yt(r,s)

⟩′,

−→µ = 〈µ1, . . . , µs〉′ e−→Υ t =

⟨εt(r,1), . . . , εt(r,s)

⟩′(′ denota a matriz transposta).

A ordem é P = dp/se, p = max1≤m≤s(pm) e dxe denota o menor inteiromaior ou igual a x. Por simplicidade, neste artigo, será assumido que pm = p,m = 1, . . . , s. As matrizes coe�cientes s× s têm entradas

(Φ0)i,j =

1, se i = j;0, se i < j;

φ(i)i−j , se i > j.

; (Φk)i,j = φ(i)ks+i−j , se 1 ≤ k ≤ P .

Portanto, pelas características de processos vetoriais (Lütkepohl (2005)),{Yt(r,m)} é um processo estacionário se todas as raízes x∗ da equação

0 = det[Φ0 − Φ1x− · · · − ΦPx

P]

= det[Φ0λ

P − Φ1λP−1 − · · · − ΦP

] (4)

estão situadas fora do círculo unitário, onde x é um escalar e λ = 1/x, istoé, se x∗1, . . . , x

∗q , são as q raízes distintas da Equação 4, então {Yt(r,m)} é

estacionário quando |x∗j | > 1, j = 1, . . . , q e (x∗1, . . . , x∗q) ∈ Cq. Note que a

condição |x∗j | > 1 equivale à |λ∗j | < 1, ou seja, a condição de estacionariedadepara o processo {Yt(r,m)} é que as raízes |λ∗j | < 1, j = 1, . . . , q.

3. Procedimento usual para estimação dos parâmetros e

identificação do modelo

O ajuste de um modelo ao conjunto de dados deve seguir três etapas im-portantes: identi�cação, estimação e diagnóstico do ajuste. Este artigo abor-dará as etapas de identi�cação e estimação do processo PAR. A seguir serãosumarizadas as teorias referentes à estimação das funções PeACV e PeACFe dos parâmetros do modelo PAR e à identi�cação da ordem do modelo. Osdesempenhos dos métodos de estimação e identi�cação foram investigados,através de ensaios de Monte Carlo, e os resultados estão apresentados naSeção 6.

3.1. Estimação das Funções PeACV e PeACF. Como mencionado naseção anterior, as funções PeACV e PeACF são ferramentas muito impor-tantes, tanto na estimação quanto na identi�cação do modelo. Entretanto,ambas funções, em geral, não são conhecidas e precisam ser estimadas. Seja{yt(r,m)}, r = 1, . . . , N e m = 1, . . . , s, uma série temporal de tamanho n ob-servada de um processo periodicamente estacionário. Neste artigo, {yt(r,m)}é uma série gerada pelo processo PAR. Por simplicidade, será assumido quen/s = N ∈ Z, isto é, a série possui s períodos com N observações cada.McLeod (1994) sugere estimar as funções PeACF e PeACV utilizando as


funções PeACF e PeACV amostrais dadas, respectivamente, por

ρ̂(m)y (h) =

γ̂(m)y (h)√

γ̂(m)y (0)γ̂

(m−h)y (0)

e

γ̂(m)y (h) =

1

N

N∑r=r∗

{(yt(r,m) − y(m)

)(yt(r,m)−h − y(m−h)

)},

onde y(m) = 1N

∑Nr=1 yt(r,m) é a média amostral do período m, e r∗ é o

menor r tal que t(r,m) − h > 0. Os estimadores das funções PeACV ePeACF mantêm as mesmas propriedades dos estimadores das funções ACVe da ACF para séries temporais não-periódicas (Priestley (1981)).

3.2. Estimação do Modelo. Suponha que a ordem correta do modelo PARseja conhecida a priori. Existem vários procedimentos de estimação dos pa-râmetros do modelo PAR. McLeod (1994) sugere estimar os parâmetros φ(m)

j ,j = 1, . . . , p e m = 1, . . . , s através das equações de Yule-Walker periódicasque são dadas por

pm∑i=1

φ(m)i γ

(m−i)Y (k − i) = γ

(m)Y (k), k = 1, . . . , pm. (5)

Esse procedimento de estimação consiste em substituir, na Equação 5,γ(m)Y (h) pelas suas estimativas γ̂(m)

Y (h), e as soluções dos sistemas de equaçõesresultantes, que são facilmente implementadas computacionalmente, são asestimativas dos parâmetros.

Note que, em geral, a ordem do modelo PAR não é conhecida e isto requerum procedimento de identi�cação da ordem do mesmo. A metodologia paratratar essa problemática é sumarizada na próxima subseção.

3.3. Identi�cação do Modelo. A identi�cação da ordem do modelo PARpode ser feita através da análise visual de suas funções de autocorrelaçãoperiódica (PeACF) e de autocorrelação parcial periódica (PePACF) e/ouatravés de critérios de informação tais como os de Akaike (AIC) e Bayesiano(BIC) (Akaike (1974) e Schwarz (1978)) que podem ser fatorados obtendoum critério separado para cada período. Os critérios AIC e BIC periódicos(fatorados) são calculados da seguinte forma

AIC =

s∑m=1

AICm, AICm = N ln(σ̂2m)

+ 2p, m = 1, . . . , s

e

BIC =

s∑m=1

BICm, BICm = N ln(σ̂2m)

+ p ln (N) , m = 1, . . . , s

onde σ̂2m = SQR/N , e SQR é a soma de quadrados dos resíduos. O pro-cedimento de obtenção da ordem do modelo PAR consiste em estimar os


parâmetros dos modelos concorrentes através da Equação 5 e adotar o mo-delo que minimize o AIC ou o BIC.

McLeod (1994), através de experiências práticas, conclui que a análisegrá�ca das funções PeACF e PePACF, freqüentemente, indica a mesma or-dem sugerida pelos critérios AIC e BIC. Na seção 6 é investigado, através deestudos de Monte Carlo, o desempenho dos critérios AIC e BIC sob várioscenários.

4. Impacto do outlier aditivo na estrutura de correlação

periódica

Como descrito na Seção 2, outliers podem afetar a estrutura de corre-lação de uma série e, conseqüentemente, a estimação dos parâmetros e aidenti�cação da ordem do modelo. Nesta seção são derivados alguns resul-tados referentes ao efeito de AO na estrutura de correlação de um processoperiódico, {Zt(r,m)}, apresentado na de�nição 2.3.

PROPOSIÇÃO 1. Suponha que {Zt(r,m)} seja um processo periódico comAO, tal como descrito pela Equação 2, então a função PeACV de {Zt(r,m)}é dada por

γ(m)Z (h) =

γ(m)Y (h) +

∑li=1 ωip

2i , se h = 0;

γ(m)Y (h), se h 6= 0.

Algumas conseqüências da Proposição 1 são apresentadas no seguintecorolário.

COROLÁRIO 1. Seja γ(m)Z (h) a PeACV apresentada na Proposição 1 e

suponha que 0 < γ(m)Z (0) <∞. Logo, tem-se que

(1) |ρ(m)Z (h)| ≤ |ρ(m)

Y (h)|;(2) lim

ωi→∞ρ(m)Z (h) = 0, para qualquer i = 1, . . . , l.

A prova desse corolário é obtida imediatamente através do resultado apre-sentado na Proposição 1.

OCorolário 1 mostra que o processo periódico com outlier aditivo, {Zt(r,m)},possui PeACF, ρ(m)

Z (h) menor do que a PeACF de {Yt(r,m)}, ρ(m)Y (h), isto

signi�ca que o AO introduz perda de memória em {Zt(r,m)}.O próximo resultado apresenta o efeito que a ocorrência de um outlier

aditivo acarreta na PeACV Amostral.

PROPOSIÇÃO 2. Seja {zt(r,m)}, r = 1, . . . , N e m = 1, . . . , s uma sérietemporal observada do modelo apresentado na de�nição 2.3, com l = 1, esuponha também que ocorre apenas um outlier no tempo t(m0, r0). Então a


função PeACV amostral de {zt(r,m)} é dada por

γ̂(m)z (h) = γ̂(m)

y (h)± ω

N

[yt(r0,m0)−h − y(m0−h)

]I(m0)(m)

± ω

N

[yt(r0,m0)+h − y(m0+h)

]I(m0)(m−h)

+ω2

Nδ(h)I

(m0)(m) I

(m0)(m−h) + op(N

−1)[I(m0)(m) + I

(m0)(m−h)

]onde

I(v)(u) =

{1, se u ≡ v mod(s);0, c.c.

, δ(h) =

{1, se h = 0;0, se h 6= 0.

e u ≡ v mod(s) signi�ca que u é congruente a v módulo s, isto é, u−v = ks,para algum k ∈ Z.

Os resultados obtidos na Proposição 2 mostram um acréscimo signi�cativona PeACV amostral com a presença de um outlier aditivo na série. Essesresultados também fornecem evidências da redução da correlação devido aoacréscimo na variância amostral pelo termo adicional ω2/N quando m ≡m0 mod(s) e m− h ≡ m0 mod(s).

As Proposições 1 e 2 mostram que outliers aditivos podem afetar a inferên-cia realizada sobre processos periódicos. Os resultados revelam que existemperda de memória na estrutura de correlação serial, o que acarreta na sub-estimação dos parâmetros do modelo PAR. Nesse contexto, faz-se necessárioa utilização de métodos robustos de estimação de modelos para séries tem-porais com outliers. A próxima seção dará enfoque nesse assunto, propondoalgumas metodologias robustas.

5. Procedimento robustos para estimação dos parâmetros e

identificação do modelo

A seção anterior mostrou o efeito do AO na estimação dos parâmetros domodelo. As subseções seguintes sugerem metodologias de estimação robustasdas funções PeACV e PeACF para séries temporais periódicas que possuemoutliers aditivos.

5.1. Estimação das funções PeACV e PeACF através da funçãorobusta de escala Qn′(·). Ma & Genton (2000) propuseram o estimadorrobusto da função de autocovariância dada por

γ̃(h) =1

4

[Q2

n′−h(u+ v)−Q2n′−h(u− v)

](6)

e o estimador robusto da função de autocorrelação de�nida por

ρ̃(h) =Q2

n′−h(u+ v)−Q2n′−h(u− v)

Q2n′−h(u+ v) +Q2

n′−h(u− v), (7)


onde u e v são os vetores das n′ − h observações iniciais e das n′ − h ob-servações �nais, respectivamente. O estimador apresentado acima é baseadona igualdade

Cov(X,Y ) =1

4ab[V ar(aX + bY )− V ar(aX − bY )] (8)

onde a escolha a = 1√V ar(X)

e b = 1√V ar(Y )

é recomendada por Gnanadesikan

& Kettenring (1972).A função robusta de escala Qn′(·), proposta por Rousseeuw & Croux

(1992, 1993), é a k-ésima estatística de ordem das(n′

2

)distâncias {|zi−zj |; i <

j} e pode ser escrita como

Qn′(z) = c{|zi − zj |; i < j}(k) (9)

onde z = (z1, z2, . . . , zn′)′, c é uma constante utilizada para garantir a con-

sistência (c = 2.2191 para dados provenientes de distribuição normal) e

k =

⌊(n′2 )+2

4

⌋+ 1.

A função apresentada na Equação 9 pode ser calculada usando o algoritmoproposto por Croux & Rousseeuw (1992), o qual é computacionalmente e�-ciente. Rousseeuw & Croux (1993) mostraram que o ponto de ruptura as-sintótico de Qn′(·) é de 50%, isto é, a propriedade de robustês de Qn′(·) éobservada, mesmo quando a metade do conjunto de observações seja con-taminada por outliers.Qn′(·) tem sido utilizada em diferentes contextos de modelagens, tais como

análise de regressão e estatística espacial (ver Ma & Genton (2000) e suasreferências).

É possível mostrar que −1 ≤ ρ̃(h) ≤ 1 para todo h. Ma & Genton (2000)estabeleceram as propriedades assintóticas de γ̃(h), tais como consistência enormalidade. Os autores também mostraram que γ̃(h) tem ponto de rupturatemporal máximo de 25% e propuseram uma transformação que assegura àfunção γ̃(h) ser positiva-de�nida, porém a robustês do estimador transfor-mado ainda não foi estabelecida. A função γ̃(h) foi utilizada por Ma &Genton (2000) para a modelagem de processos ARMA.

Neste artigo, a quantidade γ̃(h) será utilizada na estimação do modeloPAR(p). Com base na função de autocovariância robusta γ̃(h) (Equação6), a função de Autocovariância Periódica Robusta (PeACVR) amostral éde�nida da seguinte forma

γ̃(m)y (h) =

1

4

[Q2

N−r∗(u(m) + v(m))−Q2

N−r∗(u(m) − v(m))

](10)

onde u(m) e v(m) são, respectivamente, os vetores(yt(r∗,m)−h, . . . , yt(N,m)−h

)e(yt(r∗,m), . . . , yt(N,m)

)ambos com N − r∗ elementos, e r∗ é o menor inteiro

r tal que t(r,m)− h > 0.De maneira análoga, baseando-se na função de autocorrelação robusta

(Equação 7), a função de Autocorrelação Robusta (PeACFR) amostral é


de�nida como sendo

ρ̃(m)y (h) =

Q2N−r∗(u

(m) + v(m))−Q2N−r∗(u

(m) − v(m))

Q2N−r∗(u

(m) + v(m)) +Q2N−r∗(u

(m) − v(m)). (11)

É possível mostrar que −1 ≤ ρ̃(m)y (h) ≤ 1 para todo h. A função PeACVR

amostral γ̃(m)y (h), dada pela Equação 10, será utilizada para a estimação dos

parâmetros do modelo PAR sendo esse o tópico apresentado na subseção 5.2.

5.2. Estimação robusta do modelo. Este procedimento de estimação ro-busta consiste em substituir, na Equação 5, as funções γ(m)

y (h) por γ̃(m)y (h)

obtendo o seguinte sistema de equações:p∑

i=1

φ(m)i γ̃(m−i)y (k − i) = γ̃(m)

y (k), k = 1, . . . , p. (12)

As estimativas robustas φ̃(m)i dos parâmetros do modelo PAR são as

soluções dos sistemas de equações dados pela Equação 12. Note que, nessesistema, ambos os lados das equações são múltiplas de c2 e, portanto, osistema de equações resultantes não dependerá da constante c.

5.3. Identi�cação do modelo. Por �m, a identi�cação do modelo deveser feita utilizando os critérios AIC ou BIC, aplicando as estimativas obtidasatravés do método robusto citado nas seções anteriores.

O modelo a ser escolhido é o que minimiza o critério utilizado, da mesmamaneira que o procedimento apresentado na seção 3.

A próxima seção apresenta os resultados empíricos, utilizando ensaios deMonte Carlo, que investigam o comportamento dessa metodologia em diver-sos cenários.

6. Estudos de Monte Carlo

Nesta seção as metodologias discutidas anteriormente são analizadas em-piricamente através de ensaios de Monte Carlo, e os resultados são apresen-tados nas subseções 6.1 e 6.2. Na primeira, os critérios de identi�cação AICe BIC são investigados com objetivo de veri�car a proporção de identi�caçãocorreta do modelo. A performance dos procedimentos de estimação do mo-delo é a motivação do estudo apresentado na segunda subseção. Para efeitode comparação é considerada a utilização da metodologia clássica de esti-mação das autocorrelações. No estudo de simulação, considerou-se o modeloPAR(1) com ou sem outliers, onde s = 4 e n = 400. Os outliers possuem asmagnitudes ω = 4 e 7, com probabilidade de ocorrência p = 0.01.

Foram realizadas 10000 replicações de Monte Carlo utilizando a linguagemde programação FORTRAN 90. Os grá�cos foram confeccionados na plata-forma R, versão 2.5.1. Os parâmetros apresentados na seguinte tabela foramselecionados de forma a garantir a estacionariedade periódica dos processose que os valores de λ, raiz da Equação 4, estejam próximos de zero (Modelo1) e próximos da raiz unitária (Modelo 2).


Tabela 1: Parâmetros dos modelos PAR(1) simulados

Parâmetros Modelo 1 Modelo 2

φ(1)1 0.9 1.5φ(2)1 0.8 0.8φ(3)1 0.7 1.2φ(4)1 0.6 0.5λ 0.3024 0.7200

6.1. Frequência de acerto dos critérios AIC e BIC. Estudos similaresrelacionados ao comportamento dos critérios AIC e BIC em séries temporaisperiódicas não contaminadas podem ser encontrados em Franses & Paap(1994). Os autores consideraram que as variâncias e as ordens do modelonão dependiam do período. Neste artigo a suposição de variância constantenão é feita, porém, por simplicidade, as ordens dos modelos estimados foram�xadas para todos os períodos.

Para a investigação empírica da performance dos critérios AIC e BIC,foram estimados modelos de ordens 1 a 4. Na Tabela 2 são apresentadasas frequências de seleção do modelo correto dos critérios AIC e BIC para osprocessos 1 e 2, sendo esses contaminados ou não por outliers. Os critériosAIC e BIC foram calculados utilizando os estimadores das autocovariânciaspelos métodos clássico e robusto (Subseção 5.1).

Tabela 2: Freqüência de acerto dos critérios AIC e BIC estimadas para séries temporaisde tamanho n = 400 geradas pelos modelos 1 e 2, na ausência ou na presença de outliers,

com probabilidade de ocorrência p = 0.01

Usual RobustoModelo Mag. AIC BIC AIC BIC

ω = 0 0.7738 0.9973 0.9119 0.9990ω = 4 0.5601 0.9645 0.8406 0.9877

1 ω = 7 0.2222 0.7701 0.7762 0.9529ω = 0 0.7921 0.9977 0.9981 1.0000ω = 4 0.3185 0.8185 0.9097 0.9770

2 ω = 7 0.0437 0.2466 0.5666 0.7492

Inicialmente é analisado o caso em que as séries são não contaminadas, ouseja, ω = 0. Nesse caso, quando o valor de λ, raiz da Equação 4, está próximode zero (Modelo 1) nota-se que as freqüências de acerto dos critérios BICsão altas e competitivas, porém percebe-se que o critério AIC possui umafreqüência de acerto maior quando utiliza-se a metodologia de estimação daautocovariância robusta. Entretanto, quando o valor de λ está próximo daraiz unitária (Modelo 2), percebe-se que as freqüencias de acerto dos critériosAIC e BIC aumentam com a utilização em ambas metodologias.


Para ω > 0, isto é, a série é contaminada, percebe-se uma redução sig-ni�cativa das freqüências de acerto dos critérios AIC e BIC no caso dametodologia clássica. Esse resultado não é surpreendente e pode ser jus-ti�cado pela discussão teórica dada na seção 4. A redução da proporção deacerto é mais signi�cativa quando o valor de λ está próximo da raiz unitária(Modelo 2). A metodologia robusta mantêm a freqüência de acerto em nívelelevado mesmo na presença de outliers, exceto no Modelo 2 (λ próximo daraiz unitária) quando ω = 7, porém, neste caso, a proporção da acertos aindaé maior do que a metodologia clássica.

Portanto, em geral, nota-se que os critérios AIC e BIC calculados uti-lizando a metodologia sugerida não são afetados por outliers com a mesmaintensidade que a metodologia usual.

6.2. Estimação dos parâmetros. Nesta subseção será investigado o de-sempenho do estimador proposto utilizando como critérios de comparação amédia e o Erro Quadrático Médio (EQM) das estimativas. Essas quantidadessão apresentadas nas Tabelas 3 (Modelo 1) e 4 (Modelo 2).

Tabela 3: Médias e EQM das estimativas dos parâmetros do Modelo 1, para séries detamanho n = 400, na ausência ou na presença de outliers, com probabilidade de ocorrência

p = 0.01

Usual Robustoω Parâmetros Est. EQM Est. EQM

φ(1)1 = 0.9 0.8866 0.0062 0.8803 0.0078φ(2)1 = 0.8 0.7963 0.0043 0.7946 0.0061φ(3)1 = 0.7 0.6968 0.0042 0.6943 0.0061

0 φ(4)1 = 0.6 0.5952 0.0047 0.5908 0.0067φ(1)1 = 0.9 0.8192 0.0175 0.8659 0.0093φ(2)1 = 0.8 0.7506 0.0094 0.7857 0.0069φ(3)1 = 0.7 0.6565 0.0079 0.6829 0.0067

4 φ(4)1 = 0.6 0.5573 0.0084 0.5807 0.0075φ(1)1 = 0.9 0.7267 0.0578 0.8767 0.0086φ(2)1 = 0.8 0.6803 0.0312 0.7907 0.0067φ(3)1 = 0.7 0.5988 0.0235 0.6904 0.0066

7 φ(4)1 = 0.6 0.5038 0.0222 0.5879 0.0068

Primeiramente é analizado o caso em que a série não é contaminada, ouseja, ω = 0. Através da Tabela 3 e da Tabela 4 nota-se que, quando ovalor de λ está próximo de zero (Modelo 1) ambas metodologias, em média,são extremamente competitivas, sendo que o EQM do estimador usual éligeiramente menor que o do robusto. Entretanto, quando o valor de λestá próximo da raiz unitária (Modelo 2), percebe-se que o estimador usualcontinua fornecendo estimativas que aproximam-se muito dos valores reais,enquanto a estimativa robusta perdem qualidade.


Tabela 4: Médias e EQM das estimativas dos parâmetros do Modelo 2, para séries detamanho n = 400, na ausência ou na presença de outliers, com probabilidade de ocorrência

p = 0.01

Usual Robustoω Parâmetros Est. EQM Est. EQM

φ(1)1 = 1.5 1.4756 0.0040 1.3361 0.0302φ(2)1 = 0.8 0.7940 0.0012 0.7858 0.0025φ(3)1 = 1.2 1.1944 0.0016 1.1630 0.0038

0 φ(4)1 = 0.5 0.4915 0.0011 0.4271 0.0070φ(1)1 = 1.5 1.4136 0.0163 1.3136 0.0392φ(2)1 = 0.8 0.7795 0.0022 0.7808 0.0030φ(3)1 = 1.2 1.1654 0.0042 1.1528 0.0052

4 φ(4)1 = 0.5 0.4835 0.0016 0.4244 0.0076φ(1)1 = 1.5 1.3078 0.0672 1.3177 0.0382φ(2)1 = 0.8 0.7527 0.0059 0.7822 0.0031φ(3)1 = 1.2 1.1154 0.0162 1.1541 0.0055

7 φ(4)1 = 0.5 0.4698 0.0031 0.4238 0.0077

No caso em que a série é contaminada, isto é, ω > 0, nota-se no Modelo 1que as estimativas usuais, em geral, afastam-se dos valores reais dos parâme-tros e o EQM aumenta, enquanto a metodologias robusta fornece estimativasque permanecem praticamente inalteradas. No Modelo 2, para afetar signi-�cativamente as estimativas dos parâmetros obtidas através da metodologiaclássica é necessário que o(s) outlier(s) possuam magnitudes elevadas.

A presente investigação empírica mostra que, em geral, quando tem-se acerteza de que a série não possui outliers, �ca evidente que a metodologiausual é mais apropriada para a estimação do modelo. Entretanto, se a sérieé contaminada por outliers deve-se utilizar a metodologia robusta, pois essanão sofre danos signi�cativos provocados por outliers.

7. Conclusões

Este artigo investiga o impacto que outliers aditivos provocam na es-timação e identi�cação do modelo PAR(p). O estudo é feito através daestrutura de autocovariância periódica do modelo. Os resultados teóricosmostram que a presença de outlier acarreta perda de memória do processo,resultado também observado na estimativa da função de autocorrelação pe-riódica. Como conseqüência, as estimativas dos parâmetros são afetadas pelaperda de memória (diminuição em valor absoluto) da correlação existente noprocesso. Os estudos empíricos evidenciaram os resultados teóricos encontra-dos. Diante dessa problemática, foi proposto uma metodologia de estimaçãoalternativa para as funções de autocovariância e autocorrelação periódicaso qual mostrou ser, em geral, robusta na presença de outliers. A investi-gação empírica foi feita em diferentes cenários destacando, efeito do tamanho


amostral, do número de outliers, da magnitude dos mesmos em processos pe-riodicamente estacionários e próximos da região de não-estacionariedade. Ametodologia proposta mostrou ser uma alternativa robusta no contexto deestimação do modelo PAR com outliers aditivos.

Agradecimentos

Os autores agradecem ao apoio �nanceiro da FAPEMIG e do PIBIC-CNPq.


Referências

Akaike, H. (1974), À new look at the statistical model indenti�cation', IEEE Transactionson Automatic Control 19, 716�23.

Bloom�eld, P., Hurd, H. L. & Lund, R. B. (1994), `Periodic correlation in stratosphericozone data', Journal of Time Series Analysis 15, 127�150.

Box, G. E. P. & Jenkins, G. M. (1976), Time Series Analysis: Forecasting and Control,Holden day.

Brockwell, P. J. & Davis, R. A. (1991), Time Series: Theory and Methods, Springer verlag.Bustos, O. H. & Yohai, V. J. (1986), `Robust estimates for ARMA models', Journal of

the American Statistical Association 81, 155�68.Chan, W. (1992), À note on time series model speci�cation in the presence outliers',

Journal of Applied Statistics 19, 117�124.Chan, W. (1995), Òutliers and �nancial time series modelling: a cautionary note', Math-

ematics and Computers in Simulation 39, 425�430.Chang, I., Tiao, G. C. & Chen, C. (1988), Èstimation of time series parameters in presence

of outliers', Technometrics 30, 1936�204.Chen, C. & Liu, L. (1993), `Joint estimation of model parameters and outlier e�ects in

time series', Journal of the American Statistical Association 88, 284�297.Croux, C. & Rousseeuw, P. J. (1992), `Time-e�cient algorithms for two highly robust

estimators of scale', Computational Statistics 1, 1�18.Denby, L. & Martin, R. D. (1979), `Robust estimation of the �rst order autoregressive

parameter', Journal of the American Statistical Association 74, 140�46.Deutsch, S. J., Richards, J. E. & Swain, J. (1990), È�ects of a single outlier on ARMA

identi�cation', Communications in Statistics: Theory and Methods 19, 2207�2227.Franses, P. H. & Paap, R. (1994), `Model selection in periodic autoregressions', Oxford

Bulletin of Economics and Statistics 56, 421�39.Franses, P. H. & Paap, R. (1999), `Forecasting with periodic autoregressive time series

models', In: Econometric Institute Report 9927/A .Gardner, W. A. & Franks, L. E. (1975), Ìntervention analysis with applications to eco-

nomic and enviromental problems', IEEE Transactions on Information Theory 21, 4�14.

Gnanadesikan, R. & Kettenring, J. R. (1972), `Robust estimates, residuals, and outlierdetection with multiresponse data', Biometrics 28, 81�124.

Hurd, H. L. & Gerr, N. L. (1991), `Graphical methods for determining the presence ofperiodic correlation', Journal of the American Statistical Association 12, 337�50.

Ledolter, J. (1989), `The e�ect of additive outliers on the forecast from ARMA models',International Journal of Forecasting 5, 231�240.

Lund, R. B. & Basawa, I. V. (1997), `Modeling and inference for periodically correlatedtime series, To Appear, In: S. Ghosh, ed.', Asymptotics, Nonparametrics and TimeSeries .

Lund, R. B. & Basawa, I. V. (2000), `Recursive prediction and likelihood evaluation forperiodic ARMA models', Journal of Time Series Analysis 20, 75�93.

Lund, R. B., Hurd, H. L., Bloom�eld, P. & Smith, R. L. (1995), `Climatological time serieswith periodic correlation', Journal of Climate 8, 2787�2809.

Lütkepohl, H. (2005), New Introduction to Multiple Time Series Analysis, Springer verlag.Ma, Y. & Genton, M. (2000), `Highly robust estimation of the autocovariance function',

Journal of Time Series Analysis 21, 663�684.Martin, R. D. & Yohai, V. J. (1985), `Robustness in time series and estimating ARMA

models', In: Handbook of Statistics pp. 119�55.McLeod, A. I. (1993), `Parsimony, model adequacy and periodic correlation in time series

forecasting', International Statistical Review 61, 387�93.


McLeod, A. I. (1994), `Diagnostic checking periodic autoregression models with applica-tion', Journal of Time Series Analysis 15, 221�33.

Priestley, M. B. (1981), Spectral Analysis and Time Series, Academic press.Rousseeuw, P. J. & Croux, C. (1992), Èxplicit scale estimators with high breakdown

point', In: Statistical Analyses and Related Methods pp. 77�92.Rousseeuw, P. J. & Croux, C. (1993), Àlternatives to the median absolute deviation',

Journal of the American Statistical Association 88, 1273�1283.Schwarz, G. (1978), Èstimating the dimension of a model', The Annals of Statistics 6, 461�

64.Tiao, G. C. & Grupe, M. R. (1980), `Hidden periodic autoregressive-moving average mod-

els in time series data', Biometrika 67, 365�73.Vecchia, A. V. & Ballerini, R. (1991), `Testing for periodic autocorrelations in seasonal

time series data', Biometrika 78, 53�63.


Anexo

Nesta seção são discutidas as demonstrações das Proposições 1 e 2 apre-sentadas na seção 4.

Prova da proposição 1. A expressão de γ(m)Z (h) é dada por

γ(m)Z (h) = Cov

(Zt(r,m), Zt(r,m)−h

)= Cov

(Yt(r,m) +

k∑i=1

ωiV(i)

t(r,m), Yt(r,m)−h +

k∑j=1

ωjV(j)

t(r,m)−h

)

= Cov(Yt(r,m), Yt(r,m)−h

)+ Cov

(Yt(r,m),

k∑j=1

ωjV(j)

t(r,m)−h

)

+ Cov

(k∑

i=1

ωiV(i)

t(r,m), Yt(r,m)−h

)+ Cov

(k∑

i=1

ωiV(i)

t(r,m),

k∑j=1

ωjV(j)

t(r,m)−h

)

= γ(m)Y (h) +

k∑j=1

ωjCov(Yt(r,m), V

(j)

t(r,m)−h

)+

k∑i=1

ωiCov(V

(i)

t(r,m), Yt(r,m)−h

)

+

k∑i=1

k∑j=1

ωiωjCov(V

(i)

t(r,m), V(j)

t(r,m)−h

)

=

γ(m)Y (h) +

k∑i=1

ω2i pi, se h = 0

γ(m)Y (h), se h 6= 0

Prova da proposição 2. Se ocorre apenas 1 outlier no ano r0 e no período m0

então, primeiramente, tem-se que a média amostral é dada por

z(m) =1

N

N∑r=1

(zt(r,m)

)=

1

N

N∑r=1

(yt(r,m) + I

[T ]

[t(r,m)]

)=

{y(m) ±

ω

N, se m = m0

y(m), se m 6= m0

onde,

I[T ]

[t(r,m)]=

{0, se r 6= r0 ou m 6= m0

1, se r = r0 e m = m0


Para veri�car a expressão de γ̂(m)Z (h) = R

(m)Z (h) é necessário analisar 4

casos, m ≡ m0 mod(s) ou m− h ≡ m0 mod(s) ou ambos m ≡ m0 mod(s)e m − h ≡ m0 mod(s) ou m e m − h não congruentes a m0 módulo s. A

expressão de R(m)Z (h) é então obtida da seguinte maneira:

• Se m e m− h são não congruentes a m0 módulo s, então

R(m)z (h) =

N∑r=r∗

{(zt(r,m) − z(m)

) (zt(r,m)−h − z(m−h)

)}=

N∑r=r∗

{(yt(r,m) ± ωI [T ]

[t(r,m)]− y(m)

)(yt(r,m)−h ± ωI [T ]

[t(r,m)−h] − y(m−h)

)}=

N∑r=r∗

{(yt(r,m) − y(m)

)(yt(r,m)−h − y(m−h)

)}= R(m)

y (h)

• Se m ≡ m0 mod(s), então

R(m)z (h) = R(m0)

z (h) =1

N

N∑r=r∗

{(zt(r,m) − z(m)

) (zt(r,m)−h − z(m−h)

)}=

1

N

N∑r=r∗

{[(yt(r,m) − y(m)

)± ωI [T ]

[t(r,m)]∓ ω

N

] [(yt(r,m)−h − y(m−h)

)± ωI [T ]

[t(r,m)−h]

]}= R(m0)

y (h)± 1

N

N∑r=r∗

{ωI

[T ]

[t(r,m)−h]

[yt(r,m) − y(m)

]}± 1

N

N∑r=r∗

{ω[yt(r,m)−h − y(m−h)

] [I[T ]

[t(r,m)]− 1

N

]}

+1

N

N∑r=r∗

{ω2I

[T ]

[t(r,m)−h]

[I[T ]

[t(r,m)]− 1

N

]}

= R(m0)y (h)± 1

N

N∑r=r∗


] [I[T ]

[t(r,m)]− 1

N

]}

= R(m0)y (h)± ω

N

[yt(r0,m0)−h − y(m0−h)

]∓ ω

N2

N∑r=r∗

[yt(r,m)−h − y(m−h)

]= R(m0)

y (h)± ω

Nyt(r0,m0)−h ∓

ω

Ny(m0−h) ±

ω

N

(1− r∗

N

)y(m0−h) ∓

ω

N2

N∑r=r∗

[yt(r,m0)−h

]


• Se m− h ≡ m0 mod(s) (ou ainda, m ≡ m0 + h mod(s)), então

R(m)z (h) = R(m0)

z (h) =1

N

N∑r=r∗

{(zt(r,m) − z(m)

) (zt(r,m)−h − z(m−h)

)}=

1

N

N∑r=r∗

{[(yt(r,m) − y(m)

)± ωI [T ]

[t(r,m)]

] [(yt(r,m)−h − y(m−h)

)± ωI [T ]

[t(r,m)−h] ∓ω

N

]}= R(m0)

y (h)± 1

N

N∑r=r∗

{ω[yt(r,m) − y(m)

] [I[T ]

[t(r,m)−h] −1

N

]}

± 1

N

N∑r=r∗

{ωI

[T ]

[t(r,m)]


]}+

1

N

N∑r=r∗

{ω2I

[T ]

[t(r,m)]

[I[T ]

[t(r,m)−h] −1

N

]}

= R(m0)y (h)± 1

N

N∑r=r∗


] [I[T ]

[t(r,m)−h] −1

N

]}

= R(m0)y (h)± ω

N

[yt(r0,m0)+h − y(m0+h)

]∓ ω

N2

N∑r=r∗

[yt(r,m) − y(m)

]= R(m0)

y (h)± ω

Nyt(r0,m0)+h ∓

ω

Ny(m0+h) ±

ω

N

(1− r∗

N

)y(m0+h) ∓

ω

N2

N∑r=r∗

[yt(r,m0)+h

]


• Se m ≡ m0 mod(s) e m− h ≡ m0 mod(s) (m ≡ m0 + h mod(s)), então

R(m)z (h) = R(m0)

z (h) =1

N

N∑r=r∗

{(zt(r,m) − z(m)

) (zt(r,m)−h − z(m−h)

)}=

1

N

N∑r=r∗

{[(yt(r,m) − y(m)

)± ωI [T ]

[t(r,m)]∓ ω

N

] [(yt(r,m)−h − y(m−h)

)± ωI [T ]

[t(r,m)−h] ∓ω

N

]}= R(m0)

y (h)± 1

N

N∑r=r∗


] [I[T ]

[t(r,m)−h] −1

N

]}

± 1

N

N∑r=r∗


] [I[T ]

[t(r,m)]− 1

N

]}

+1

N

N∑r=r∗

{ω2

[I[T ]

[t(r,m)]− 1

N

] [I[T ]

[t(r,m)−h] −1

N

]}

= R(m0)y (h)± ω

N

[yt(r0,m0)+h − y(m0+h)

]∓ ω

N2

N∑r=r∗

[yt(r,m) − y(m)

]± ω

N

[yt(r0,m0)−h − y(m0−h)

]∓ ω

N2

N∑r=r∗


]+

1

N

N∑r=r∗

ω2

{I[T ]

[t(r,m)]I[T ]

[t(r,m)−h] −1

N

(I[T ]

[t(r,m)]+ I

[T ]

[t(r,m)−h]

)+

1

N2

}

= R(m0)y (h)± ω

Nyt(r0,m0)+h ∓

ω

Ny(m0+h) ±

ω

N

(1− r∗

N

)y(m0+h) ∓

ω

N2

N∑r=r∗

[yt(r,m0)+h

]± ω

Nyt(r0,m0)−h ∓

ω

Ny(m0−h) ±

ω

N

(1− r∗

N

)y(m0−h) ∓

ω

N2

N∑r=r∗

[yt(r,m0)−h

]+

ω2

Nδ(h)− 2ω2

N2+ω2

N2

(1− r∗

N

)= R(m0)

y (h)± ω

N

[yt(r0,m0)+h + yt(r0,m0)−h − y(m0+h) − y(m0−h)

]+ω2

Nδ(h)

± ω

N

(1− r∗

N

)[y(m0+h) + y(m0−h)

]∓ ω

N2

N∑r=r∗

[yt(r,m0)+h + yt(r,m0)−h

]− ω2

N2

(1 +

r∗

N

)

Documents

ESTIMAÇÃO ROBUSTA EM PROCESSOS PERIÓDICOS AUTO …na estimativa dos parâmetros de modelos Auto-Regressivos Integrados de Médias Móveis (ARIMA). Por exemplo, Ledolter (1989) mostrou