In [1]:
options(OutDec = ",")
#==============================================================================
# Descricao
#==============================================================================
# Os dados a seguir correspondem aos indices diarios da BOVESPA de 3 de
# janeiro de 2011 a 10 de setembro de 2021.
#
# Os dados foram obtidos do sitio http://br.financas.yahoo.com/
#
# Concentraremos nos indices (precos) de fechamento ("Close")
#
#==============================================================================
# Carregando pacotes
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(tseries)) # Teste de Jarque-Bera: 'jarque.bera.test'
silence(library(psych)) # Estatistica descritiva: 'describe'
silence(library(latticeExtra)) # Para construir dois graficos em um so'
In [2]:
#==============================================================================
# Lendo os dados do arquivo BVSP.csv
#==============================================================================
BVSP <- read.csv("../dados/BVSP.csv",header=TRUE,sep=",")
print(head(BVSP))
# Copiando os indices de fechamento como a serie de precos
preco <- BVSP$Close
dias <- as.Date(BVSP$Date) # Dias
n <- length(preco) # Numero de observacoes
Date Open High Low Close Adj.Close Volume 1 2011-01-03 69310 70471 69305 69962 69962 1862400 2 2011-01-04 69962 70318 69560 70318 70318 2427200 3 2011-01-05 70311 71173 69802 71091 71091 2309200 4 2011-01-06 71093 71167 70469 70579 70579 2546000 5 2011-01-07 70580 70783 69718 70057 70057 1761000 6 2011-01-10 70056 70133 69666 70127 70127 1610800
In [3]:
#==============================================================================
# Grafico da serie de precos
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(15,5,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
plot(dias,preco/1000,lwd=2,type="l",col="black",main="",xlab="Dias",
ylab="Preço (milhares)",ylim=c(20,140),
xlim=c(as.Date("2011-01-01"),as.Date("2021-12-31")))
In [4]:
#==============================================================================
# Calculando o retorno de capitalizacao continua, ou simplesmente log-retorno
#==============================================================================
rt <- diff(log(preco)) # Calcula a primeira diferenca
In [5]:
# Grafico da serie de log retornos: r_t = log(1+R_t) = log(P_t)-log(P_{t-1})
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.4,lab=c(15,5,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
plot(dias[2:n],rt,lwd=2,col="black",main="",xlab="Dias",ylab="Log Retorno",
type="l",xlim=c(as.Date("2011-01-01"),as.Date("2021-12-31")),
ylim=c(-0.16,0.16))
In [6]:
#==============================================================================
# A serie de precos e log retornos em um unico grafico
#==============================================================================
obj1 <- xyplot(rt ~ dias,type="l",lwd=2,col=1,xlab="Dias",
ylab="Log Retorno")
obj2 <- xyplot(preco ~ dias,type="l",lwd=2,col="red",ylab="Preços")
doubleYScale(obj1,obj2,add.ylab2=TRUE,use.style=FALSE)
In [7]:
#==============================================================================
# Grafico da serie de log retornos ao quadrado
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(15,5,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
plot(dias[2:n],rt^2,lwd=2,col=1,main="",xlab="Dias",
ylab="Log Retorno ao Quadrado",type="l",
xlim=c(as.Date("2011-01-01"),as.Date("2021-12-31")))
In [8]:
#==============================================================================
# ACF dos log retornos
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(rt,lwd=2,col=1,lag.max=16,xlab="Defasagem",main="",ylim=c(-0.2,1))
In [9]:
#==============================================================================
# ACF dos log retornos ao quadrado
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(rt^2,lwd=2,col=1,lag.max=16,xlab="Defasagem",main="",ylim=c(-0.2,1))
In [10]:
options(digits=4)
r.describe <- describe(rt)
r.media <- r.describe[,3]
r.desvpad <- r.describe[,4]
r.mediana <- r.describe[,5]
r.assimetria <- r.describe[,11]
r.exc.curtose <- r.describe[,12]
STATS <- cbind(r.media,r.desvpad,r.mediana,r.assimetria,
r.exc.curtose)
colnames(STATS) <- c("media","desvpad","mediana","assimetria","exc.curt")
rownames(STATS) <- "IBOVESPA"
print(STATS)
media desvpad mediana assimetria exc.curt IBOVESPA 0,0001878 0,01616 0,0003192 -0,8946 12,23
In [11]:
#==============================================================================
# Histograma dos log retornos
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(12,5,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
hist(rt,lwd=2,col=0,main="",xlab="Log Retornos",ylab="Densidade",prob=T,
nclass=50,xlim=c(-0.17,0.17))
xseq <- seq(r.media-7*r.desvpad,r.media+7*r.desvpad,length=1000)
yseq <- dnorm(xseq,r.media,r.desvpad)
lines(xseq,yseq,lwd=3,col="red")
text(-0.09,25,"Certamente não normal",cex=1.2)
In [12]:
#==============================================================================
# Boxplot dos log retornos
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
boxplot(rt,lwd=2,col="red",main="",ylab="Log Retornos",pch=16,
ylim=c(-0.2,0.2))
In [13]:
#==============================================================================
# Normal QQ-plot dos log retornos
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
qqnorm(rt,lwd=2,col="red",main="",xlab="Quantis teóricos",
ylab="Quantis amostrais",pch=16,xlim=c(-4,4),ylim=c(-0.2,0.2))
qqline(rt,lwd=2,col="black",xlim=c(-4,4),xpd=F)
In [14]:
#==============================================================================
# Testando a hipotese de distribuicao normal para os log retornos
#==============================================================================
# Lembre-se de verificar de novo a assimetria e a curtose
print(STATS[1,4:5])
print("========================")
# Nivel de significancia dos testes
alpha <- 0.05
N <- length(rt) # Tamanho da serie de log retornos
# Testando a assimetria zero
print("Teste de assimetria")
tSr <- STATS[1,4]/sqrt(6/N) # Razao t-Student
print(abs(tSr) > qnorm(1-alpha/2)) # Testanto (TRUE=Rejeitar, FALSE=Nao rejeitar)
valorpSr <- 2*(1-pnorm(abs(tSr))) # Valor p
print(valorpSr)
print("========================")
# Testando a excesso de curtose zero
print("Teste de curtose")
tKr <- STATS[1,5]/sqrt(24/N) # Razao t-Student
print(abs(tKr) > qnorm(1-alpha/2)) # Testanto (TRUE=Rejeitar, FALSE=Nao rejeitar)
valorpKr <- 2*(1-pnorm(abs(tKr))) # Valor p
print(valorpKr)
print("========================")
# Testando a normalidade - teste de Jarque Bera
print("Teste de Jarque-Bera: normalidade")
JB <- tSr^2 + tKr^2
print(JB > qchisq(1-alpha,2)) # Testanto (TRUE=Rejeitar, FALSE=Nao rejeitar)
print(jarque.bera.test(rt)) # Verificar o valor p (mesma conclusao)
assimetria exc.curt -0,8946 12,2347 [1] "========================" [1] "Teste de assimetria" [1] TRUE [1] 0 [1] "========================" [1] "Teste de curtose" [1] TRUE [1] 0 [1] "========================" [1] "Teste de Jarque-Bera: normalidade" [1] TRUE Jarque Bera Test data: rt X-squared = 16677, df = 2, p-value <2e-16
In [15]:
#==============================================================================
# Teste para efeitos ARCH
#==============================================================================
# Passo 1 em nossa analise
#==============================================================================
# Vimos pela ACF dos retornos uma coorelacao serial fraca.
# Teste de Ljung-Box para os retornos
# Lembre-se da hipotese nula (H_0): independencia dos retornos
ell <- floor(log(N)) # ell =~ log(N) como valor mais pratico
print(ell)
for(i in 1:ell)
print(Box.test(rt, lag = i, type = "Ljung-Box"))
# Conclusao: rejeitamos a hipotese nula de independencia dos retornos.
# Precisamos modelar a media. Utilizaremos um modelo AR.
[1] 7 Box-Ljung test data: rt X-squared = 26, df = 1, p-value = 4e-07 Box-Ljung test data: rt X-squared = 31, df = 2, p-value = 2e-07 Box-Ljung test data: rt X-squared = 31, df = 3, p-value = 7e-07 Box-Ljung test data: rt X-squared = 32, df = 4, p-value = 2e-06 Box-Ljung test data: rt X-squared = 41, df = 5, p-value = 8e-08 Box-Ljung test data: rt X-squared = 55, df = 6, p-value = 4e-10 Box-Ljung test data: rt X-squared = 75, df = 7, p-value = 1e-13
In [16]:
#==============================================================================
# Passo 2 em nossa analise
#==============================================================================
# Ajuste de um modelo AR
r.ar <- ar.mle(rt)
p <- r.ar$order
a <- r.ar$res[(p+1):N] # residuos nao padronizados
n2 <- length(a)
print(r.ar) # O melhor modelo foi um AR(10) segundo o AIC
Call:
ar.mle(x = rt)
Coefficients:
1 2 3 4 5 6 7 8 9 10
-0,085 0,032 0,006 -0,013 0,045 -0,054 0,069 -0,026 0,025 0,041
Order selected 10 sigma^2 estimated as 0,000254
In [17]:
#==============================================================================
# Verificando a ACF dos log retornos
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(a,xlab="Defasagem",main="",lwd=2,col=1,lag.max=15,ylim=c(-0.2,1))
In [18]:
#==============================================================================
# Calculando varios momentos amostrais dos residuos padronizados
# do modelo AR(10)
#==============================================================================
options(digits=4)
sigma2 <- r.ar$var # variancia do erro estimada
sigma <- sqrt(sigma2) # desvio padrao do erro estimado
resid <- a/sigma # residuos padronizados
r.describe <- describe(resid)
r.media <- r.describe[,3]
r.desvpad <- r.describe[,4]
r.mediana <- r.describe[,5]
r.assimetria <- r.describe[,11]
r.exc.curtose <- r.describe[,12]
STATS <- cbind(r.media,r.desvpad,r.mediana,r.assimetria,
r.exc.curtose)
colnames(STATS) <- c("media","desvpad","mediana","assimetria","exc.curt")
rownames(STATS) <- "IBOVESPA"
print(STATS)
# Se os erros forem normais, entao os residuos devem ter
# media e mediana 0 (zero), desvio padrao 1 (um),
# assimetria 0 (zero) e excesso de curtose 0 (zero).
# Possivelmente a presenca de nao normalidade.
media desvpad mediana assimetria exc.curt IBOVESPA -0,0002323 1,002 0,0007451 -0,7339 10,34
In [19]:
#==============================================================================
# Passo 3 em nossa analise
#==============================================================================
# Teste de Ljung-Box para os residuos (resid) do modelo AR(10)
# Lembre-se da hipotese nula (H_0): independencia dos residuos
for(i in (p+1):(p+ell+1))
print(Box.test(resid, lag = i, type = "Ljung-Box",fitdf = p))
# Agora, os p-valores sao todos grandes.
# Conclusao: nao rejeitamos a hipotese nula de independencia dos residuos
Box-Ljung test data: resid X-squared = 1,4, df = 1, p-value = 0,2 Box-Ljung test data: resid X-squared = 1,5, df = 2, p-value = 0,5 Box-Ljung test data: resid X-squared = 5,1, df = 3, p-value = 0,2 Box-Ljung test data: resid X-squared = 7,2, df = 4, p-value = 0,1 Box-Ljung test data: resid X-squared = 7,3, df = 5, p-value = 0,2 Box-Ljung test data: resid X-squared = 10, df = 6, p-value = 0,1 Box-Ljung test data: resid X-squared = 10, df = 7, p-value = 0,2 Box-Ljung test data: resid X-squared = 13, df = 8, p-value = 0,1
In [20]:
#==============================================================================
# Passo 4 em nossa analise
#==============================================================================
# Teste para efeitos ARCH
# Teste de Ljung-Box para os residuos ao quadrado do modelo AR(10)
# Lembre-se da hipotese nula (H_0): independencia dos residuos ao quadrado
resid2 <- resid^2
for(i in 1:ell)
print(Box.test(resid2, lag = i, type = "Ljung-Box"))
# Os p-valores sao todos pequenos.
# Conclusao: rejeitamos a hipotese nula de independencia dos residuos
# ao quadrado.
Box-Ljung test data: resid2 X-squared = 424, df = 1, p-value <2e-16 Box-Ljung test data: resid2 X-squared = 882, df = 2, p-value <2e-16 Box-Ljung test data: resid2 X-squared = 1159, df = 3, p-value <2e-16 Box-Ljung test data: resid2 X-squared = 1351, df = 4, p-value <2e-16 Box-Ljung test data: resid2 X-squared = 1504, df = 5, p-value <2e-16 Box-Ljung test data: resid2 X-squared = 1530, df = 6, p-value <2e-16 Box-Ljung test data: resid2 X-squared = 1738, df = 7, p-value <2e-16
In [21]:
#==============================================================================
# Verificando a ACF dos retornos ao quadrado
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(resid2,xlab="Defasagem",main="",lwd=2,col=1,lag.max=16,ylim=c(-0.2,1))
