1
Mestrado em Análise de Dados e Sistemas de Apoio à Decisão
Disciplina: Extracção e conhecimento de dados I
Trabalho nº4: Séries Temporais
Data: 26 de Fevereiro de 2005
Aluno: Elisabeth Silva Fernandes
Nº 050414012
2
Índice
Introdução e objectivo do trabalho………………………………………………………..3
Análise preliminar dos dados…………………………………………………………..4
Modelos para séries não estacionárias………………………………………………..12
Modelo não-
lineares…………………………………………………………………...14
Redes Neuronais……………………………………………………………………….19
Projection pursuit
regression………………………………………………………….21
Modelo de
Kernel………………………………………………………………………22
Conclusão………………………………………………………………………………....2
5
3
Introdução e objectivo do trabalho
No âmbito da disciplina de Extracção e Conhecimento de Dados foi proposta a realização
de um trabalho cujo objectivo é o de tratar uma série temporal de um activo financeiro,
utilizando métodos fornecidos nas aulas, e prever os valores de fecho do índice S&P
500, nos dias 27, 28 de Fevereiro e 1, 2 e 3 de Março.
4
Análise preliminar dos dados
No gráfico seguinte é apresentada a evolução da séries temporal do valor do fecho do
S&P 500 desde 3 de Janeiro de 1950. Pode observar-se mudanças repentinas do valor
de fecho deste índice, pode ainda constatar-se que a variância evolui ao longo do
tempo (volatilidade) e que existe irreversibilidade no tempo.
0 2000 4000 6000 8000 10000 12000 14000
05
00
10
00
15
00
Index
da
do
s
Figura 1
Uma vez que esta série exibe comportamentos incompatíveis com a formulação de
um processo linear apresento na secção seguinte modelos não-lineares.
5
Um primeiro problema que se coloca é a escolha dos dados a usar para obtenção dos
modelos. O que a literatura indica é que deve usar-se um conjunto de dados que
tenham um comportamento parecido com a totalidade da série, logo decidi utilizar
aproximadamente os últimos 1000 registos (que corresponde, aproximadamente, à
zona dentro do círculo azul). Parece-me que estes dados têm um comportamento
parecido com o da série completa.
Os valores de fecho deste índice que vão ser considerados são de 6 de Março de 2002
até 23 de Fevereiro de 2006.
data<-read.table('C:/Documents and Settings/Utilizador/Ambiente de
trabalho/trabalho4-series/treino.txt',header=TRUE)
É possível observar o comportamento desta série na figura seguinte:
0 200 400 600 800 1000
80
09
00
10
00
11
00
12
00
13
00
Index
fe
ch
o
Figura 2
Mais uma vez se verificam alterações repentinas nos valores, volatilidade e
irreversibilidade no tempo. Verifica-se ainda alguma tendência.
Algumas estatísticas descritivas dos valores de fecho deste índice:
6
Min. 1st Qu. Median Mean 3rd Qu. Max.
776.8 983.5 1111.0 1079.0 1185.0 1294.0
O método das diferenças é uma operação simples que transforma uma série não-
estacionária numa série estacionária.
> h.returns<-function(x,h=1){
diff(x,lag=h)/x[1:(length(x)-h)]
}
> h.returns(data[,’Close’])
> plot(h.returns(data[,’Close’],h=1))
0 200 400 600 800 1000
-0
.0
4-0
.0
20
.0
00
.0
20
.0
40
.0
6
Index
h.re
tu
rn
s(d
ad
os, h
=
1
)
Figura 3- Série estacionário
A figura 4 revela uma distribuição simétrica dos valores, com clara presença de
eventos raros o que é de esperar visto que estamos a tratar uma série económica.
7
-0.04
-0.02
0.00
0.02
0.04
0.06
Figura 4
Uma outra análise que pode ser feita é comparar o histograma da nova série com a
curva gaussiana.
hist(h.returns(data[,’Close’],h=1),prob=T)
lines(density(h.returns(data[,’Close’],h=1)))
Histogram of h.returns(dados, h = 1)
h.returns(dados, h = 1)
De
nsity
-0.04 -0.02 0.00 0.02 0.04 0.06
01
02
03
04
0
Figura 5
Se estimarmos os parâmetros da normal e se sobrepuser a curva ao gráfico anterior,
fazendo:
8
x<-seq(-0.06,0.06,length=100)
y<-
dnorm(x,mean(h.returns(data[,’Close’],h=1)),sd(h.returns(data[,’Clo
se’],h=1)))
lines(x,y,col="red")
Obtenho o seguinte gráfico:
Histogram of h.returns(dados, h = 1)
h.returns(dados, h = 1)
De
nsity
-0.04 -0.02 0.00 0.02 0.04 0.06
01
02
03
04
0
Verifica-se que o histograma se aproxima muito da distribuição normal.
A função de Autocorrelação fornece-nos uma forma de ver o quanto os valores da
série temporal estão correlacionados com os valores anteriores.
> FECHO<-h.returns(data[,’Close’],h=1)
> acf(FECHO)
9
0 5 10 15 20 25 30
0.0
0.2
0.4
0.6
0.8
1.0
Lag
AC
F
Series FECHO
Figura 6: h=1
Pela análise deste gráfico podemos concluir que existe alguns valores com correlação
significativa.
0 5 10 15 20 25 30
0.0
0.2
0.4
0.6
0.8
1.0
Lag
AC
F
Series FECHO
Figura 7: h=5
Autocorrelações parciais
pacf(fecho)
10
0 5 10 15 20 25 30
0.0
0.2
0.4
0.6
0.8
1.0
Lag
Pa
rtia
l A
CF
Series fecho
Filtro Linear – Médias móveis
O termo média móvel é utilizado porque à medida que a próxima observação se torna
disponível, a média das observações é recalculada, incluindo essa observação no conjunto
de observações e desprezando a observação mais antiga.
plot(data[,’Close’])
fecho.1<-filter(data[,’Close’],filter=rep(1/30,30))
lines(fecho.1,col="red")
fecho.2<-filter(data[,’Close’],filter=rep(1/90,90))
lines(fecho.2,col="blue")
11
0 2 0 0 4 0 0 6 0 0 8 0 0 1 0 0 0
80
09
00
10
00
11
00
12
00
13
00
Ind e xfe
ch
o
Figura 8
0 200 400 600 800 1000
80
09
00
10
00
11
00
12
00
13
00
Index
fe
ch
o
Figura 9
O filtro que melhor se ajusta aos dados é o filtro 1.
12
Modelos para séries não estacionárias
Nesta secção apresento um modelo obtido para esta série temporal usando modelos
ARIMA(p,d,q).
Uma vez que este tipo de séries não é estacionária faz-se a transformação log à série
e esta passa a ser estacionária. Após esta tranformação já se pode utilizar um modelo
ARIMA(p,d,q).
data<-read.table('C:/Documents and Settings/Utilizador/Ambiente de
trabalho/trabalho4-series/treino.txt',header=TRUE)
> SPteste<-data[999:1002,5]
O modelo que obtive é um ARIMA(1,0,2) e os passos segui para calcular o erro
quadrático médio das previsões da semana de 22 de Fevereiro até 24 de Fevereiro
foram os seguintes:
ARIMA<-arima0(log(data[,5]), order = c(1,0,2))
tsdiag(ARIMA, gof.lag=20)
Pelo teste de Ljung-Box-Pierce posso afirmar que este modelo se adequa à série em
estudo.
13
Standardized Residuals
Time
0 200 400 600 800 1000
-4
02
40 5 10 15 20 25 30
0.0
0.4
0.8
Lag
AC
F
ACF of Residuals
5 10 15 20
0.0
0.4
0.8
p values for Ljung-Box statistic
lag
p value
predict(arima0(log(SPtreino[,5]), order = c(1,0,2)), n.ahead = 5)
A previsão para a semana de teste é a seguinte:
pred<-exp(c(7.159758, 7.159168 ,7.158585, 7.158003))
> pred
[1] 1286.600 1285.841 1285.091 1284.344
sum((pred-SPteste[,5])^2)/4
[1] 23.13420
14
Modelo não-lineares
Modelo ARCH- Autocoregressive Conditionally heteroscedastic
O primeiro modelo que fornece uma forma de modelar a volatilidade é o modelo
ARCH proposto por Engle em 1982.
Este modelo assume que a variância dependo dos quadrados dos erros passados, o
modelo é da forma:
ttty εσ=
2
110
2
−+=
ttyαασ
onde t
ε é ruído branco e 01
>α .
O modelo mais simples é o ARCH(1) que é um ruido branco com variância
condicional não constante, e essa variância condicional depende do valor anterior.
Neste caso, o modelo pode ser reescrito como um AR(1) estacionário para 2
ty com
erros não normais que têm média zero e variância não constante.
15
Se o processo ARCH(1) for estacionário pode mostrar-se que o coeficiente de curtose
é dado por
( )
2
2
31
13
α
ε
−
−
=k
que é sempre maior do que 3 ( a curtode da distribuição normal). Os processos
ARCH(1) têm caudas mais pesadas do que a distribuição normal e são portanto
adequados para modelar séries temporais com esta característica.
O modelo ARCH(m) é dado por:
22
110
2
...mtmtt
yy−−
+++= ααασ
Identificação
A variância condicional dos erros comporta-se como um processo autorregressivo,
logo deve esperar-se que os resíduos de um modelo ARMA ajustado a uma série
temporal observada também sigam este padrão característico. Se o modelo ajustado
for adequado então a FAC e a FACP dos resíduos devem indicar um processo
puramente aleatório, no entanto se a FAC dos quadrados dos resíduos tiver um
decaimento característico de uma autorregressão isto é uma indicação de que o
modelo ARCH pode ser apropriado. A ordem m do modelo pode ser identificada
através da FACP dos quadrados dos resíduos.
Modelos GARCH
Uma generalização dos modelos ARCH consiste em assumir que a variância
condicional se comporta como um processo ARMA, isto é, depende também dos seus
valores passados.
O modelo GARCH(m,r) é da forma:
ttty εσ=
∑ ∑= =
−−++=
m
j
r
j
jtjjtjty
1 1
22
0
2
σβαασ
16
Aplicação na série temporal de fecho do S&P500
O conjunto de dados utilizado para obter este modelo são os valores de 6 de
Março de 2002 até 17 de Fevereiro de 2006.
> data<-read.table('C:/Documents and Settings/Utilizador/Ambiente de
trabalho/trabalho4-series/treino.txt',header=TRUE)
> SPtreino1<- data[1:998,5]
> SPteste<-data[999:1002,5]
Pelas consultas feitas em bibliografia da área conclui que um modelo que se
adequa a esta série temporal é o seguinte:
GARCH(1,1)
2
11
2
110
2
−−++=
ttty σβαασ
G<-garch(x=diff(SPtreino1),c(1,1))
plot(G)
O modelo obtido foi o seguinte:
17
Coefficient(s):
a0 a1 b1
0.64616 0.03839 0.95440
Aplicando o Garch à série directamente obtenho o seguinte modelo:
Coefficient(s):
a0 a1 b1
1.477e+04 9.851e-01 2.095e-03
Valores previstos do dia 21 de Fevereiro até 24 de Fevereiro.
1284.705
1280.547
1290.053
1285.256
mse<-(sum((prev-SPteste[,5])^2)/4)
MSE=43.07905
Figura 10 – Gráfico da série suavizada e dos
resíduos.
Figura 11- Histograma da série e dos resíduos.
18
Figura 12 – Gráfico dos quartis.Figura 13 – FAC da série e dos resíduos
Pelos gráficos anteriores é possível verificar que este método modela razoavelmente
a série em estudo, uma vez que não é rejeitada a normalidade dos resíduos e a FAC
dos resíduos a partir do primeiro lag tem todos os valores dentro das bandas de
confiança.
Efectuei outras experiências mas pareceu-me ser este o melhor modelo de entre este
tipo de modelos. De seguida apresento um dos modelos que também me pareceu
razoável: ARCH(2)
A<-garch(diff(SPtreino1), order = c(0,2))
Como é possível observar nos gráficos seguintes os Q-Qplot aproximam-se de
uma recta e o histograma dos resíduos parece aproximar-se de uma Normal, logo este
modelo ajusta-se razoavelmente bem aos dados.
19
Figura 14
Figura 15
Figura 16
Figura 17
Previsões para os dias 21 de Fevereiro até 24 de Fevereiro:
1284.807
1280.736
1289.650
1285.474
MSE= 41.17191
O modelo que tem maior MSE é o Garch(1,1), logo de entre estes escolheria o
Arch(2) para prever valores futuros.
Redes Neuronais
20
Uma vez que as séries económicas são (geralmente) não lineares e como as
redes neuronais lidam facilmente com este tipo de casos, a aplicação de Redes
Neuronais para prever valores futuros deste tipo de séries tem sido muito usada.
O conjunto de dados utilizado para obter a rede neuronal, SPtreino , são os
valores de 6 de Março de 2002 até 17 de Fevereiro de 2006.
Os valores utilizados para teste, SPteste, são os valores de 20 de Fevereiro até 24 de
Fevereiro.
SPtreino<-data[1:998,]
SPteste<-data[999:1002,]
O objectivo neste método é obter o melhor modelo de Rede Neuronal, para tal
utilizamos as primeiras 250 observações para obter as diferentes variantes do
modelo, e as restantes são usadas para decidir qual é o “melhor” modelo. Finalmente,
este “melhor” modelo será novamente obtido mas com todos os dados (SPtreino).
train<-data[1:250,]
sel<-data[251:500,]
test<-data[751:1002,]
train.best<-data
library(nnet)
alt<-expand.grid(size=c(10,20),decay=c(0.01,0.001),theil=0)
prevs.ant<-c(train[nrow(train),"Close"],sel[1:(nrow(sel)-1),"Close"])
for(a in 1:nrow(alt)){
nn<-
nnet(Close~.,train[,-1],size=alt[a,"size"],decay=alt[a,"dec
ay"],maxit=1000,linout=T)
prevs<-predict(nn,sel)
alt[a,"theil"]<-sqrt(sum(((sel[,"Close"]-
prevs)/prevs.ant)^2)/sum(((sel[,"Close"]-prevs.ant)/prevs.ant)^2))
}
best<-which.min(alt[,"theil"])
21
A melhor rede neuronal é a que tiver menor valor da estatística de Theil.
nn<-
nnet(Close~.,data=train.best[,-1],size=alt[best,"size"],decay=alt
[best,"decay"],maxit=1000,linout=T)
prevs<-predict(nn,test)
Previsões para os valores da semana de 21 de Fevereiro até 24 de Fevereiro.
1285.3627
1290.4732
1287.4952
1288.3014
MSE= 2.907016
22
Projection pursuit regression
O projection Pursuit é um tipo de modelo aditivo com a forma:
= ∑∑==
a
j
iij
F
i
iXfxr
1
,
1
)( α
Isto é, trata-se da soma de funções de combinações lineares dos atributos.
Os passos seguidos para a obtenção deste modelo são os seguintes:
library(modreg)
pp<-ppr(V5~.,data=SPtreino[,-1],nterms=5)
pp.preds<-predict(pp,SPteste)
As previsões obtidas para a semana de teste são:
pp.preds
1288.697
1294.194
1289.575
1289.561
mse<-sum((SPteste[,5]-pp.preds)^2)/4
> mse
[1] 9.410206
23
Modelo de Kernel
Um outro modelo que decide experimentar é o modelo de Regressão de Kernel, este
modelo faz uma média pesada à volta da vizinhança de cada caso de teste.
Uma função implementada no R para este método é kernapply(), que usei, para ver
o que se obtêm apresento os seguintes gráficos e experiências.
k1 <- kernel("daniell", 15)
x1 <- kernapply(data[,5], k1)
plot(data[,5])
lines(x1, col = "red")
0 200 400 600 800 1000
80
09
00
10
00
11
00
12
00
13
00
Index
da
ta
[, 5
]
Figura 18
Usando o seguinte valor kd para o parâmetro da dimensão de kernel, o modelo obtido
ajusta-se melhor à série em estudo.
kd <- kernel("daniell", c(3,3))
x1 <- kernapply(data[,5], kd)
24
0 200 400 600 800 1000
80
09
00
10
00
11
00
12
00
13
00
Index
da
ta
[, 5
]
Figura 19
Uma outra função implementada no R é a função loess()no quadro seguinte
apresento os passos efectuados:
Usando 30% da amostra para largura da banda.
> ll0<-loess(Close~.,data=SPtreino[,2:5],degree=0,span=0.3)
> preds.ll0<-predict(ll0,SPteste[,2:5])
> ll1<-loess(Close ~.,data=SPtreino[,2:5],degree=1,span=0.3)
> preds.ll1<-predict(ll1,SPteste[,2:5])
> ll2<-loess(Close ~.,data=SPtreino[,2:5],degree=2,span=0.3)
> preds.ll2<-predict(ll2,SPteste[,2:5])
preds.ll0
[1] 1242.691 1242.525 1243.093 1242.873
> preds.ll1
[1] 1286.186 1292.139 1287.472 1289.680
> preds.ll2
[1] 1286.342 1291.332 1289.350 1289.880
col="red",ylim=range(c(preds.ll1,preds.ll2),na.rm=T))
Cálculo do erro quadrático médio:
> mes.ll0<-mean((preds.ll0-SPteste[,5])^2)
> mes.ll1<-mean((preds.ll1-SPteste[,5])^2)
> mes.ll2<-mean((preds.ll2-SPteste[,5])^2)
> mes.ll0
[1] 2076.796
25
> mes.ll1
[1] 2.600855
> mes.ll2
[1] 3.848167
26
Conclusão
Após experimentar os vários métodos, apresentados anteriormente, os dois modelos
que melhor se ajustam à série do fecho do índice S&P500 são:
- Rede neuronal para size=20 e decay=0.001 em que o MSE foi de 2.9.
- Modelo de Kernel para um polinómio de grau 1 em que o MSE foi de 2.6.
Previsão dos valores de fecho do índice S&P500 para os dias 27, 28 de Fevereiro,
1,2 e 3 de Março
27-Fev- 1288.818
28-Fev-1288.278
1-Mar- 1287.755
2-Mar- 1287.233
3-Mar-1286.712
27
Bibliografia
[1]-“Data Mining with R”- Luis Torgo;
[2]-“Time Series Analysis and Its Applications”- Robert H. Shumway- Springer;
[3]-“A Course in Time Series Analysis”- Daniel Peña- Wiley-Interscience
Publication;
[4]-“Non-linear Time Series”- Howell Tong- Oxford Science Publications.