apostila (1)

Prévia do material em texto

Universidade Federal de Mato Grosso
Instituto de Ciências Exatas e da Terra
Departamento de Estatística
Análise de Séries Temporais
Anderson Castro Soares de Oliveira
Juliano Bortolini
2015
SUMÁRIO
1 Introdução 4
2 Análise Exploratória de Séries Temporais 7
2.1 Estacionaridade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.1 Diferença . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.2 Teste de Dickey-Fuller (DF) . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.3 Teste de Dickey-Fuller Aumentado (ADF) . . . . . . . . . . . . . . . . . . . 13
2.2 Função de Autocorrelação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Função de Autocorrelação Parcial . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Modelos de Decomposição . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4.1 Tendência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.1.1Teste de Tendência . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.1.1.1 Teste de Cox-Stuart . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.1.2Removendo Tendência . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4.1.2.1 Método de regressão . . . . . . . . . . . . . . . . . . . . . . 24
2.4.1.2.2 Médias Móveis . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.1.2.3 Diferença . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.2 Sazonalidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.2.1Periodograma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4.2.2Teste de Fisher . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4.2.3Removendo Sazonalidade . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.4.2.3.1 Método de regressão . . . . . . . . . . . . . . . . . . . . . . 32
2.5 Tranformação de Séries Temporais . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3 Modelos de Suavização Exponencial 38
3.1 Suavização Exponencial Simples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2 Suavização Exponencial de Holt (SEH) . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3 Suavização Exponencial de Holt-Winters (SEHW) . . . . . . . . . . . . . . . . . . 43
3.3.1 Suavização Exponencial de Holt-Winters Aditivo (SEHWA) . . . . . . . . . 44
3.3.2 Suavização Exponencial de Holt-Winters Multiplicativo (SEHWM) . . . . . 45
4 Modelos de Box & Jenkins 47
4.1 Modelos para Séries Estacionárias . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.1.1 Modelos Auto-Regressivos (AR) . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.1.2 Modelos de Médias Móveis (MA) . . . . . . . . . . . . . . . . . . . . . . . . 55
4.1.3 Modelos Auto-Regressivos de Médias Móveis (ARMA) . . . . . . . . . . . . 61
4.2 Modelos para Séries Não Estacionárias . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2.1 Modelo auto-regressivo integrado de médias móveis (ARIMA) . . . . . . . . 62
4.2.2 Modelo sazonal auto-regressivo integrado de médias móveis SARIMA . . . . 68
4.3 Metodologia de Box & Jenkins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.3.1 Identificação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.2 Estimação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.3 Verificação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.3.3.1 Função de autocorrelação dos resíduos . . . . . . . . . . . . . . . . . 82
4.3.3.2Teste de Ljung-Box . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.3.4 Previsões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.4 Seleção de Modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.4.1 Critério de informação de Akaike (AIC) . . . . . . . . . . . . . . . . . . . . . 89
4.4.2 Erro quadrático médio de previsão (EQMP) . . . . . . . . . . . . . . . . . . 89
4.4.3 Erro Percentual Médio Absoluto (MAPE) . . . . . . . . . . . . . . . . . . . 89
1
Introdução
Uma série temporal é um conjunto de dados coletados ao longo do tempo. Uma característica
definidora de séries temporais é a ordenação dos dados no tempo. Para muitos conjuntos de
dados, por exemplo, as alturas de um grupo de indivíduos, não importa que ordem os dados são
obtidos ou listados. Mas para os dados de séries temporais, a ordem é crucial. A ordenação dos
dados impõe uma certa estrutura de dados, que caracteriza sua dependência e alterar a ordem
pode mudar o significado dos dados.
Exemplo 1.1: São exemplos de séries temporais:
• monitoramento da frequência cardíaca de uma pessoa;
• leituras horárias de temperatura do ar;
• preço de fechamento diário de ações de uma empresa;
• dados de precipitação mensal;
• números de vendas anuais.
Uma série temporal pode ser caracterizada como:
• Discreta - as observações são tomadas em um conjunto discreto, ou seja, em intervalos fixos
de tempo, cuja representação é dada por Yt.
• Contínua - as observações são obtidas continuamente, isto é, a todo instante ao longo do
tempo, cuja representação é Y (t).
Uma série temporal contínua pode ser discretizada, considerando seus valores registrados em
certos intervalos de tempo. Assim, se consider um intervalo de tempo contínuo [0, T ] pode-se
obter N pontos em que
N =
T
∆t
Isto implica em agregar os valores série num intervalo de tempo. Por exemplo observar a
quantidade de chuva medida diariamente.
Introdução 5
A variável observada Y pode assumir valores discretos ou contínuos, isso não implica que se
Y é contínua a série temporal é continua, a definição do tipo de série refere-se a forma que as
observações são obtidas ao longo do tempo.
Como os dados de séries temporais não são aleatoriamente coletados a partir de uma amostra,
faz-se necessária a utilização de métodos estatísticos específicos para a sua análise. Desta forma,
a análise de séries temporais explora a dependência dos dados. Para explorar essa dependência
duas abordagens podem ser utilizadas: abordagem no domínio do tempo e da abordagem no
domínio da freqüência.
A abordagem no domínio do tempo é geralmente motivada pela suposição de que correlação
entre pontos adjacentes no tempo é melhor explicada em termos de uma dependência do valor
atual com os valores passados. Esta abordagem centra-se na modelagem de algum valor futuro
de uma série temporal como uma função paramétrica dos valores atuais e passados.
Por outro lado, a abordagem no domínio da freqüência assume as principais características
de interesse em análises de séries temporais referem-se a variações periódicas ou sistemáticas
encontradas naturalmente na maioria dos dados. Estas variações periódicas são muitas vezes
causados por fenômenos biológicos, físicos ou ambientais de interesse.
Os principais objetivos da análise de séries temporais são:
• Descrição: consiste na utilização de métodos de análise de séries temporais para descrever
a série temporal por meio de medidas descritivas ou gráficos. Neste contexto descreve-se a
componentes presentes na série.
• Previsão: consiste na utilização dos valores observados para prever os valores futuros de
uma série temporal.
A analise de serie temporal pressupõe que exista um processo estocástico gerador da série,
ou seja, que, a cada possível realização aleatória da variável, esteja associada uma probabilidade
de ocorrência da observação. Assim, os modelos utilizados para descrever séries temporais são
processos estocásticos, isto é, processos controlados por leis probabilísticas.
Uma serie temporal representa uma trajetória entre todas as possíveis. Assim para cada t
fixo, Yt é uma variável aleatória. para cada t, teremos uma média µt e uma variância σ2t . Sua
análise supõe que cada ocorrência Y1, Y2, ....Yt da serie seja obtida aleatoriamente com base em
uma estrutura probabilística. Assim, a série temporal Yt é uma sequência de variáveisaleatórias
conjuntamente distribuídas, que descreve a estrutura probabilística da variável no tempo, e partir
deste função pode-se inferir sobre a probabilidade de ocorrência de um outro futuro valor.
Definição 1.1: Um processo estocástico é uma família Y = {Y (t), t ∈ T}, tal que para cada
t ∈ T , Y (t) é uma variável aleatória definida num espaço de probabilidades (Ω,A, P ).
Seja {Y (t), t ∈ T} um processo estocástico, então:
• A função valor médio é definida por E(Yt) = µt, para qualquer t ∈ T .
• A função variância é definida por V ar(Yt) = σ2t , para qualquer t ∈ T .
• A função covariância é definida por Cov(Yt, Ys) = γ(Yt, Ys), para quaisquer t, s ∈ T
Introdução 6
Definição 1.2: Um processo estocástico é dito estacionário se é invariante no tempo, ou
seja, se a suas característica se mantém constantes ao longo do tempo
A estacionariedade é a suposição mais importante para um processo estocástico. A idéia
básica de estacionariedade é que as leis de probabilidade que atuam no processo não mudam com
o tempo, isto é, o processo mantém o equilíbrio estatístico.
Definição 1.3: Um processo estocástico é dito fracamente estacionário ou estacionário de
2a ordem se, e somente se:
• A média é constante para todo os tempos, isto é, E(Yt) = µt = µ, para qualquer t ∈ T .
• A covariância entre Yt e Ys depende somente da diferença temporal k = |t − s| e não dos
tempos t e s. Assim:
Cov(Yt, Ys) = Cov(Yt, Yt+k) = γ(k)
Num processo fracamente estacionário, fazendo k = 0 temos Cov(Yt, Yt+0) = V ar(Yt) = γ(0),
ou seja, a variância do processo, assim como a média, também é constante. Note que tanto a
média quanto a variância precisam ser finitos.
Definição 1.4: Um processo estocástico é dito gaussiano se para qualquer t1, t2, ..., tn ∈ T
as variáveis aleatórias Yt1 , Yt2 , ..., Ytn tem distribuição normal n-variada.
Definição 1.5: Um processo estocástico {et, t ∈ T} é dito ruido branco se:
• A média e a variância são constantes para todo os tempos, isto é, E(et) = µt = µ e
V ar(et) = σ
2, para qualquer t ∈ T .
• As variáveis aleatórias et são não correlacionadas, isto é, Cov(et, es) = 0 ∀ t 6= s
2
Análise Exploratória de Séries Temporais
A análise exploratória de séries temporais envolve não só uma estimativa média e desvio
padrão, mas também :
• caracterização da série por meio de identificação de padrões não aleatórios na série tempo-
ral,
• verificação da correlação entre as observações e
• verificação da estabilidade da variância.
As ferramentas tais como a função de autocorrelação são importantes para a exibição da
maneira em que o passado continua a afetar o futuro. Outras ferramentas, como o periodograma,
são úteis quando os dados contem oscilações em freqüências específicas.
O primeiro passo para qualquer estudo de séries temporais sempre envolve um exame cuida-
doso dos dados por meio de um gráfico da série. A representação gráfica de uma série temporal
é feita através do gráfico de linha, sendo no eixo horizontal o indicador de tempo e, no eixo
vertical, a variável a ser representada.
Exemplo 2.1 (Precipitação anual da cidade Los Angeles (Cryer e Chan, 2008)): A Figura
2.1 apresenta a série temporal dos valores anuais de precipitação registados na cidade de em
Los Angeles, Califórnia (EUA), ao longo de mais de 100 anos. A série mostra uma variação
considerável na quantidade de chuva ao longo dos anos. Poucos anos são de precipitação baixa,
alguns de alta, e muitos são de valores intermediários. O ano de 1883 foi um ano excepcionalmente
úmido para Los Angeles, enquanto 1988 foi bastante seco.
# install.packages("TSA") # instalar pacote "TSA" - Time Series Analysis.
library(TSA) # carregar pacote "TSA"
data(larain) # carregar dados "larain"
# ?larain # descrição dos dados "larain"
plot(larain,ylab="Polegadas",xlab="Ano",type="l") # gráfico da série temporal
Análise Exploratória de Séries Temporais 8
Ano
Po
leg
ad
as
1880 1900 1920 1940 1960 1980
10
20
30
40
Figura 2.1: Precipitação anual, em polegadas, da cidade de Los Angeles, Califórnia (EUA), nos
anos de 1878 a 1992.
Exemplo 2.2 (Temperaturas médias mensais da cidade Dubuque (Cryer e Chan, 2008)):
As temperaturas médias mensais (em graus Fahrenheit), de janeiro de 1964 a dezembro de 1975,
registrado em Dubuque, Iowa (EUA), são mostradas na Figura 2.2.
# library(TSA) # carregar pacote "TSA"
data(tempdub) # carregar dados "tempdub"
# ?tempdub # descrição dos dados "tempdub"
plot(tempdub,ylab="Temperatura",xlab="Ano",type="o")
Ano
Te
mp
era
tur
a
1964 1966 1968 1970 1972 1974 1976
10
20
30
40
50
60
70
Figura 2.2: Temperaturas médias mensais (grau Fahrenheit) de Dubuque, Iowa (EUA), de janeiro
de 1964 a dezembro de 1975.
O gráfico da série temporal permite fazer uma análise prévia, no intuito de compreender o
mecanismo gerador da série, que auxilia a formulação de um modelo matemático que explica os
dados. Em geral avalia-se a existência de três possíveis componentes determinísticas, a tendência,
sazonalidade e ciclo, e uma componente aleatória.
Análise Exploratória de Séries Temporais 9
A construção de gráfico na análise de séries temporais é uma ferramenta importante. Por meio
dele é possível identificar características inerentes aos dados, como variabilidade, observações
atípicas, sazonalidade, tendência, dentre outras. Mas, como os procedimentos visuais nem sempre
são confiáveis, existem testes para confirmar a presença destas componentes.
A componente aleatória, ou variações aleatórias, refere-se aos valores da série que flutuam
em torno de uma média constante. A existência das componentes determinísticas fazem com que
não se perceba a componente aleatória.
A tendência é caracterizada como aquele movimento regular e continuo de longo prazo, re-
fletindo um movimento ascendente ou descendente em longo período de tempo.
Figura 2.3: Exemplo de série temporal com tendência crescente (A) e decrescente (B).
A sazonalidade pode ser entendida como variações periódicas ou cíclicas que se repetem em
intervalos relativamente constantes de tempo. A sazonalidade pode ocorrer em um periodo anual,
mensal, diário, etc.
Figura 2.4: Exemplo de série temporal com sazonalidade (A) e com sazonalidade e tendência (B)
Análise Exploratória de Séries Temporais 10
O Ciclo é caracterizado pela existência de variações ascendentes e descendentes, porém, em
intervalos não regulares de tempo. São as flutuações de longo prazo em torno da curva de
tendência.
Figura 2.5: Exemplo de série temporal com Ciclo
Exemplo 2.3: Na figura 2.6 é apresentado a produção anual de pescado no Brasil no período
de 1950 a 2012. Em que é possível observar que no decorrer dos anos houve um aumento gradual
da produção, indicando a existência de tendência.
Figura 2.6: Produção anual de pescado no Brasil no periodo de 1950 a 2012
Exemplo 2.4: Na figura 2.7 produção mensal de álcool nível 80 no período de janeiro de
1991 a dezembro de 2005. Em que é possível uma oscilação da serie, indicando a existência de
sazonalidade.
Análise Exploratória de Séries Temporais 11
Figura 2.7: dução mensal de álcool nível 80 - índice (média 2002=100), no período de janeiro de
1991 a dezembro de 2005
2.1 Estacionaridade
Uma característica importante de uma série temporal é que os dados são correlacionados, no
sentido que Yt pode ser influenciada por Yt−1, Yt−2, .... Uma das suposições mais frequentes que
se faz a respeito de uma série temporal é de que ela é estacionária, ou seja, ela se desenvolve
aleatoriamente no tempo ao redor de uma média constante, refletindo alguma forma de equilíbrio
estável. De um modo geral uma série será estacionária se não houver mudanças sistemáticas na
média e na variância. Em outras palavras, as propriedades de uma parte das observações são
muito parecidos com os de qualquer outra parte.
Na maioria das séries é observado um comportamento não estacionário, provocado pela pre-
sença de componentes determinísticas. Séries temporais sazonais ou com tendência são exemplos
de séries temporais com comportamento nãoestacionário. Assim, a maior parte dos métodos
que trata de séries temporais não-estacionárias está baseada em técnicas para remover ou filtrar
a parte não-estacionária, deixando apenas a parte que pode ser tratada como estacionária.
2.1.1 Diferença
Definição 2.1: O operador de retardo B é uma função matemática aplicada para um valor
de uma série é obtido a mesma série retardada de um período. Isto é:
BYt = Yt−1
B2Yt = Yt−2
BnYt = Yt−n
Análise Exploratória de Séries Temporais 12
Definição 2.2: O operador diferença é dada por
∆ = Yt − Yt−1 = (1−B)Yt
.
Um outro procedimento para eliminar componentes determinísticas é por meio de diferen-
ças. Este método consiste em tomar diferenças sucessivas da série original, até obter uma série
estacionária.
A primeira diferença de uma série é definida por:
∆Yt = Yt − Yt−1
A segunda diferença é dada por:
∆2Yt = ∆(∆Yt) = ∆(Yt − Yt−1) =
= ∆Yt −∆Yt−1 = (Yt − Yt−1)− (Yt−1 − Yt−2)
= Yt − 2Yt−1 + Yt−2
Assim, a n-ésima diferença é dada por:
∆nYt = ∆(∆
n−1Yt) = (1−B)nYt
2.1.2 Teste de Dickey-Fuller (DF)
O teste de de Dickey-Fuller verifica se uma série é não estacionária. Este teste se baseia na
idéia da existência de uma raiz unitária.
Considerando o seguinte modelo:
Yt = ρYt−1 + et − 1 ≤ ρ ≤ 1
em que et e um termo de erro de ruido branco.
Este modelo será estacionário se |ρ| < 1. Se ρ = 1 (raiz unitaria), o processo gerador da serie
Yt e o passeio aleatório e a serie e não-estacionaria. Este modelo pode ser rescrito da seguinte
forma:
et = Yt − ρYt−1
= Yt − ρBYt
= (1− ρB)Yt
Este modelo será não estacionário se a raiz de (1 − ρB) for 1, por isso a denominação raiz
unitária. Assim, as hipóteses geradas neste teste são:{
H0 : ρ = 1 Série não estacionária
H1 : |ρ| < 1 Série é estacionária
Análise Exploratória de Séries Temporais 13
O teste de Dickey-Fuller trabalha com o modelo na forma de diferenças:
Yt − Yt−1 = ρYt−1 − Yt−1 + et
∆Yt = (ρ− 1)Yt−1 + et
= λYt−1 + et
Assim, as hipóteses podem ser reescritas da seguinte forma:{
H0 : λ = 0 Série não estacionária
H1 : λ < 0 Série é estacionária
Para testar essas hipóteses utiliza-se a estatística:
τc =
λ̂− 1
Sλ
em que:
λ̂ =
n∑
t=2
YtYt−1
n∑
t=2
Y 2t−1
Sλ =
√√√√√√√
n∑
t=2
(Yt − α̂Yt−1)2
(n− 2)
n∑
t=2
Y 2t−1
Rejeita-se H0 se valor− p < α ou τc > τ , em que τ é obtido de uma tabela de Dickey-Fuller
2.1.3 Teste de Dickey-Fuller Aumentado (ADF)
Um problema com o teste de Dickey-Fuller é o fato de considerar que existe dependência
apenas em uma observação anterior e os erros et são ruido branco, assim quando os erros são
correlacionados pode haver uma distorção na conclusão do teste, levando a rejeição da hipótese
nula erroneamente. Para contornar esse problema Dickey e Fuller sugeriram uma extensão no
teste trabalhando com um modelo que contém mais de um diferença defasada, este teste foi
denominado teste de Dickey-Fuller aumentado (ADF).
No teste ADF é considerado o seguinte modelo:
Yt = µ+ ρ1Yt−1 + ρ2Yt−2 + ρpYt−p + et
∆Yt = θ0 + λYt−1 +
p∑
i=1
∆βiYt−i + et
em que
• λ =
(
p∑
i=1
ρi − 1
)
• βi =
(
−
p∑
j=1
ρj+1
)
• et e um termo de erro de ruido branco.
As hipóteses do teste ADF são as mesmas do teste original.
Análise Exploratória de Séries Temporais 14
Exemplo 2.5: A série temporal de temperaturas médias mensais da cidade Dubuque (Cryer
e Chan, 2008) é considerada não estacionária de acordo com o teste ADF (valor-p>0,05).
library(TSA) # carregar pacote "TSA"
library(fUnitRoots) # carregar pacote "tseries"
data(tempdub) # carregar dados "tempdub"
k=trunc((length(tempdub)-1)^(1/3)) #determina numero de lags para o teste ADF
adfTest(tempdub,lag=) #Teste de Dickey-Fuller Aumentado (ADF)
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 5
STATISTIC:
Dickey-Fuller: -0.7689
P VALUE:
0.369
Exemplo 2.6: A série temporal de precipitação anual da cidade Los Angeles (Cryer e Chan,
2008) é considerada estacionária de acordo com o teste ADF (valor-p<0,05).
library(TSA) # carregar pacote "TSA"
library(fUnitRoots) # carregar pacote "fUnitRoots"
data(larain) # carregar dados "larain"
k=trunc((length(larain)-1)^(1/3)) #determina numero de lags para o teste ADF
adfTest(larain,lag=k) #Teste de Dickey-Fuller Aumentado (ADF)
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 4
STATISTIC:
Dickey-Fuller: -4.0359
P VALUE:
0.01
2.2 Função de Autocorrelação
A autocorrelação refere-se à correlação de uma série de tempo com o seu próprio valores pas-
sados e os valores futuros. Autocorrelação também é chamado às vezes "correlação defasada"ou
"correlação serial", que refere-se a correlação entre os membros de uma série de números dispostos
no tempo.
Análise Exploratória de Séries Temporais 15
Autocorrelação positiva pode ser considerado uma forma específica de "persistência", uma
tendência para um sistema de permanecer no mesmo estado, de uma observação para a próxima.
Por exemplo, a probabilidade de amanhã ser chuvoso é maior se hoje é chuvoso do que se hoje é
seco.
Definição 2.3: Lag ou defasagem temporal é um período de tempo entre duas observações
de uma série temporal.
Se um processo estocástico é estacionário então sua média e variância são constantes, assim
a covariância ou autocovariância pode ser definida em termos dos lags.
Definição 2.4: Seja {Yt, t ∈ T} um processo estacionário real discreto de média zero, sua
função de autocovariância (facv) é definida por γτ = E[Yt, Yt+τ ].
Proposição 2.1: A função de autocovariância γτ deve satisfazer as seguintes propriedades:
1. γ0 > 0 ;
2. γ−τ = γτ ;
3. |γτ | ≤ γo;
4. γτ é não negativa definida, tal que
n∑
j=1
n∑
k=1
ajakγτj−τk ≥ 0, para quaisquer números reais
a1, . . . , an, e τ1, . . . , τn de T.
Proposição 2.2: A função de autocorrelação do processo é definida por
ρτ =
γτ
γ0
, τ ∈ T
e possui as mesmas propriedades da função de autocovariancia γτ , exeto que agora temos ρ0 = 1.
A função de autocorrelação amostral rk é definida como:
rk =
ck
c0
=
N−k∑
t=1
(Yt − Ȳ )(Yt+k − Ȳ )
N∑
t=1
(Yt − Ȳ )2
Em geral a função de autocorrelação é apresentada num gráfico k versus ck, com linhas
verticais que expressam a magnitude da autocorrelação (Figura 2.8).
Para determinar o numero máximo de lag K a ser utilizado na construção do gráfico de
autocorrelação utiliza-se a expressão
K = 10log(N)
em que N representa o tamanho da serie.
A partir do gráfico da fac pode-se identificar se uma série é aleatória ou não. Em uma série
aleatória, os valores defasados das séries são não correlacionadas e é esperado que rk ∼= 0. Para
N suficientemente grande e sob a hipótese que pk = 0, a distribuição de rk é aproximadamente
Análise Exploratória de Séries Temporais 16
Figura 2.8: Representação da função de autocorrelação amostral
normal, com média igual a zero e variância dada por:
σ̂2(rk) =
1
n
Baseado no resultado do acima, obtém-se um método de avaliar se a autocorrelação rk é
significativamente diferente de zero, em que verifica-se se os valores da autocorrelação estão
foram de um intervalo intervalo de ± 2√
n
.
Exemplo 2.7: Na figura 2.9 observa-se que nos lags 1, 2 e 3 os valores da autocorrelação
estão fora do intervalo de confiança, assim considera-se que existe autocorrelação apenas nestes
lags.
Figura 2.9: Representação da função de autocorrelação amostral
Análise Exploratória de Séries Temporais 17
A inspeção do gráfico da função de autocorrelação permite indentificar alguns padrões:
• Séries aleatórias: Se uma série é completamente aleatória, os valores de rk são estatistica-
mente iguais a zero.
Figura 2.10: Exemplo de uma função de autocorrelação de uma série aleatória
• Séries com alternâncias: Se uma série temporal alterna, com sucessivas observações em
diferentes lados da média geral, então o gráfico da função da autocorrelação também tende
a alternar, entre valores positivos e negativos de autocorrelação.
Figura 2.11: Exemplo de uma função de autocorrelação de uma série com alternâncias
• Série com sazonalidade: Se a sérieapresenta sazonalidade, então o gráfico da função da
autocorrelação também exibirá uma oscilação na mesma freqüência.
Análise Exploratória de Séries Temporais 18
Figura 2.12: Exemplo de uma função de autocorrelação de uma série com sazonalidade
• Séries com tendência: Se a série contém uma tendência, então os valores de de rk não caem
para zero exceto para valores de "lag"(intervalo de tempo) muito altos. Isto ocorre porque
uma observação de um lado da média geral tende a ser seguida por um grande número de
observações do mesmo lado da média por causa da tendência.
Figura 2.13: Exemplo de uma função de autocorrelação de uma série com tendência
Análise Exploratória de Séries Temporais 19
Figura 2.14: Exemplo de uma função de autocorrelação de uma série com tendência e sazonali-
dade
Exemplo 2.8: A série temporal de temperaturas médias mensais da cidade Los Angeles
apresenta sazonalidade, sua função da autocorrelação apresenta uma oscilação (Figura 2.15)
library(TSA) # carregar pacote "TSA"
data(tempdub) # carregar dados "tempdub"
acf(tempdub)
Figura 2.15: Função de autocorrelação da série temporal de temperaturas médias mensais da
cidade Los Angeles
Exemplo 2.9: A série temporal de precipitação da cidade Los Angeles pode ser considerada
aleatória pois sua de autocorrelação não apresenta nenhum lag significativo (Figura 2.16)
Análise Exploratória de Séries Temporais 20
library(TSA) # carregar pacote "TSA"
data(larain) # carregar dados "larain"
acf(larain)
Figura 2.16: Função de autocorrelação da série temporal de precipitação mensais da cidade Los
Angeles
2.3 Função de Autocorrelação Parcial
A autocorrelação parcial mede a autocorrelação entre Yt e Yt+k sem considerar as a autocor-
relação das observações Yt+1 a Yt+k−1.
Definição 2.5: A função de autocorrelaçao parcial (facp) do processo é definida por:
φkk =
|P ∗k |
|Pk|
em que
• Pk é a matriz de autocorrelação dada por
1 ρ1 ρ2 · · · ρk−1
ρ1 1 ρ1 · · · ρk−2
ρ2 ρ1 1 · · · ρk−3
...
...
...
. . .
...
ρk−1 ρk−2 ρk−3 · · · 1

• P ∗k é a matriz Pk com a última coluna substituída pelo vetor de autocorrelação {ρ1, ρ2, ..., ρk}
Análise Exploratória de Séries Temporais 21
Assim, temos:
φ11 = 1
φ22 =
∣∣∣∣∣ 1 ρ1ρ1 ρ2
∣∣∣∣∣∣∣∣∣∣ 1 ρ1ρ1 1
∣∣∣∣∣
=
ρ2 − ρ21
1− ρ21
φ33 =
∣∣∣∣∣∣∣
1 ρ1 ρ1
ρ1 1 ρ2
ρ2 ρ1 ρ3
∣∣∣∣∣∣∣∣∣∣∣∣∣∣
1 ρ1 ρ2
ρ1 1 ρ1
ρ2 ρ1 1
∣∣∣∣∣∣∣
=
(1− ρ21)ρ3 + ρ1ρ22 − 2ρ1ρ2 + ρ31
1− ρ22 + 2ρ21ρ2 − 2ρ21
Da mesma forma que a fac a função de autocorrelação parcial é apresentada num gráfico
k versus φkk, com linhas verticais que expressam a magnitude da autocorrelação parcial. Para
avaliar se a autocorrelação parcial φkk é significativamente diferente de zero, verifica-se se os
valores estão foram de um intervalo intervalo de ± 2√
n
.
2.4 Modelos de Decomposição
Considerando as observações {Yt, t ∈ T} de uma série temporal. Um modelo de decomposição
consiste em escrever Yt como uma soma de três componentes
Yt = Tt + St + at
em que
• Tt é a componente de tendência
• St é a componente cíclica ou sazonal
• at é a componente aleatória
Um série temporal também pode ser representada por meio de um modelo multiplicativo
Yt = TtStat
Removendo-se as componentes Tt e St resta apenas a componente aleatória ou residual at.
Uma suposição que se faz é que at seja um ruído branco, então E[atas] = 0 para t 6= s. Mas
supor apenas que at seja um processo estacionário com média zero e variância σ2.
Assim, o principal interesse é obter uma série estacionária, ou seja, livre de tendência e sazona-
lidade. Quando estas componentes estão presentes na série pode-se obter uma serie estacionária
Análise Exploratória de Séries Temporais 22
estimando-se a tendência e a sazonalidade, e subtraindo suas estimativas da série original
Yt = Tt + St + at série não estacionária
at = Yt − T̂t − Ŝt série estacionária
2.4.1 Tendência
A Tendência em uma série de tempo é uma mudança lenta, gradual, de alguma propriedade
da série em todo o intervalo de observação. Algumas séries apresentam tendência forte e outras
fraca.
Supondo um modelo em que apenas a componente de tendência esteja presente, temos:
Yt = Tt + at
em que
• Tt é a componente de tendência
• at é a componente aleatória
2.4.1.1 Teste de Tendência
2.4.1.1.1 Teste de Cox-Stuart
O teste de Cox-Stuart é aplicável a uma grande variedade de situações e pode ser utilizado
para verificar a presença de tendência. O método proposto baseia-se na distribuição binomial.
Considere um conjunto de observações Y1, Y 2, ..., YN . Considere:
c = N2 se N par
c = N+12 se N impar
Agrupe as observações em pares (Y1, Y1+c), (Y2, Y2+c), ..., (YN , YN+c).
Assim, queremos testar:{
H0 : P (Yi < Yi+c) = P (Yi > Yi+c), ∀i Não exsite tendência
H1 : P (Yi < Yi+c) 6= P (Yi > Yi+c), ∀i Existe tendência
A cada par de observações é associado o sinal
+ se Yi < Yi+c
0 se Yi = Yi+c
− se Yi > Yi+c
Seja T o numero total de pares com sinal + e n o numero de pares onde Yi 6= Yi+c.
Supondo H0 verdadeira T ∼ Binomial(n, 12) assim, rejeita-se H0
• n < 30 se T ≥ n− t, em que t encontrado numa tabela binomial.
• n ≥ 30 se T ≥ 12 + zα2
√
n
4
Análise Exploratória de Séries Temporais 23
Exemplo 2.10: Considere parte da série temporal de produção de pescado no Brasil (Tabela
2.1), temos N = 20, assim c = 10
Tabela 2.1: Produção anual de pescado no Brasil no período de 1993 a 2012
t Yt t Yt
1 676,441 11 985,4
2 701,251 12 1015,9
3 652,91 13 1008,0
4 693,172 14 1050,8
5 732,261 15 1072,2
6 710,704 16 1157,6
7 744,598 17 1241,6
8 839,296 18 1265,5
9 935,946 19 1433,3
10 1003,26 20 1551,2
Assim, considerando os pares (tabela 2.2) temos T = 10 e n = 10, desta forma rejeita-se H0.
Tabela 2.2: Sinal dos pares da série Produção anual de pescado no Brasil no período de 1993 a
2012
Pares Yi Yi + c Sinal
1 676,4 985,4 +
2 701,3 1015,9 +
3 652,9 1008 +
4 693,2 1050,8 +
5 732,3 1072,2 +
6 710,7 1157,6 +
7 744,6 1241,6 +
8 839,3 1265,5 +
9 935,9 1433,3 +
10 1003,3 1551,2 +
T 10
require(xlsx)
require(snpar)
dados=read.xlsx("pescado.xlsx",1)
attach(dados)
y=as.ts(y)
cs.test(e)
Exact Cox-Stuart trend test
data: e
S = 5, p-value = 0.7539
alternative hypothesis: data have a monotonic trend
Análise Exploratória de Séries Temporais 24
2.4.1.2 Removendo Tendência
Assim, se a tendência for removida pode-se obter uma série estacionária. Quatro abordagens
alternativas podem ser utilizadas para remover a tendência:
• Método de regressão
• Médias Móveis
• Diferença
2.4.1.2.1 Método de regressão
A análise de regressão podem ser utilizado para estimar os parâmetros do modelo para ten-
dência. Podem ser utilizado vários tipos de modelos: polinomiais, exponencial, logística, etc. O
objetivo é estimar a tendência como uma função do tempo:
Tt = f(t)
Se for utilizado um modelo polinomial:
Tt = β0 + β1t+ ...+ βmt
m
Exemplo 2.11: Considere parte série temporal de produção de pescado no Brasil (Tabela
2.1), é possível observar que para esse um polinômio de primeiro grau (regressão linear) pode ser
ajustada. Assim,
Tt = β0 + β1t
Assim, ajustando uma regressão linear temos:
T = 527, 3 + 4, 5t
Utilizando o modelo estimado pode-se prever os valores da série com tendência. Assim, a
partir da diferença entre os valores observados e estimados é obtida a série sem tendência dado
pelo erro de previsão, que deverá ser uma série estacionária:
at = Yt − Ŷt
Na tabela 2.3 é apresentado os valores estimado e o erro de previsão da série
Análise Exploratória de Séries Temporais 25
Tabela 2.3: Valores observados, previsto e erro de previsão da série de produção anual de pescado
no Brasil no período de 1993 a 2012
t Yt Ŷt at t Yt Ŷt at
1 676,4 569,8 106,6 11 985,4 994,8 -9,4
2 701,3 612,3 89 12 1015,9 1037,3 -21,4
3 652,9 654,8 -1,9 13 1008 1079,8 -71,8
4 693,2 697,3 -4,1 14 1050,8 1122,3 -71,5
5 732,3 739,8 -7,5 15 1072,2 1164,8 -92,6
6 710,7 782,3 -71,6 16 1157,6 1207,3 -49,7
7 744,6 824,8 -80,2 17 1241,6 1249,8 -8,2
8 839,3 867,3 -28 18 1265,5 1292,3 -26,8
9 935,9 909,8 26,1 19 1433,3 1334,898,5
10 1003,3 952,3 51 20 1551,2 1377,3 173,9
#Estimando a tendência
n=length(y)
t=seq(1,n,1)
model=lm(y~t)
summary(model)
Call:
lm(formula = y ~ t)
Residuals:
Min 1Q Median 3Q Max
-92.617 -55.167 -8.818 32.304 173.886
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 527.326 34.201 15.42 8.13e-12 ***
t 42.499 2.855 14.89 1.47e-11 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 73.63 on 18 degrees of freedom
Multiple R-squared: 0.9249, Adjusted R-squared: 0.9207
F-statistic: 221.6 on 1 and 18 DF, p-value: 1.465e-11
e=resid(model)
e=as.ts(e)
plot(y,xlab="Tempo",ylab="Produção Pesqueira em Milhões de Toneladas")
plot(e)
cs.test(e)
Análise Exploratória de Séries Temporais 26
Exact Cox-Stuart trend test
data: e
S = 5, p-value = 0.7539
alternative hypothesis: data have a monotonic trend
Figura 2.17: Produção anual de pescado no Brasil no período de 1993 a 2012, série com e sem
tendência
2.4.1.2.2 Médias Móveis
As médias móveis é um tipo de suavização de uma série temporal que permite estimar a
tendência. A sua idéia básica é obter um subconjunto de tamanho fixo das primeiras observações
da série, assim a primeira média móvel será dada por pela média desse subconjunto, em seguida é
descolado o subconjunto para frente e calculado novamente sua média. Este processo é repetido
até que todas as observações da série sejam utilizados.
A forma mais simples de obter uma média móvel é por meio da expressão:
Y ∗t =
1
2n+ 1
n∑
j=−n
Zt+j
em que 2n+ 1 representa o tamanho do subconjunto.
Assim, se considerando n = 1, temos um subconjunto de 3 observações da série. Se n = 2
temos um subconjunto de 5 observações da série.
Uma limitação do método de médias móveis é que não é possível fazer estimativas para o
observações nos instantes t = 1, ..., n e t = N − n+ 1.
A partir da diferença entre os valores observados e as médias móveis é obtida a série sem
Análise Exploratória de Séries Temporais 27
tendência que deverá ser uma série estacionária:
at = Yt − Y ∗t
Exemplo 2.12: Na tabela 2.4 é apresentado os valores observados, médias móveis e erro de
previsão da série de produção anual de pescado no Brasil no período de 1993 a 2012.
Tabela 2.4: Valores observados, médias móveis e erro de previsão da série de produção anual de
pescado no Brasil no período de 1993 a 2012
t Yt Y ∗t at
1 676,4 − −
2 701,3 676,9 24,4
3 652,9 682,5 -29,6
4 693,2 692,8 0,4
5 732,3 712,1 20,2
6 710,7 729,2 -18,5
7 744,6 764,9 -20,3
8 839,3 839,9 -0,6
9 935,9 926,2 9,7
10 1003,3 974,9 28,4
11 985,4 1001,5 -16,1
12 1015,9 1003,1 12,8
13 1008 1024,9 -16,9
14 1050,8 1043,7 7,1
15 1072,2 1093,5 -21,3
16 1157,6 1157,1 0,5
17 1241,6 1221,6 20,0
18 1265,5 1313,5 -48,0
19 1433,3 1416,7 16,6
20 1551,2 − −
##Médias Móveis
n=1
j=2*n+1
y1=filter(y,rep(1,j)/j)
e=y-y1
cs.test(e)
Exact Cox-Stuart trend test
data: e
S = 3, p-value = 0.5078
alternative hypothesis: data have a monotonic trend
plot(y,xlab="Tempo",ylab="Produção Pesqueira em Milhões de Toneladas")
plot(e)
Análise Exploratória de Séries Temporais 28
Figura 2.18: Produção anual de pescado no Brasil no período de 1993 a 2012, série com e sem
tendência
2.4.1.2.3 Diferença
Um outro procedimento para eliminar a tendência é por meio de diferenças. Este método
consiste em tomar diferenças sucessivas da série original, até obter uma série estacionária.
Em geral toma-se uma ou duas diferenças para eliminar a tendência.
Exemplo 2.13: Na tabela 2.5 é apresentado os valores observados e a primeira diferença da
série de produção anual de pescado no Brasil no período de 1993 a 2012.
Tabela 2.5: Valores observados e a primeira diferença da série de produção anual de pescado no
Brasil no período de 1993 a 2012
t Yt Yt − Y t− 1 t Yt Yt − Y t− 1
1993 676,4 - 2003 985,4 -17,9
1994 701,3 24,9 2004 1015,9 30,5
1995 652,9 -48,4 2005 1008 -7,9
1996 693,2 40,3 2006 1050,8 42,8
1997 732,3 39,1 2007 1072,2 21,4
1998 710,7 -21,6 2008 1157,6 85,4
1999 744,6 33,9 2009 1241,6 84
2000 839,3 94,7 2010 1265,5 23,9
2001 935,9 96,6 2011 1433,3 167,8
2002 1003,3 67,4 2012 1551,2 117,9
##Diferença
y1=diff(y,1)
cs.test(y1)
Análise Exploratória de Séries Temporais 29
Exact Cox-Stuart trend test
data: y1
S = 2, p-value = 0.1797
alternative hypothesis: data have a monotonic trend
plot(y,xlab="Tempo",ylab="Produção Pesqueira em Milhões de Toneladas")
plot(y1)
Figura 2.19: Produção anual de pescado no Brasil no período de 1993 a 2012, série com e sem
tendência
2.4.2 Sazonalidade
Da mesma forma que a tendência, a sazonalidade (ou periodicidade) constitui uma outra
forma de não-estacionaridade e deve ser estimada e retirada da série. Considera-se como sazonais
os fenômenos que ocorrem com regularidade no tempo. O que se observa em séries sazonais é
que ocorrem relacoções:
• Entre observações para meses sucessivos em um ano particular
• entre as observações do mesmo mês em anos sucessivos.
As séries sazonais são caracterizadas por apresentarem correlação alta em "lag sazonais", isto
é, lags que são múltiplos do período.
A sazonalidade por ser:
• Deterministica: pode ser predita perfeitamente a partir de meses anteriores
• Estocástica: a componente sazonal varia com o tempo.
Sendo Ŝt a estimativa de St, a série sazonalmente ajustada é dada por:
Análise Exploratória de Séries Temporais 30
• Modelo Aditivo
Y SAt = Yt − Ŝt
• Modelo Multiplicativo
Y SAt =
Yt
Ŝt
2.4.2.1 Periodograma
O periodograma é uma descrição dos valores observados numa realização de uma série através
da sobreposição de ondas sinusoidais com várias freqüências. A aplicação prática mais óbvia desta
decomposição é a de servir de instrumento à identificação de componentes cíclicas ou periódicas
Priestley (1989) define o periodograma do processo estacionário (at) como:
Ip(fi) =
2
n
(∑
t=1
)natcos
2πi
n
)2(∑
t=1
)natsen
2πi
n
)2
em que:
• fi é a freqüência analisada em que 0 < fi < 12
• Ip(fi) é a intensidade da freqüência fi
A periodicidade de período é dada por 1fi e pode ser observada pela existência de picos na
freqüência fi = 1n .
No gráfico do periodograma, o periodo 1fi é representada no eixo das ordenadas e a intensidade
da frequência Ip(fi) no das abscissas. Na maioria das vezes, o pico de maior intensidade é o
componente periódico.
2.4.2.2 Teste de Fisher
O teste de Fisher é utilizado para confirmar a existência de sazonalidade em uma série
temporal. Segundo Priestley (1989), esse teste utiliza os valores do periodograma e detecta
grandes periodicidades.
As hipóteses a serem testadas são:{
H0 : Não existe sazonalidade
H1 : Existe sazonalidade
O passos para realização do teste são:
• Roda-se o periodograma, utilizando um pacote estatístico.
• Toma-se a maior periodicidade encontrada no periodograma, max(Ip).
• Calcula-se a estatística g
g =
max(Ip)∑
p=1
)N/2Ip
Análise Exploratória de Séries Temporais 31
em que Ip é o valor do periodograma no período p e N é o número de observações da série
• Calcula-se a estatística do Teste de Fisher, zα dada por:
zα = 1−
(α
n
) 1
n−1
em que n = N2 e α é o nível de significância do teste
• Se g > zα rejeita-se H0, ou seja, a série apresenta periodicidade p.
Pode-se testar o segundo maior periodo de Ip. Usa-se a estatística
g =
I ′′p∑
p=1
)N/2Ip −max(Ip)
sendo I ′′p o segundo maior valor do periodograma.
Nesse caso utiliza-se n = N−12 no calculo de zα.
Exemplo 2.14: Considere parte série temporal de produção de pescado no Brasil (Tabela
2.6).
Tabela 2.6: Produção mensal de álcool nível 80 no período de janeiro de 1991 a dezembro de
1995
Mês Ano
1991 1992 1993 1994 1995
jan 46,19 38,55 36,61 20,10 34,25
fev 31,19 17,63 19,12 12,76 22,49
mar 19,06 7,02 2,66 6,92 11,04
abr 2,38 3,36 3,72 4,93 2,17
mai 101,98 55,46 109,98 123,25 69,47
jun 209,45 200,23 201,54 221,63 205,64
jul 257,16 257,90 249,53 257,34 239,72
ago 269,92 262,06 255,99 268,40 268,87
set 268,54 208,06 221,44 260,94 266,50
out 241,34 230,37 216,55 220,45 216,90
nov 142,30 174,90129,34 120,15 154,13
dez 62,21 80,60 43,40 53,17 76,75
Na figura 2.20 é apresentado o periodograma da série, em que verifica-se a existência de um
período bem definido.
A maior periodicidade ocorre com p = 12 e apresenta intensidade da freqüência igual a
573816,7 e soma das intensidades é dada por 607113,1, assim:
g =
573816, 7
607113, 1
= 0, 9451
O tamanho da série é N = 60, assim n = 30, logo considerando α = 0, 05 então zα será dado
Análise Exploratória de Séries Temporais 32
Figura 2.20: Periodograma para produção mensal de álcool nível 80 no período de janeiro de
1991 a dezembro de 1995
por:
zα = 1−
(
0, 05
30
) 1
30−1
= 0, 1979
Desta forma como g > zα rejeita-se H0, ou seja, a série apresenta periodicidade de 12 meses.
require(xlsx)
dados=read.xlsx("alcool.xlsx",2)
attach(dados)
y=ts(alcool,start=1991,frequency=12)
##Sazonalidade
P=periodograma(y)
Fisher.test(P)
g zalpha valorp periodo
0.9451562 0.1979496 8.162462e-36 12
2.4.2.3 Removendo Sazonalidade
2.4.2.3.1 Método de regressão
A análise de regressão podem ser utilizado para estimar os parâmetros do modelo para sazo-
nalidade. É utilizada um modelo de regressão não linear dado por:
St = µ+
k∑
i=1
[
αicos
(
2πit
p
)
+ βjsen
(
2πit
p
)]
em que pode k ser escolhido entre os valores 1 < k < p
Análise Exploratória de Séries Temporais 33
Exemplo 2.15: Considerando a série apresentada na tabela 2.20), é possivel obter as se-
guintes estimativas para um modelo considerando k = 3, Assim, retirando α2 o modelo para
Tabela 2.7: Estimativas do modelo de regressão para Sazonalidade para produção mensal de
álcool nível 80 no período de janeiro de 1991 a dezembro de 1995
Parâmetro Estimativa valor-p
µ 130,26 < 0, 0001
α1 -49,79 < 0, 0001
β1 -129,03 < 0, 0001
α2 3,54 0,2072
β2 15,67 < 0, 0001
α3 -20,31 < 0, 0001
β3 -5,60 0,0488
sazonalidade é dada por:
St = 130, 26−49, 79cos
(
2πt
12
)
−129, 03sen
(
2πt
12
)
−15, 67sen
(
2πt
6
)
−20, 31cos
(
2πt
4
)
−5, 60sen
(
2πt
4
)
Na tabela 2.8 é apresentado os valores estimado e o erro de previsão da série
Tabela 2.8: Valores observados, previsto e erro de previsão da série produção mensal de álcool
nível 80 no período de janeiro de 1991 a dezembro de 1995
Meses Anos
1991 1992 1993 1994 1995
Ŷt at Ŷt at Ŷt at Ŷt at Ŷt at
jan 30,6 15,59 30,6 7,95 30,6 6,01 30,6 -10,5 30,6 3,65
fev 27,5 3,69 27,5 -9,87 27,5 -8,38 27,5 -14,74 27,5 -5,01
mar 6,83 12,23 6,83 0,19 6,83 -4,17 6,83 0,09 6,83 4,21
abr 9,53 -7,15 9,53 -6,17 9,53 -5,81 9,53 -4,6 9,53 -7,36
mai 89,7 12,28 89,7 -34,24 89,7 20,28 89,7 33,55 89,7 -20,23
jun 200,36 9,09 200,36 -0,13 200,36 1,18 200,36 21,27 200,36 5,28
jul 257,05 0,11 257,05 0,85 257,05 -7,52 257,05 0,29 257,05 -17,33
ago 260,15 9,77 260,15 1,91 260,15 -4,16 260,15 8,25 260,15 8,72
set 253,7 14,84 253,7 -45,64 253,7 -32,26 253,7 7,24 253,7 12,8
out 223,86 17,48 223,86 6,51 223,86 -7,31 223,86 -3,41 223,86 -6,96
nov 143,69 -1,39 143,69 31,21 143,69 -14,35 143,69 -23,54 143,69 10,44
dez 60,16 2,05 60,16 20,44 60,16 -16,76 60,16 -6,99 60,16 16,59
###Estimando a sazonalidade
p=12
n=length(y)
t=seq(1,n,1)
model=nls(y~mu+a1*cos(2*pi*t/12)+a2*cos(4*pi*t/12)+a3*cos(6*pi*t/12)+b1*sin(2*pi*t/12)+
b2*sin(4*pi*t/12)+b3*sin(6*pi*t/12),start=list(mu=10,a1=-1,a2=10,a3=1,b1=10,b2=10,b3=1))
summary(model)
Análise Exploratória de Séries Temporais 34
Parameters:
Estimate Std. Error t value Pr(>|t|)
mu 130.262 1.962 66.391 < 2e-16 ***
a1 -49.785 2.775 -17.942 < 2e-16 ***
a2 3.543 2.775 1.277 0.2072
a3 -20.312 2.775 -7.320 1.38e-09 ***
b1 -129.030 2.775 -46.502 < 2e-16 ***
b2 15.665 2.775 5.646 6.59e-07 ***
b3 -5.595 2.775 -2.016 0.0488 *
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 15.2 on 53 degrees of freedom
Number of iterations to convergence: 1
Achieved convergence tolerance: 7.307e-08
model1=nls(y~mu+a1*cos(2*pi*t/12)+a3*cos(6*pi*t/12)+b1*sin(2*pi*t/12)+b2*sin(4*pi*t/12)+
b3*sin(6*pi*t/12),start=list(mu=10,a1=1,b1=10,b2=10,a3=1,b3=1))
summary(model1)
Parameters:
Estimate Std. Error t value Pr(>|t|)
mu 130.262 1.973 66.007 < 2e-16 ***
a1 -49.785 2.791 -17.839 < 2e-16 ***
b1 -129.030 2.791 -46.233 < 2e-16 ***
b2 15.665 2.791 5.613 7.06e-07 ***
a3 -20.312 2.791 -7.278 1.46e-09 ***
b3 -5.595 2.791 -2.005 0.05 .
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 15.29 on 54 degrees of freedom
Number of iterations to convergence: 1
Achieved convergence tolerance: 4.354e-08
e=resid(model1)
e=as.ts(e,start=1991,frequency=12))
plot(y,xlab="Tempo",ylab="Produção Pesqueira em Milhões de Toneladas")
plot(e,xlab="Tempo",ylab=expression(a[t]))
P=periodograma(e)
Fisher.test(P)
g zalpha valorp periodo
Análise Exploratória de Séries Temporais 35
0.1370854 0.1979496 0.4170492 20
Na figura 2.21 é apresentado o periodograma para o erro do modelo modelo de regressão para
sazonalidade para produção mensal de álcool nível 80 no período de janeiro de 1991 a dezembro
de 1995, em que observa-se vários picos não caracterizando a presença de sazonalidade.
Figura 2.21: Periodograma para o erro do modelo modelo de regressão para Sazonalidade para
produção mensal de álcool nível 80 no período de janeiro de 1991 a dezembro de 1995
Figura 2.22: Produção mensal de álcool nível 80 no período de janeiro de 1991 a dezembro de
1995, série com e sem sazonalidade
Análise Exploratória de Séries Temporais 36
2.5 Tranformação de Séries Temporais
Se a variância da séria aumenta ao longo do tempo, a serie pode apresentar heteroscedas-
ticidade, e conveniente fazer uma transformação a fim de tornar sua variância constante. Para
verificar a necessidade de transformação utiliza-se um gráfico da média versus amplitude. Para
obter esse gráfico deve-se calcular a média e amplitude considerando intervalos regulares da série.
Se a série for sazonal deve-se o intervalo deve ser do tamanho do periodo de sazonalidade, nas
demais séries pode-se utilizar intervalos de 8 a 12 observações. A idéia básica desse gráfico é
colocar em um eixo uma medida de posição e, no outro, uma da variação. Se a amplitude for
independente da média, não há necessidade de transformação. Se for diretamente proporcional
, a transformação sugerida é a logarítmica natural.
Se houver uma tendência na série e do tamanho do efeito sazonal parece aumentar com o
tempo, então os efeitos são multiplicativos, nesse caso é aconselhável transformar os dados. Para
tornar os efeitos aditivos a transformação logaritmica é a mais indicada.
Outra possibilidade são as transformações da família Box-Cox dada por:
yt =
{
xλt −1
λ λ 6= 0
ln(xt λ = 0
Exemplo 2.16: Considere a geração mensal de eletricidade dos Estados Unidos (em milhões
de quilowatts-hora) de todos os tipos: carvão, gás natural, nuclear, petróleo, e de vento, no
periodo entre 01/1973 - 12/2005 (Figura 2.23), em que observa-se que ao longo do tempo tem-se
um aumento da variância.
Figura 2.23: Série temporal da geração mensal de eletricidade dos Estados Unidos (em milhões
de quilowatts-hora) de todos os tipos: carvão, gás natural, nuclear, petróleo, e de vento, no
periodo entre 01/1973 - 12/2005.
Na figura 2.24 tem-se o gráfico da média versus amplitude, em que verifica-se um padrão
Análise Exploratória de Séries Temporais 37
linear entre a média e a amplitude. Assim, deve-se fazer uma transformação nos dados.
Figura 2.24: gráfico da média versus amplitude da série temporal da geração mensal de eletrici-
dade dos Estados Unidos (em milhões de quilowatts-hora) de todos os tipos: carvão, gás natural,
nuclear, petróleo, e de vento, no periodo entre 01/1973 - 12/2005.
Utilizando a transformação de boxcox com λ = −0.25 tem uma série com a variância estabi-
lizada (Figura 2.25)
Figura 2.25: Gráfico da média versus amplitude da série temporal transformada (λ = −0.25) da
geração mensal de eletricidade dos Estados Unidos (em milhões de quilowatts-hora) de todos os
tipos: carvão, gás natural, nuclear, petróleo, e de vento, no periodoentre 01/1973 - 12/2005.
3
Modelos de Suavização Exponencial
Os modelos de suavização exponencial são uma grande classe de métodos de previsão que
se baseiam na idéia de que observações passadas contem informações sobre o padrão da série
temporal. O propósito dos métodos é distinguir um padrão de comportamento de qualquer
outro ruído que possa estar contído nas observações da série e então usar esse padrão para
prever valores futuros da série.Em outras palavras, o maior peso é dado a observação mais
recente.Estes modelos tem como principais vantagens a simplicidade, baixo custo e razoável
precisão e a principal desvantagem é a dificuldade em determinar os valores mais apropriados
das constantes de suavização.
3.1 Suavização Exponencial Simples
O mais simples dos métodos de suavização exponencial é naturalmente chamado de "suavi-
zação exponencial simples"(SES). Este método é adequado para series temporais estacionarias,
ou seja, sem tendência e sazonalidade.
Seja Y1, Y2, ..., YN uma série temporal estacionária. O modelo de suavização exponencial
simples pode ser descrito como:
Ȳt = γYt + (1− γ)Ȳt−1, Ȳ0 = Y1, t = 1, . . . , N,
ou
Ȳt = γ
t−1∑
j=0
(1− γ)jYt−j + (1− γ)tȲ0, t = 1, . . . , N,
em que Ȳt é denominado valor exponencialmente suavizado e γ é a constante de suavização, tal
que 0 ≤ γ ≤ 1.
Fazendo a expansão da expressão acima temos
Ȳt = γYt + γ(1− γ)Yt−1 + γ(1− γ)2Yt−2 + . . .
assim, o modelo é uma é uma média ponderada dos valores passados em que é atribuído os pesos
maiores às observações mais recentes.
Modelos de Suavização Exponencial 39
A previsão dos h valores futuros é dada pelo último valor exponencialmente suavizado, ou
seja
Ŷt(h) = Ȳt, ∀h,
Ŷt(h) = γYt + (1− γ)Ŷt−1(h+ 1).
A última expressão pode ser interpretada como uma equação de atualização de previsão, quando
tem-se uma nova observação.
O erro de previsão at é dado por
at = Yt − Ȳt.
Supondo at ∼ N(0, σ2a), podemos construir um intervalo de confiança assintótico para Yt+h,
dado por (
Ŷt(h)− zασa
√
γ
2− γ
; Ŷt(h) + zασa
√
γ
2− γ
)
em que zα é o quantil de uma distruibuição normal padrão com nível de significância α.
A constante de suavização γ pode ser obtido pela mininização da soma de quadrado do
erro. Assim, quanto menor for o valor γ mais estaveis serão as previsões finais, uma vez que a
utilizacao de baixo valor de γ implica que pesos maiores sao dados as observações passadas e,
consequentemente, qualquer flutuação aleatória, no presente, exercera um peso menor no calculo
da previsão. Quanto mais aleatória a serie estudada, menores serão os valores de γ.
Exemplo 3.1: Aplicando a SES aos valores da série produção mensal de álcool nível 80 no
período de janeiro de 1991 a dezembro de 1995 livre de sazonalidade pelo método da regressão,
série at da Tabela 2.8, verifica-se que o valor de γ que minimiza a soma de quadrados do erro é
γ = 0, 17.
A Tabela 3.1 apresenta as doze últimas observações e os respectivos valores suavizados, e a
Figura 3.1 as séries original e suavizada.
Tabela 3.1: Valores das séries produção mensal de álcool nível 80 no ano de 1995 livre de tendência
pelo método da regressão (Yt) e suavizada pelo método SES (Ȳt).
Período Yt Ȳt
jan/95 3,65 -2,52
fev/95 -5,01 -1,47
mar/95 4,21 -2,08
abr/95 -7,36 -1,01
mai/95 -20,23 -2,09
jun/95 5,28 -5,17
jul/95 17,33 -3,39
ago/95 8,72 -5,76
set/95 12,80 -3,30
out/95 -6,96 -0,56
nov/95 10,44 -1,65
dez/95 16,59 0,40
Modelos de Suavização Exponencial 40
As previsões para os próximos meses, com origem em dezembro de 1995, são dadas por
Ŷ60(h) = Ȳ60 = 0, 40, h = 1, 2, · · · ,
que podem ser atualizadas a cada nova observação por
Ŷt(h) = 0, 17Yt + (1− 0, 17)Ŷt−1(h+ 1), h = 1, 2, · · · .
Assim, considerando o valor Y61 = 3, 13 como a observação real para o mês de janeiro de 1996,
a nova previsão para o mês de fevereiro de 1996 é
Ŷ61(1) = 0, 17Y61 + (1− 0, 17)Ŷ60(2) = 0, 17× 3, 13 + 0, 83× 0, 40 = 0, 86.
# Carregando os dados
require(xlsx)
dados=read.xlsx("alcool.xlsx",1)
attach(dados)
y=ts(alcool,start=1991,frequency=12)
par(cex.lab=1.5,cex.axis=1.5)
plot(y,xlab="Tempo",ylab="produção mensal de álcool nível 80")
# Estimando a sazonalidade
p=12
t=seq(1,60,1)
y2 <- y[1:60]
model=nls(y2~mu+a1*cos(2*pi*t/12)+a2*cos(4*pi*t/12)
+a3*cos(6*pi*t/12)+b1*sin(2*pi*t/12)+b2*sin(4*pi*t/12)
+b3*sin(6*pi*t/12),
start=list(mu=10,a1=-1,a2=10,a3=1,b1=10,b2=10,b3=1))
summary(model)
model1=nls(y2~mu+a1*cos(2*pi*t/12)+a3*cos(6*pi*t/12)
+b1*sin(2*pi*t/12)+b2*sin(4*pi*t/12)+b3*sin(6*pi*t/12)
,start=list(mu=10,a1=1,b1=10,b2=10,a3=1,b3=1))
e=resid(model1)
e=as.ts(e,start=1991,frequency=12)
par(cex.lab=1.5,cex.axis=1.5)
plot(e,xlab="Tempo",ylab=expression(a[t]))
# Suavização exponencial simples (SES)
require("stats")
ses <- HoltWinters(e, beta = FALSE, gamma = FALSE) # pacote stats
ses$alpha
ses$fitted # série suavizada
ses.prev <- predict(ses, n.ahead = 12, prediction.interval=TRUE, level = 0.95)
plot(ses, ses.prev, xlab="Tempo",ylab="", main = "")
Modelos de Suavização Exponencial 41
legend("bottomright", legend = c(
expression(bar(Y)[t]), expression(Y[t]), expression(IC[0.95])),
pch = NA,col = c("red","black","blue"),
lty = 1, bty="n", lwd = 1.5, cex = 1)
Tempo
0 10 20 30 40 50 60 70
−4
0
−2
0
0
20
40
Yt
Yt
IC0.95
Figura 3.1: Valores das séries produção mensal de álcool nível 80 no período de janeiro de 1991
a dezembro de 1995 livre de tendência pelo método da regressão (Yt) e suavizada pelo método
SES (Ȳt), e previsão para o ano de 1996, com seus respectivos intervalos de confiança.
3.2 Suavização Exponencial de Holt (SEH)
Os modelos de Suavização Exponencial de Holt (SEH) são utilizados em séries temporais que
apresentam tendência de crescimento linear. Este método envolve uma equação de previsão e
duas equações de suavização, uma para o nível e uma para a tendência:
Ȳt = γYt + (1− γ)
(
Ȳt−1 + T̄t−1
)
T̄t = β
(
Ȳt−1 + Ȳt−1
)
+ (1− β)T̄t−1
em que γ e β são constantes de suavização, tal que 0 ≤ γ ≤ 1 e 0 ≤ β ≤ 1.
A previsão dos h valores futuros é dada por:
Ŷ (h) = Ȳt + hT̄t
Os modelos de SEH necessitam de valores iniciais para as variáveis Ȳ0 e T̄0. Para o cálculo
dos valores iniciais pode-se utilizar uma regressão linear simples. As constantes γ e β do modelo
são obtidas através do cálculo da minimização da soma dos erros quadráticos de previsão.
Exemplo 3.2: Considerando a produção anual de pescado no Brasil no período de 1993 a
2012 (série Yt na Tabela 2.5) na aplicação da SEH, tem-se que os valores de γ e β que minimizam
a soma de quadrados de previsão são γ = 1 e β = 0, 29.
Modelos de Suavização Exponencial 42
A Tabela 3.2 apresenta os valores da série original (Yt), nível (Ȳt), tendência (T̄t) e suavizada
(Ŷt).
Tabela 3.2: Valores das séries produção anual de pescado no Brasil no período de 1993 a 2012
(Yt), nível (Ȳt), tendência (T̄t) e suavizada pelo método SEH (Ŷt).
t Yt Ȳt T̄t Ŷt
1993 676,40 - - -
1994 701,30 - - -
1995 652,90 701,30 24,90 726,20
1996 693,20 652,90 3,70 656,60
1997 732,30 693,20 14,29 707,49
1998 710,70 732,30 21,46 753,76
1999 744,60 710,70 9,01 719,71
2000 839,30 744,60 16,21 760,81
2001 935,90 839,30 38,91 878,21
2002 1003,30 935,90 55,59 991,49
2003 985,40 1003,30 59,01 1062,31
2004 1015,90 985,40 36,77 1022,17
2005 1008,00 1015,90 34,95 1050,85
2006 1050,80 1008,00 22,56 1030,56
2007 1072,20 1050,80 28,41 1079,21
2008 1157,60 1072,20 26,39 1098,59
2009 1241,60 1157,60 43,45 1201,05
2010 1265,50 1241,60 55,18 1296,78
2011 1433,30 1265,50 46,13 1311,63
2012 1551,20 1433,30 81,32 1514,62
As previsões para os anos subseqüentes a 2012 são apresentadas na Figura 3.2, e são calculadas
por
Ŷ2012(h) = Ȳ2012 + hT̂2012, h = 1, 2, · · · .
O intervalo de confiança de 95% de nível de confiança para as previsões e a série original
também são apresentadas na Figura 3.2.
# Carregando os dados
require(xlsx)
dados=read.xlsx("brasil.xlsx",2)
attach(dados)
y=ts(y, start = 1993, freq = 1)
par(cex.lab=1.5,cex.axis=1.5,mfrow=c(1,1))
plot(y,xlab="Tempo",ylab="ProduçãoPesqueira em Milhões de Toneladas")
# Suavização exponencial de Holt (SEH)
seh <- HoltWinters(y, gamma = FALSE) # pacote stats
seh
seh$fitted # série suavizada
prev.seh <- predict(seh, n.ahead = 5, prediction.interval=T, level = 0.95)
plot(seh, prev.seh, xlab="Tempo",
Modelos de Suavização Exponencial 43
ylab="Produção Pesqueira em Milhões de Toneladas", main="")
legend("bottomright", legend = c(
expression(bar(Y)[t]), expression(Y[t]), expression(IC[0.95])),
pch = NA,col = c("red","black","blue"),
lty = 1, bty="n", lwd = 1.5, cex = 1)
TempoPr
od
uç
ão
 P
es
qu
ei
ra
 e
m
 M
ilh
õe
s 
de
 T
on
el
ad
as
1995 2000 2005 2010 2015
10
00
15
00
20
00
Yt
Yt
IC0.95
Figura 3.2: Valores das séries produção anual de pescado no Brasil no período de 1993 a 2012
(Yt), suavizada pelo método SEH (Ŷt) e previsões com intervalo de confiança de 95% de confiança
para os anos de 2013 a 2017.
3.3 Suavização Exponencial de Holt-Winters (SEHW)
Os modelos de Holt-Winters são utilizados em séries temporais que apresentam tendência e
sazonalidade. Este método envolve três equações com três parâmetros de suavização que são
associados a cada componente da série: nível, tendência e sazonalidade.
Existem dois tipos de modelos Holt-Winters um aditivo (SEHWA) e outro multiplicativo
(SEHWM).O método aditivo é preferido quando as variações sazonais são mais ou menos cons-
tante através da série, enquanto que o método multiplicativo é preferido quando as variações
sazonais estão mudando proporcional ao nível da série. Cada um dos métodos possui três equa-
ções em relação ao padrão da série: nível, tendência e sazonalidade e, cada uma das três equações
possui uma constante de suavização diferente.
As constantes de suavização são obtidas por meio da minimização da soma de quadrado dos
erros. A dificuldade na determinação dessas constantes é uma das desvantagens do processo de
suavização Holt-Winters. Como vantagem pode ser citada a aplicabilidade do processo, adequada
a análise de séries com comportamento mais geral. Além disso, o modelo Holt-Winters não exige
a normalidade dos erros, nem que a série seja um processo estacionário (ruído branco).
Modelos de Suavização Exponencial 44
3.3.1 Suavização Exponencial de Holt-Winters Aditivo (SEHWA)
O modelo aditivo é dado pela seguinte equação:
Yt = µt + Tt + St + at
As três equações de suavização são apresentadas a seguir uma para o nível, uma para a
tendência e uma para sazonalidade.
S̄t = δ
(
Yt − Ȳt
)
+ (1− δ)S̄t−s
Ȳt = γ
(
Yt − S̄t−s
)
+ (1− γ)
(
Ȳt−1 + T̄t−1
)
T̄t = β
(
Ȳt−1 + Ȳt−1
)
+ (1− β)T̄t−1
em que γ e β são constantes de suavização, tal que 0 ≤ γ ≤ 1, 0 ≤ β ≤ 1 e 0 ≤ δ ≤ 1, e s é o
período da sazonalidade.
O valor das constantes de suavização é determinado pelo vetor (δ, γ, β) que minimiza a soma
de quadrados dos erros de previsão.
A previsão dos h valores futuros é dada por:
Ŷ (h) =
(
Ȳt + hT̄t
)
S̄t+h−s
Exemplo 3.3: Considerando os níveis de dióxido de Carbono (CO2) no norte do Canadá,
entre os anos de 1994 a 2004 (Figura 3.3), na aplicação da SEHWA, tem-se que os valores de γ,
β e δ que minimizam a soma de quadrados de previsão são γ = 0, 36, β = 0 e δ = 0, 35.
As previsões para o ano de 2005 são apresentadas na Figura 3.3.O intervalo de confiança de
95% de nível de confiança para as previsões e a série original também são apresentadas na Figura
3.3.
Detalhes da série de níveis de CO2 no norte do Canadá estão disponíveis no Exemplo 4.12.
require(TSA)
?co2
data(co2)
ts.plot(co2)
sehwa <- HoltWinters(co2, seasonal = "additive") # pacote stats
sehwa$alpha # equivalente a gamma na apostila
sehwa$beta # equivalente a beta na apostila
sehwa$gamma # equivalente a delta na apostila
sehwa$fitted # série suavizada
prev.sehwa <- predict(sehwa, n.ahead = 12, prediction.interval=T, level = 0.95)
plot(sehwa, prev.sehwa, xlab="Tempo",
ylab=expression(CO[2]), main="")
legend("bottomright", legend = c(
expression(bar(Y)[t]), expression(Y[t]), expression(IC[0.95])),
pch = NA,col = c("red","black","blue"),
Modelos de Suavização Exponencial 45
Tempo
C
O
2
1996 1998 2000 2002 2004 20063
50
36
0
37
0
38
0
Yt
Yt
IC0.95
Figura 3.3: Séries de níveis mensais de dióxido de carbono (CO2) no norte do Canadá, de janeiro
de 1994 a dezembro de 2004 (Yt), suavizada pelo método SEHWA (Ŷt) e previsões com intervalo
de confiança de 95% de confiança para os meses de 2005.
lty = 1, bty="n", lwd = 1.5, cex = 1)
3.3.2 Suavização Exponencial de Holt-Winters Multiplicativo (SEHWM)
O modelo mutiplicativo é dado pela seguinte equação:
Yt = µtSt + Tt + at
As três equações de suavização são apresentadas a seguir uma para o nível, uma para a
tendência e uma para sazonalidade.
S̄t = δ
(
Yt
Ȳt
)
+ (1− δ)S̄t−s
Ȳt = γ
(
Yt
S̄t−s
)
+ (1− γ)
(
Ȳt−1 + T̄t−1
)
T̄t = β
(
Ȳt−1 + Ȳt−1
)
+ (1− β)T̄t−1
em que γ e β são constantes de suavização, tal que 0 ≤ γ ≤ 1, 0 ≤ β ≤ 1 e 0 ≤ δ ≤ 1, e s é o
período da sazonalidade.
O valor das constantes de suavização é determinado pelo vetor (δ, γ, β) que minimiza a soma
de quadrados dos erros de previsão.
A previsão dos h valores futuros é dada por:
Ŷ (h) = Ȳt + hT̄t + S̄t+h−s
Exemplo 3.4: Considerando o total mensal (em milhares) de passageiros de voos internaci-
onais, entre os anos de 1949 a 1960 (Figura 3.4), na aplicação da SEHWM, tem-se que os valores
Modelos de Suavização Exponencial 46
de γ, β e δ que minimizam a soma de quadrados de previsão são γ = 0, 27, β = 0, 33 e δ = 0, 87.
As previsões para o ano de 1961 são apresentadas na Figura 3.4.
O intervalo de confiança de 95% de nível de confiança para as previsões e a série original
também são apresentadas na Figura 3.4.
Tempo
P
as
sa
ge
iro
s 
In
te
rn
ac
io
na
is
1950 1952 1954 1956 1958 1960 1962
10
0
30
0
50
0
70
0
Yt
Yt
IC0.95
Figura 3.4: Séries de totais mensais de passageiros internacionais, de janeiro de 1949 a dezembro
de 1960 (Yt), suavizada pelo método SEHWM (Ŷt) e previsões com intervalo de confiança de
95% de confiança para os meses de 1961.
?AirPassengers
data(AirPassengers)
ts.plot(AirPassengers)
sehwm <- HoltWinters(AirPassengers, seasonal = "multiplicative") # pacote stats
sehwm$alpha # equivalente a gamma na apostila
sehwm$beta # equivalente a beta na apostila
sehwm$gamma # equivalente a delta na apostila
sehwm$fitted # série suavizada
prev.sehwm <- predict(sehwm, n.ahead = 12, prediction.interval=T, level = 0.95)
plot(sehwm, prev.sehwm, xlab="Tempo",
ylab="Passageiros Internacionais", main="")
legend("bottomright", legend = c(
expression(bar(Y)[t]), expression(Y[t]), expression(IC[0.95])),
pch = NA,col = c("red","black","blue"),
lty = 1, bty="n", lwd = 1.5, cex = 1)
4
Modelos de Box & Jenkins
Uma metodologia bastante utilizada na análise de modelos de séries temporais é conhecida
como modelos de Box & Jenkins. Embora a maioria dos modelos utilizados por esses dois
autores já eram bastante conhecidos no segundo métade do século XX, bem antes, portanto, de
suas pesquisas. A justificativa para o uso do nome Box & Jenkins associado a esses processos
está na contribuição representativa que deram ao estudo de series temporais, através de uma
minuciosa e fundamentada integração desses processos, que lhes permitiu um tratamento analítico
de inferencia estatística nas previsões de valores futuros das variáveis dinâmicas.
Os modelos Box & Jenkins podem ser classificados em estacionários e não estacionários
4.1 Modelos para Séries Estacionárias
Modelos estacionários são aqueles que assumem que o processo está em equilíbrio, ou seja que
a série temporal desenvolve-se no tempo com médias e variância constantes, sem a necessidade
de uma transformação para estabiliza-la.
Suponha que exista um sistema (Figura 4.1) cuja entrada é um ruído branco at, no qual
possui um filtro que que gera uma sequência de valores observados seguindo um padrão, que
corresponde à série temporal Yt.
Figura 4.1: Representação de uma série temporal como saída de um filtro linear
Formalmente,tem-se:
Yt = µ+ at + ψ1at−1 + ψ2at−2 + ...
= µ+ (1 + ψ1B + ψ2B
2 + ...)at
= µ+ ψ(B)at
em que ψ(B) é a função de transferência do filtro e µ é um parâmetro que determina o nível de
série.
Modelos de Box & Jenkins 48
Fazendo:
Ỹt = Yt − µ
Assim, tem-se
Ỹt = ψ(B)at =
∞∑
j=0
ψjat−j
em que com ψ0 = 1
Desta forma o processo linear Ỹt, é escrito como uma soma ponderada de ruídos branco do
presente e do passado. Esse processo consiste em uma sequência de variáveis aleatórias não-
correlacionadas com média zero e variância constante, isto é,
E[at] = 0 ∀ t
V ar[at] = σ
2 <∞ ∀ t
E[asat] = 0 ∀ s 6= t
Se a série
∞∑
j=0
ψjat−j for convergente então Yt é estacionária e µ é a média do processo. Caso
contrário Yt é não estacionária e µ não tem significado especifico.
Podemos escrever Yt é uma forma alternativa, como uma soma ponderada de valores passados
de Ỹt mais um ruído branco at:
Ỹt = π1Ỹt−1 + π2Ỹt−2 + ...+ at
=
∞∑
j=1
πj Ỹt−j + at
Ỹt −
∞∑
j=1
πj Ỹt−j = at1− ∞∑
j=1
πjB
j
 Ỹt = at
π(B)Ỹt = at
em que π(B) é o operador dador por π(B) = 1− π1B − π1B2 − .... Temos que
Ỹt = ψ(B)at
Assim,
π(B)Ỹt = at
π(B)ψ(B)at = at
Como at = at, temos que
π(B)ψ(B) = 1 π(B) = ψ−1(B)
Modelos de Box & Jenkins 49
Um processo liner será:
• Estacionário se a série ψ(B) convergir para |B| ≤ 1
• Invertível se a série π(B) convergir para |B| ≤ 1
4.1.1 Modelos Auto-Regressivos (AR)
Os modelos auto-regressivos se baseiam na idéia de que a observação presente da série Yt
pode ser explicada como uma função de p observações passadas, Zt−1, Zt−2, ..., Zt−p, em que p,
denominado ordem do modelo, determina o número de passos entre as observações passadas e a
observação presente
Ummodelo auto-regressivo de ordem p, denotado por AR(p), é descrito pela expressão abaixo:
Yt = µ+ φ1Yt−1 + φ2Yt−2 + ...+ φpYt−p + at
Fazendo
Ỹt = Yt − µ
Pode-se escrever o modelo em termos da operador diferença:
Ỹt = φ1Ỹt−1 + φ2Ỹt−2 + ...+ φpỸt−p + at
Ỹt − φ1Ỹt−1 − φ2Ỹt−2 − ...− φpỸt−p = at
(1− φ1B − φ2B2 − ...− φpBp)Ỹt = at
φ(B)Ỹt = at
Assim, temos que φ(B) = π(B).
O modelo autoregressivo de ordem 1 e 2 pode ser representado pelas expressões
AR(1) ⇒ Ỹt = φ1Ỹt−1 + at
AR(2) ⇒ Ỹt = φ1Ỹt−1 + φ2Ỹt−2 + at
Como o polinômio φ(B) = π(B) a condição de estacionaridade do AR(p) estabelece que
todas as raízes da equação devem cair fora do círculo unitário, ou seja |B| > 1. E, sendo um
processo de ordem finita, então temos que AR(p) é sempre invertível.
Exemplo 4.1: Seja um modelo auto-regressivo de ordem 2
Yt = 0, 5Yt−1 − 0, 3Yt−2 + et
O modelo é invertível, vamos verificar se é estacionário
φ(B) = 1− 0, 5B + 0, 3B2
1− 0, 5B + 0, 3B2 = 0
B1 = −2, 83 e B2 = 1, 77
Modelos de Box & Jenkins 50
Como |B1| > 1 e |B2| > 1 então o modelo é estacionário.
Exemplo 4.2: Seja um modelo auto-regressivo de ordem 2
Yt = 1, 61Yt−1 − 0, 56Yt−2 + et
O modelo é invertível, vamos verificar se é estacionário
φ(B) = 1− 1, 61B + 0, 56B2
1− 1, 61B + 0, 56B2 = 0
B1 = 0, 91 e B2 = 1, 97
Como |B1| < 1 e |B2| > 1 então o modelo não é estacionário.
Pela condição de estacionaridade as autocovariâncias nao dependem do t, e sim da ordem p.
Assim quanto maior o valor de p, ou seja, quanto maior a distancia entre as observações, menor
a autocovariância.
A identificação do modelo é feita por meio das funções de autocorrelação (FAC) e autocor-
relação parcial (FACP). A (FAC) decai exponencialmente, alternando ou não de sinal, em geral
é uma mistura de exponenciais e ondas senóides amortecidas e a autocorrelação parcial (FACP)
apresenta um corte rápido no lag significativo, indicando a ordem p do modelo.
Exemplo 4.3: Na figura 4.2 temos de uma série temporal com suas funções de autocorrelação
e autocorrelação parcial. Verifica-se que a FAC decai rapidamente para zero, a partir do lag 2 é
não significativo, e a FACP tem apenas o primeiro lag significativo. Assim temos a caracterização
de um processo auto-regressivo, e como apenas o primeiro lag é significativo na FACP temos um
processo AR(1).
Figura 4.2: Serie temporal, funções de autocorrelação e autocorrelação parcial de um processo
AR(1)
Exemplo 4.4: Na figura 4.3 temos de uma série temporal com suas funções de autocorrelação
e autocorrelação parcial. Verifica-se que a FAC decai exponencialmente e partir do lag 10 é não
Modelos de Box & Jenkins 51
significativo, e a FACP tem apenas lag 1 e 2 significativo. Assim temos a caracterização de um
processo auto-regressivo, como o lag 1 e 2 são significativos na FACP temos um processo AR(2).
Figura 4.3: Serie temporal, funções de autocorrelação e autocorrelação parcial de um processo
AR(2)
Exemplo 4.5 (Ozônio em Azuza, Califórnia (EUA) (Morettin e Toloi, 2006)): Vamos co-
niderar a série de valores mensais concentração de Ozônio (O3) na cidade Azuza do estado da
Califórnia (EUA), de janeiro de 1956 a dezembro de 1970, apresentada na Figura 4.4.
Período
O
3
1960 1965 1970
2
4
6
8
Figura 4.4: Concentração mensal de Ozônio na cidade Azuza, Califórnia (EUA), de janeiro de
1956 a dezembro de 1970.
Esta série apresenta evidências de comportamento sazonal (Figura 4.5), que é confirmada
pelo teste de Fisher para o período 12 (valor − p < 0, 01).
Modelos de Box & Jenkins 52
Período
O
3
1968.0 1968.5 1969.0 1969.5 1970.0 1970.5 1971.0
2
3
4
5
6
7
8
9
J
F
M
A
M
J
J
A
S
O
N
D J
F
M
A M
J
J
A
S
O
N
D
J
F
M
A
M
J
J
A
S
O
N
D
Figura 4.5: Concentração mensal de Ozônio na cidade Azuza, Califórnia (EUA), de janeiro de
1968 a dezembro de 1970.
Considerando a sazonalidade como determinística, ajusta-se um modelo de regressão não
linear
St = µ+
2∑
i=1
[
αicos
(
2πit
12
)
+ βisen
(
2πit
12
)]
+ �t,
com �t a componente aleatória do modelo. As estimativas dos parâmetros, com seus respectivos
erros padrões, valores do teste t e valor − p são dadas na Tabela 4.1.
Tabela 4.1: Estimativas dos parâmetros do modelo de regressão não linear, com seus respectivos
erros padrões, valores do teste t e valor − p.
Parâmetro Estimativa Erro padrão tc valor − p
µ 5,08 0,0725 70,07 < 0, 01
α1 -2,01 0,1026 -19,59 < 0, 01
α2 -0,26 0,1026 -2,53 0, 01
β1 -1,55 0,1026 -15,11 < 0, 01
β2 0,52 0,1026 5,07 < 0, 01
A Figura 4.6 apresenta a série livre de sazonalidade determinística. Pelo teste ADF a série
sazonalmente ajusta é estacionária (valor − p < 0, 01). Pelo comportamento das FAC e FACP,
Figuras 4.7 (A) e (B) respectivamente, recomenda-se preliminarmente um modelo AR(1).
O modelo ajustado para a concentração mensal de Ozônio é
Yt = 5, 08− 2, 01cos
(
2πt
12
)
− 0, 26cos
(
4πt
12
)
− 1, 55sen
(
2πt
12
)
+ 0, 52sen
(
4πt
12
)
+ 0, 15Yt−1 + at,
em que at é um ruído branco, com V̂ ar(at) = 0, 90.
Modelos de Box & Jenkins 53
Período
O 3
−Ŝ
1960 1965 1970
−3
−2
−1
0
1
2
3
Figura 4.6: Série concentração mensal de Ozônio na cidade Azuza, Califórnia (EUA), de janeiro
de 1968 a dezembro de 1970, sazonalmente ajustada.
(A) (B)
0.5 1.0 1.5
−
0
.1
5
−
0
.1
0
−
0
.0
5
0
.0
0
0
.0
5
0
.1
0
0
.1
5
Lag
A
C
F
0.5 1.0 1.5
−
0
.1
5
−
0
.1
0
−
0
.0
5
0
.0
0
0
.0
5
0
.1
0
0
.1
5
Lag
P
a
rt
ia
l 
A
C
F
Figura 4.7: Função de autocorrelação (A) e autocorrelação parcial (B) da série temporal de
concentração mensal de Ozônio na cidade Azuza, Califórnia (EUA), sazonalmente ajustada.
require(xlsx)
require(TSA)
dados=read.xlsx("OZONIO.xlsx",1)
attach(dados)
O3 <- ts(Ozonio, start = 1956, freq = 12)
plot(O3, xlab = "Período", ylab = expression(O[3]))
plot(window(O3,start=c(1968,1)), xlab = "Período",
ylab = expression(O[3]))
Meses = c("J","F","M","A","M","J","J","A","S","O","N","D")
Modelos de Box & Jenkins 54
points(window(O3,start=c(1968,1)),pch=Meses)
P=periodograma(O3)
Fisher.test(P)
g zalpha valorp periodo
0.7464016 0.0807706 8.380618e-52 12
per <- 12
t <- 1:length(O3)
model=nls(O3~mu+a1*cos(2*pi*t/per)+a2*cos(4*pi*t/per)
+b1*sin(2*pi*t/per)+b2*sin(4*pi*t/per)
,start=list(mu=5,a1=-1,a2=-2,b1=2,b2=1))
summary(model)
Formula: O3 ~ mu + a1 *cos(2 * pi * t/per) + a2 * cos(4 * pi * t/per) +
b1 * sin(2 * pi * t/per) + b2 * sin(4 * pi * t/per)
Parameters:
Estimate Std. Error t value Pr(>|t|)
mu 5.07944 0.07253 70.030 < 2e-16 ***
a1 -2.00839 0.10258 -19.579 < 2e-16 ***
a2 -0.26333 0.10258 -2.567 0.0111 *
b1 -1.54999 0.10258 -15.110 < 2e-16 ***
b2 0.52539 0.10258 5.122 7.93e-07 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 0.9731 on 175 degrees of freedom
Number of iterations to convergence: 1
Achieved convergence tolerance: 7.814e-09
O3.SA <- ts(residuals(model), start = 1956, freq = 12)
plot(O3.SA, xlab = "Período", ylab = expression(O[3] - hat(S)),
type = "l")
require(fUnitRoots)
lags <- trunc((length(O3.SA)-1)^(1/3))
adfTest(O3.SA, lags = lags) # Teste de Dickey-Fuller Aumentado (ADF)
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 5
STATISTIC:
Dickey-Fuller: -5.1052
P VALUE:
Modelos de Box & Jenkins 55
0.01
acf(O3.SA, main = "")
pacf(O3.SA, main = "")
m1 <- arima(O3.SA, order = c(1,0,0))
m1
Call:
arima(x = O3.SA, order = c(1, 0, 0))
Coefficients:
ar1 intercept
0.1548 -0.0011
s.e. 0.0736 0.0835
sigma^2 estimated as 0.8985: log likelihood = -245.79, aic = 495.57
4.1.2 Modelos de Médias Móveis (MA)
Os modelos de médias móveis se baseiam na ideis de que a observação presente da série Yt
pode ser explicado como uma função do ruído branco at, ocorrido no período corrente e nos
períodos passados.
Um modelo de média móveis de ordem q denotado por MA(q), é descrito pela expressão
abaixo:
Yt = µ+ at − θ1at−1 − θ2at−2 − ...− θqat−q
Fazendo
Ỹt = Yt − µ
Pode-se escrever o modelo em termos da operador diferença:
Ỹt = at − θ1at−1 − θ2at−2 − ...− θqat−q
Ỹt = (1− θ1B − θ2B2 − ...− θqBq)at = at
Ỹt = θ(B)at
Assim, temos que θ(B) = ψ(B).
O modelo de médias móveis de ordem 1 e 2 pode ser representado pelas expressões
MA(1) ⇒ Ỹt = at − θ1at−1
MA(2) ⇒ Ỹt = at − θ1at−1 − θ2at−2
Como o polinômio θ(B) = ψ(B) e sendo um processo de ordem finita, então temos queMA(q)
é sempre invertível. O processo MA(q) será invertível se todas as raízes da equação cairem fora
do círculo unitário, ou seja |B| > 1.
Modelos de Box & Jenkins 56
Exemplo 4.6: Seja um modelo de média móveis de ordem 1
Yt = at + 0, 8at−1
O modelo é estacionário, vamos verificar se é invertível
θ(B) = 1 + 0, 8B
1 + 0, 8B = 0
B = −1, 25
Como |B| > 1 então o modelo é invertível.
A identificação do modelo MA é feita por meio das funções de autocorrelação (FAC) e
autocorrelação parcial (FACP). A FACP decai exponencialmente se todos os parâmetros forem
positivo, caso contrário forma uma senóide amortecida. A FAC apresenta um corte rápido no
lag significativo, indicando a ordem q do modelo.
Exemplo 4.7: Na figura 4.8 temos de uma série temporal com suas funções de autocorrelação
e autocorrelação parcial. Verifica-se que a FACP forma uma senóide amortecida, e a FAC tem
apenas lag 1 e 2 significativo. Assim temos a caracterização de um processo de média móveis,
como o lag 1 e 2 são significativos na FACP temos um processo MA(2).
Figura 4.8: Serie temporal, funções de autocorrelação e autocorrelação parcial de um processo
MA(2)
Exemplo 4.8 (Preço de petróleo (Cryer e Chan, 2008)): Considere a série temporal de
preço mensal do barril de petróleo bruto (em dólares) na cidade de Cushing, Oklahoma (EUA),
de janeiro de 1986 a dezembro de 2005, representada na Figura 4.9.
Visualmente, a série apresenta uma volatilidade para o preço do petróleo, ou seja, a variância
não é constante ao longo do tempo. Sendo assim, uma transformação logarítmica pode ser
adequada. O gráfico de amplitudes versus médias, Figura 4.10, confirma essa hipótese.
Modelos de Box & Jenkins 57
Ano
Pr
eç
o 
(U
S$
)
1990 1995 2000 2005
10
20
30
40
50
60
Figura 4.9: Preço mensal do barril de petróleo bruto (em dólares) na cidade de Cushing,
Oklahoma (EUA), de janeiro de 1986 a dezembro de 2005.
20 30 40 50
5
1
0
1
5
2
0
Médias
A
m
p
li
tu
d
e
s
Figura 4.10: Relação entre médias e amplitudes dos preços mensais do barril de petróleo bruto,
para grupos de 12 observações.
A Figura 4.11 apresenta a série transformada, na qual nota-se a presença da componente de
tendência, confirmada pelo teste de Cox-Stuart (valor − p < 0, 01).
A série livre de tendência, isto é, resultante da diferença de uma defasagem (lag = 1) da
transformação logarítmica, é apresentada na Figura 4.12. Pelo teste de Fisher rejeita-se a hipótese
de presença da componente sazonal (valor − p = 0, 57), para o período 2, 20, e pelo teste ADF
(valor − p < 0, 01) a série é considerada estacionária.
A função de autocorrelação (FAC) e autocorrelação parcial (FACP), Figuras 4.13 (A) e (B)
respectivamente, caracterizam um processo MA(1). Pois, a FACP forma uma onda senóide
amortecidas e a FAC é significativa apenas para a primeira defasagem (lag = 1).
Modelos de Box & Jenkins 58
Ano
lo
g(
P
re
ço
)
1990 1995 2000 2005
2.
5
3.
0
3.
5
4.
0
Figura 4.11: Logarítmo do preço mensal do barril de petróleo bruto (em dólares) na cidade de
Cushing, Oklahoma (EUA), de janeiro de 1986 a dezembro de 2005.
Ano
(1
−B
) lo
g(
Pr
eç
o)
1990 1995 2000 2005
−0
.4
−0
.2
0.
0
0.
2
0.
4
Figura 4.12: Série diferença de defasagem 1 (lag = 1) do logarítmo do preço mensal do barril
de petróleo bruto (em dólares) na cidade de Cushing, Oklahoma (EUA), de janeiro de 1986 a
dezembro de 2005, livre de tendência por diferenciação de uma defasagem
Modelos de Box & Jenkins 59
(A) (B)
0.5 1.0 1.5
−
0
.1
0
.0
0
.1
0
.2
Lag
A
C
F
0.5 1.0 1.5
−
0
.2
−
0
.1
0
.0
0
.1
0
.2
Lag
P
a
rt
ia
l 
A
C
F
Figura 4.13: Função de autocorrelação (A) e autocorrelação parcial (B) da série diferença de
defasagem 1 do logarítmo do preço mensal do barril de petróleo bruto.
O modelo MA(1) ajustado é dado por
Yt = at − 0, 29at−1,
em que at é um ruído branco, com V̂ ar(at) = 0, 006. O modelo é invertível pois,
θ(B) = 1− 0, 29B
1− 0, 29B = 0
B = 3, 45.
library("TSA") # carregar pacote "TSA"
data(oil.price) # carregar dados "oil.price"
# Considerando o período de 1/1986 a 12/2005:
op <- ts(oil.price[1:240], start = 1986, freq = 12)
plot(op, xlab = "Ano", ylab = "Preço (US)")
# -- Verificando se precisa de transformação
M.op <- matrix(op,nrow=12)
medias <- apply(M.op,2,mean)
maximos <- apply(M.op,2,max)
minimos <- apply(M.op,2,min)
amplitude <- maximos-minimos
plot(medias,amplitude,xlab="Médias", ylab="Amplitudes")
abline(lm(amplitude ~ medias))
# --
log.op <- log(op)
Modelos de Box & Jenkins 60
plot(log.op, xlab = "Ano", ylab = "log(Preço)")
cs.test(log.op)
Exact Cox-Stuart trend test
data: log.op
S = 24, p-value = 2.162e-11
alternative hypothesis: data have a monotonic trend
d.log.op <- diff(log.op)
plot(d.log.op, xlab = "Ano", ylab = "(1-B) log(Preço)")
P=periodograma(d.log.op)
Fisher.test(P[1:12,])
g zalpha valorp periodo
0.2418051 0.3924008 0.5711773 2.201835
require(fUnitRoots)
lags <- trunc((length(d.log.op)-1)^(1/3))
adfTest(d.log.op, lags = lags) # Teste de Dickey-Fuller Aumentado (ADF)
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 6
STATISTIC:
Dickey-Fuller: -6.4764
P VALUE:
0.01
acf(d.log.op, main = "")
pacf(d.log.op, main = "")
arima(d.log.op, order = c(0,0,1))
Call:
arima(x = d.log.op, order = c(0, 0, 1))
Coefficients:
ma1 intercept
0.2913 0.0036
s.e. 0.0696 0.0068
sigma^2 estimated as 0.006677: log likelihood = 259.42, aic = -514.83
Modelos de Box & Jenkins 61
4.1.3 Modelos Auto-Regressivos de Médias Móveis (ARMA)
Os modelos auto-regressivo de média móveis são uma combinação dos dois anteriores onde
Yt é descrito por seus valores passados e pelos ruídos branco corrente e passados.
Um modelo auto-regressivo de média móveis de ordem p, q denotado por ARMA(p, q), é
descrito pela expressão abaixo:
Yt = φ1Yt−1 + φ2Yt−2 + ...+ φpYt−p + at − θ1at−1 − θ2at−2 − ...− θqat−q
Pode-se escrever o modelo em termos da operador diferença:
Yt = φ1Yt−1