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 funcoes
#==============================================================================
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(fGarch)) # Para os modelos ARCH e GARCH
In [2]:
#==============================================================================
# Lendo os dados do arquivo BVSP.csv
#==============================================================================
BVSP <- read.csv("../dados/BVSP.csv",header=TRUE,sep=",")
print(head(BVSP))
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]:
#==============================================================================
# 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
In [4]:
#==============================================================================
# Grafico da serie de precos
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(15,7,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
plot(dias,preco/1000,lwd=2,col="black",main="",xlab="Dias",
ylab="Preço (milhares)",type="l",ylim=c(20,140),
xlim=c(as.Date("2011-01-01"),as.Date("2021-12-31")))
In [5]:
#==============================================================================
# Calculando o retorno de capitalizacao continua, ou simplesmente log-retorno
#==============================================================================
rt <- diff(log(preco)) # Calcula a primeira diferenca
N <- length(rt)
In [6]:
#==============================================================================
# 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.4,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,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 [7]:
#==============================================================================
# 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 [8]:
#==============================================================================
# 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 [9]:
#==============================================================================
# Calculando varios momentos amostrais
#==============================================================================
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 [10]:
#==============================================================================
# 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 [11]:
#==============================================================================
# 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 [12]:
#==============================================================================
# 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 [13]:
#==============================================================================
# Ajuste de modelos GARCH(1,1)
#==============================================================================
defaultW <- getOption("warn")
options(warn = -1)
garch.1.1 <- garchFit( rt ~ garch(1,1),data=rt,trace=F) #GARCH(1,1)
options(warn = defaultW)
summary(garch.1.1)
# Note que o teste de Jarque e Bera e o teste de Shapiro-Wilk rejeitam a
# hipotese nula de normalidade
Title:
GARCH Modelling
Call:
garchFit(formula = rt ~ garch(1, 1), data = rt, trace = F)
Mean and Variance Equation:
data ~ garch(1, 1)
<environment: 0x55fe4d214cc8>
[data = rt]
Conditional Distribution:
norm
Coefficient(s):
mu omega alpha1 beta1
4,4906e-04 1,0336e-05 8,4029e-02 8,7044e-01
Std. Errors:
based on Hessian
Error Analysis:
Estimate Std. Error t value Pr(>|t|)
mu 4,491e-04 2,642e-04 1,700 0,0892 .
omega 1,034e-05 2,237e-06 4,620 3,84e-06 ***
alpha1 8,403e-02 1,145e-02 7,341 2,12e-13 ***
beta1 8,704e-01 1,749e-02 49,757 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Log Likelihood:
7391 normalized: 2,828
Description:
Sun Sep 21 13:13:23 2025 by user: ralph
Standardised Residuals Tests:
Statistic p-Value
Jarque-Bera Test R Chi^2 823,6914 0,000000
Shapiro-Wilk Test R W 0,9825 0,000000
Ljung-Box Test R Q(10) 11,0414 0,354303
Ljung-Box Test R Q(15) 16,4086 0,355427
Ljung-Box Test R Q(20) 25,8237 0,171713
Ljung-Box Test R^2 Q(10) 27,2453 0,002381
Ljung-Box Test R^2 Q(15) 28,0668 0,021155
Ljung-Box Test R^2 Q(20) 30,0458 0,069116
LM Arch Test R TR^2 27,5924 0,006343
Information Criterion Statistics:
AIC BIC SIC HQIC
-5,654 -5,645 -5,654 -5,651
In [14]:
#==============================================================================
# Agora, analisaremos os residuos padronizados do modelo GARCH(1,1) ajustado.
# Verificacao das hipoteses
#==============================================================================
# Note que o teste de Jarque e Bera e o teste de Shapiro-Wilk rejeitam a
# hipotese nula de normalidade
# Calculo dos residuos padronizados
resid <- garch.1.1@residuals/volatility(garch.1.1)
In [15]:
#==============================================================================
# ACF dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(resid,xlab="Defasagem",main="",ylim=c(-0.2,1))
In [16]:
#==============================================================================
# Histograma dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
hist(resid,xlab="Resíduos",main="",prob=T,ylab="Densidade",col="orange")
In [17]:
#==============================================================================
# Boxplot dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
boxplot(resid,ylab="Resíduos",col="lightblue",ylim=c(-8,4))
In [18]:
#==============================================================================
# QQ-Norm dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
qqnorm(resid,main="",xlab="Percentis teóricos",xlim=c(-5,5),
ylim=c(-5,5),ylab="Percentis amostrais",col="red",pch=15)
qqline(resid,lwd=2)
In [19]:
#==============================================================================
# Histograma dos residuos padronizados ao quadrado
#==============================================================================
# Note que o teste de Ljung-Box nao rejeita a hipotese nula de autocorrelacao
# zero para a serie r_t (R) em varias escolhas de m (nivel de 5%).
# O mesmo nao pode ser dito sobre os retornos ao quadrado (R^2).
# Mesmo assim, podemos verificar a ACF dos residuos ao quadrado.
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(resid^2,xlab="Defasagem",main="",ylim=c(-0.2,1))
# Parece existir uma dependencia de defasagem 7 (e 10).
# Talvez seja algum efeito do dia da semana (utilizar regressao binaria).
# Conclusao momentanea: continuamos com o modelo GARCH(1,1) com erros normais.
In [20]:
#==============================================================================
# Calculando varios momentos amostrais
#==============================================================================
options(digits=4)
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)
print(STATS)
# Se os erros forem normais, entao os residuos devem ter media e mediana 0,
# desvio padrao 1, assimetria 0 e excesso de curtose 0.
# Conclusao: podemos considerar outras distribuicoes para o termo de erro.
r.media r.desvpad r.mediana r.assimetria r.exc.curtose [1,] -0,0214 0,9991 -0,008651 -0,419 2,615
In [21]:
#==============================================================================
# Considere por um momento o modelo GARCH(1,1) como bem ajustado.
# Previsao 10 passos a frente com intervalo de previsao de 95%.
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.1,cex.axis=1.1,lab=c(10,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
previsao <- predict(garch.1.1,n.ahead=10,plot=TRUE,conf=.95,nx=100)
In [22]:
print(previsao)
meanForecast meanError standardDeviation lowerInterval upperInterval 1 0,0004491 0,01635 0,01635 -0,03159 0,03249 2 0,0004491 0,01629 0,01629 -0,03148 0,03238 3 0,0004491 0,01624 0,01624 -0,03138 0,03228 4 0,0004491 0,01619 0,01619 -0,03128 0,03217 5 0,0004491 0,01614 0,01614 -0,03118 0,03208 6 0,0004491 0,01609 0,01609 -0,03109 0,03199 7 0,0004491 0,01605 0,01605 -0,03100 0,03190 8 0,0004491 0,01600 0,01600 -0,03091 0,03181 9 0,0004491 0,01596 0,01596 -0,03083 0,03173 10 0,0004491 0,01592 0,01592 -0,03076 0,03165
In [23]:
#==============================================================================
# Grafico das volatilidades estimadas pelo modelo GARCH(1,1)
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,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,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 [24]:
#==============================================================================
# Grafico das volatilidades estimadas pelo modelo GARCH(1,1)
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(15,8,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
plot(dias[2:n],volatility(garch.1.1),lwd=2,col="black",type="l",xlab="Dias",
ylab="Volatilidade",xlim=c(as.Date("2011-01-01"),as.Date("2021-12-31")),
ylim=c(0.01,0.08))
In [25]:
#==============================================================================
# Calculo do V@R
#==============================================================================
options(digits=10)
print(garch.1.1@fit$coef)
print("==========================================",quote=F)
print(volatility(garch.1.1)[N]^2) # sigma[N]^2
print(garch.1.1@residuals[N])
print("==========================================",quote=F)
r.pred <- as.numeric(garch.1.1@fit$coef[1])
sig2.pred <- sum(garch.1.1@fit$coef[2:4]*
c(1,garch.1.1@residuals[N]^2,volatility(garch.1.1)[N]^2))
sig.pred <- sqrt(sig2.pred)
print(r.pred)
print(sig.pred)
print("==========================================",quote=F)
# Quantil de 5%
VAR1a.ret <- r.pred + qnorm(0.05)*sig.pred
VAR1a.din <- (10^7)*VAR1a.ret
print(cbind(VAR1a.ret,VAR1a.din))
print("==========================================",quote=F)
# Quantil de 2,5%
VAR1b.ret <- r.pred + qnorm(0.025)*sig.pred
VAR1b.din <- (10^7)*VAR1b.ret
print(cbind(VAR1b.ret,VAR1b.din))
print("==========================================",quote=F)
# Quantil de 1%
VAR1c.ret <- r.pred + qnorm(0.01)*sig.pred
VAR1c.din <- (10^7)*VAR1c.ret
print(cbind(VAR1c.ret,VAR1c.din))
mu omega alpha1 beta1
4,490575253e-04 1,033603512e-05 8,402853197e-02 8,704406386e-01
[1] ==========================================
[1] 0,0002858724896
[1] -0,009811320931
[1] ==========================================
[1] 0,0004490575253
[1] 0,01634808318
[1] ==========================================
VAR1a.ret VAR1a.din
[1,] -0,02644114639 -264411,4639
[1] ==========================================
VAR1b.ret VAR1b.din
[1,] -0,03159259672 -315925,9672
[1] ==========================================
VAR1c.ret VAR1c.din
[1,] -0,03758227102 -375822,7102
In [26]:
#==============================================================================
# Calculo do V@R - em bloco
#==============================================================================
print(previsao)
print("==========================================",quote=F)
r.previsao <- previsao$meanForecast
sig.previsao <- previsao$standardDeviation
# Quantil de 5%
VAR1a.ret <- r.previsao + qnorm(0.05)*sig.previsao
VAR1a.din <- (10^7)*VAR1a.ret
print(cbind(VAR1a.ret,VAR1a.din))
print("==========================================",quote=F)
# Quantil de 2,5%
VAR1b.ret <- r.previsao + qnorm(0.025)*sig.previsao
VAR1b.din <- (10^7)*VAR1b.ret
print(cbind(VAR1b.ret,VAR1b.din))
print("==========================================",quote=F)
# Quantil de 1%
VAR1c.ret <- r.previsao + qnorm(0.01)*sig.previsao
VAR1c.din <- (10^7)*VAR1c.ret
print(cbind(VAR1c.ret,VAR1c.din))
meanForecast meanError standardDeviation lowerInterval upperInterval
1 0,0004490575253 0,01634808318 0,01634808318 -0,03159259672 0,03249071177
2 0,0004490575253 0,01629193964 0,01629193964 -0,03148255741 0,03238067246
3 0,0004490575253 0,01623817131 0,01623817131 -0,03137717341 0,03227528846
4 0,0004490575253 0,01618668450 0,01618668450 -0,03127626112 0,03217437617
5 0,0004490575253 0,01613738870 0,01613738870 -0,03117964314 0,03207775819
6 0,0004490575253 0,01609019652 0,01609019652 -0,03108714816 0,03198526321
7 0,0004490575253 0,01604502357 0,01604502357 -0,03099861080 0,03189672585
8 0,0004490575253 0,01600178843 0,01600178843 -0,03091387149 0,03181198654
9 0,0004490575253 0,01596041259 0,01596041259 -0,03083277633 0,03173089138
10 0,0004490575253 0,01592082032 0,01592082032 -0,03075517691 0,03165329196
[1] ==========================================
VAR1a.ret VAR1a.din
[1,] -0,02644114639 -264411,4639
[2,] -0,02634879848 -263487,9848
[3,] -0,02626035744 -262603,5744
[4,] -0,02617566918 -261756,6918
[5,] -0,02609458482 -260945,8482
[6,] -0,02601696058 -260169,6058
[7,] -0,02594265769 -259426,5769
[8,] -0,02587154222 -258715,4222
[9,] -0,02580348501 -258034,8501
[10,] -0,02573836153 -257383,6153
[1] ==========================================
VAR1b.ret VAR1b.din
[1,] -0,03159259672 -315925,9672
[2,] -0,03148255741 -314825,5741
[3,] -0,03137717341 -313771,7341
[4,] -0,03127626112 -312762,6112
[5,] -0,03117964314 -311796,4314
[6,] -0,03108714816 -310871,4816
[7,] -0,03099861080 -309986,1080
[8,] -0,03091387149 -309138,7149
[9,] -0,03083277633 -308327,7633
[10,] -0,03075517691 -307551,7691
[1] ==========================================
VAR1c.ret VAR1c.din
[1,] -0,03758227102 -375822,7102
[2,] -0,03745166162 -374516,6162
[3,] -0,03732657777 -373265,7777
[4,] -0,03720680154 -372068,0154
[5,] -0,03709212238 -370921,2238
[6,] -0,03698233694 -369823,3694
[7,] -0,03687724894 -368772,4894
[8,] -0,03677666898 -367766,6898
[9,] -0,03668041437 -366804,1437
[10,] -0,03658830898 -365883,0898
In [27]:
#==============================================================================
# Agora, ajustaremos o modelo GARCH(1,1) com erros t-Student
# "std" : t-Student
#==============================================================================
#GARCH(1,1)
defaultW <- getOption("warn")
options(warn = -1)
garch11.std <- garchFit( rt ~ garch(1,1),data=rt,cond.dist="std",trace=F)
options(warn = defaultW)
summary(garch11.std)
Title:
GARCH Modelling
Call:
garchFit(formula = rt ~ garch(1, 1), data = rt, cond.dist = "std",
trace = F)
Mean and Variance Equation:
data ~ garch(1, 1)
<environment: 0x55fe4fd37f50>
[data = rt]
Conditional Distribution:
std
Coefficient(s):
mu omega alpha1 beta1 shape
4,96769e-04 8,31776e-06 7,42065e-02 8,87683e-01 9,02450e+00
Std. Errors:
based on Hessian
Error Analysis:
Estimate Std. Error t value Pr(>|t|)
mu 4,9677e-04 2,5448e-04 1,9521 0,0509302 .
omega 8,3178e-06 2,1558e-06 3,8584 0,0001142 ***
alpha1 7,4207e-02 1,2222e-02 6,0718 1,265e-09 ***
beta1 8,8768e-01 1,8011e-02 49,2845 < 2,2e-16 ***
shape 9,0245e+00 1,3472e+00 6,6985 2,106e-11 ***
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Log Likelihood:
7439,131155 normalized: 2,846969443
Description:
Sun Sep 21 13:13:25 2025 by user: ralph
Standardised Residuals Tests:
Statistic p-Value
Jarque-Bera Test R Chi^2 858,1505552927 0,0000000000000
Shapiro-Wilk Test R W 0,9820154876 0,0000000000000
Ljung-Box Test R Q(10) 11,3084523833 0,3339970820185
Ljung-Box Test R Q(15) 16,8636321372 0,3270816656740
Ljung-Box Test R Q(20) 26,3304873222 0,1551722465741
Ljung-Box Test R^2 Q(10) 29,6384968380 0,0009813439448
Ljung-Box Test R^2 Q(15) 30,5998333312 0,0099331902268
Ljung-Box Test R^2 Q(20) 32,5627258677 0,0376590500165
LM Arch Test R TR^2 29,8053474797 0,0029870719429
Information Criterion Statistics:
AIC BIC SIC HQIC
-5,690111868 -5,678882908 -5,690119172 -5,686044374
In [28]:
#==============================================================================
# Agora, analisaremos os residuos padronizados do modelo GARCH(1,1) ajustado.
# Verificacao das hipoteses
#==============================================================================
# Calculo dos residuos padronizados
resid <- garch11.std@residuals/volatility(garch11.std)
graus <- garch11.std@fit$par[5]
In [29]:
#==============================================================================
# ACF dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(15,5,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(resid,xlab="Defasagem",main="",ylim=c(-0.2,1),lag.max=16)
In [30]:
#==============================================================================
# Histograma dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
hist(resid,xlab="Resíduos",main="",prob=T,ylab="Densidade",col="orange")
In [31]:
#==============================================================================
# Boxplot dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
boxplot(resid,ylab="Resíduos",col="red",ylim=c(-8,5))
In [32]:
#==============================================================================
# QQ-Norm dos residuos padronizados
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
qqplot(resid,rt(N,df=graus),main="",xlab="Percentis teóricos",
ylab="Percentis amostrais",col="red",pch=15,xlim=c(-6,6),ylim=c(-6,6))
In [33]:
#==============================================================================
# Histograma dos residuos padronizados ao quadrado
#==============================================================================
# Note que o teste de Ljung-Box nao rejeita a hipotese nula de autocorrelacao
# zero para a serie r_t (R) em varias escolhas de m (nivel de 5%).
# O mesmo nao pode ser dito sobre os retornos ao quadrado (R^2).
# Mesmo assim, podemos verificar a ACF dos residuos ao quadrado.
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(8,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(resid^2,xlab="Defasagem",main="",lag.max=16,ylim=c(-0.2,1))
# Parece existir uma dependencia de defasagem 7 (e 10).
# Talvez seja algum efeito do dia da semana (utilizar regressao binaria).
# Conclusao momentanea: continuamos com o modelo GARCH(1,1) com erros t-Student.
In [34]:
#==============================================================================
# Considere por um momento o modelo GARCH(1,1) como bem ajustado.
# Previsao 10 passos a frente com intervalo de previsao de 95%.
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(10,8,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
previsaost <- predict(garch11.std,n.ahead=10,plot=TRUE,conf=.95,nx=100)
In [35]:
#==============================================================================
# Grafico das volatilidades estimadas pelo modelo GARCH(1,1) STD
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,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,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 [36]:
#==============================================================================
# Grafico das volatilidades estimadas pelo modelo GARCH(1,1)
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.4,cex.axis=1.5,lab=c(15,8,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
plot(dias[2:n],volatility(garch11.std),lwd=2,col="black",type="l",xlab="Dias",
ylab="Volatilidade",xlim=c(as.Date("2011-01-01"),as.Date("2021-12-31")),
ylim=c(0.01,0.08))
In [37]:
#==============================================================================
# Calculo do V@R
#==============================================================================
options(digits=10)
print(garch11.std@fit$coef)
print("==========================================",quote=F)
print(volatility(garch11.std)[N]^2) # sigma[N]^2
print(garch11.std@residuals[N])
print("==========================================",quote=F)
r.pred <- as.numeric(garch11.std@fit$coef[1])
sig2.pred <- sum(garch11.std@fit$coef[2:4]*
c(1,garch11.std@residuals[N]^2,volatility(garch11.std)[N]^2))
sig.pred <- sqrt(sig2.pred)
print(r.pred)
print(sig.pred)
print("==========================================",quote=F)
# Quantil de 5%
VAR2a.ret <- r.pred + qt(0.05,df=garch11.std@fit$coef[5])*sig.pred
VAR2a.din <- (10^7)*VAR2a.ret
print(cbind(VAR2a.ret,VAR2a.din))
print("==========================================",quote=F)
# Quantil de 2,5%
VAR2b.ret <- r.pred + qt(0.025,df=garch11.std@fit$coef[5])*sig.pred
VAR2b.din <- (10^7)*VAR2b.ret
print(cbind(VAR2b.ret,VAR2b.din))
print("==========================================",quote=F)
# Quantil de 1%
VAR2c.ret <- r.pred + qt(0.01,df=garch11.std@fit$coef[5])*sig.pred
VAR2c.din <- (10^7)*VAR2c.ret
print(cbind(VAR2c.ret,VAR2c.din))
mu omega alpha1 beta1 shape
4,967693489e-04 8,317760510e-06 7,420650596e-02 8,876827728e-01 9,024503781e+00
[1] ==========================================
[1] 0,0002704784588
[1] -0,009859032754
[1] ==========================================
[1] 0,0004967693489
[1] 0,0159884252
[1] ==========================================
VAR2a.ret VAR2a.din
[1,] -0,02880276191 -288027,6191
[1] ==========================================
VAR2b.ret VAR2b.din
[1,] -0,03565659668 -356565,9668
[1] ==========================================
VAR2c.ret VAR2c.din
[1,] -0,04458811082 -445881,1082
In [38]:
#==============================================================================
# Calculo do V@R - em bloco
#==============================================================================
print(previsaost)
print("==========================================",quote=F)
r.previsaost <- previsaost$meanForecast
sig.previsaost <- previsaost$standardDeviation
# Quantil de 5%
VAR2a.ret <- r.previsaost + qt(0.05,df=garch11.std@fit$coef[5])*sig.previsaost
VAR2a.din <- (10^7)*VAR2a.ret
print(cbind(VAR2a.ret,VAR2a.din))
print("==========================================",quote=F)
# Quantil de 2,5%
VAR2b.ret <- r.previsaost + qt(0.025,df=garch11.std@fit$coef[5])*sig.previsaost
VAR2b.din <- (10^7)*VAR2b.ret
print(cbind(VAR2b.ret,VAR2b.din))
print("==========================================",quote=F)
# Quantil de 1%
VAR2c.ret <- r.previsaost + qt(0.01,df=garch11.std@fit$coef[5])*sig.previsaost
VAR2c.din <- (10^7)*VAR2c.ret
print(cbind(VAR2c.ret,VAR2c.din))
meanForecast meanError standardDeviation lowerInterval upperInterval
1 0,0004967693489 0,01598842520 0,01598842520 -0,03139986779 0,03239340649
2 0,0004967693489 0,01594381595 0,01594381595 -0,03131087309 0,03230441179
3 0,0004967693489 0,01590078871 0,01590078871 -0,03122503447 0,03221857317
4 0,0004967693489 0,01585929111 0,01585929111 -0,03114224748 0,03213578618
5 0,0004967693489 0,01581927231 0,01581927231 -0,03106241065 0,03205594935
6 0,0004967693489 0,01578068290 0,01578068290 -0,03098542543 0,03197896412
7 0,0004967693489 0,01574347491 0,01574347491 -0,03091119611 0,03190473481
8 0,0004967693489 0,01570760178 0,01570760178 -0,03083962982 0,03183316852
9 0,0004967693489 0,01567301832 0,01567301832 -0,03077063642 0,03176417512
10 0,0004967693489 0,01563968071 0,01563968071 -0,03070412845 0,03169766715
[1] ==========================================
VAR2a.ret VAR2a.din
[1,] -0,02880276191 -288027,6191
[2,] -0,02872101339 -287210,1339
[3,] -0,02864216397 -286421,6397
[4,] -0,02856611771 -285661,1771
[5,] -0,02849278140 -284927,8140
[6,] -0,02842206451 -284220,6451
[7,] -0,02835387914 -283538,7914
[8,] -0,02828813996 -282881,3996
[9,] -0,02822476417 -282247,6417
[10,] -0,02816367145 -281636,7145
[1] ==========================================
VAR2b.ret VAR2b.din
[1,] -0,03565659668 -356565,9668
[2,] -0,03555572530 -355557,2530
[3,] -0,03545843118 -354584,3118
[4,] -0,03536459595 -353645,9595
[5,] -0,03527410459 -352741,0459
[6,] -0,03518684539 -351868,4539
[7,] -0,03510270989 -351027,0989
[8,] -0,03502159280 -350215,9280
[9,] -0,03494339196 -349433,9196
[10,] -0,03486800825 -348680,0825
[1] ==========================================
VAR2c.ret VAR2c.din
[1,] -0,04458811082 -445881,1082
[2,] -0,04446231966 -444623,1966
[3,] -0,04434098950 -443409,8950
[4,] -0,04422397273 -442239,7273
[5,] -0,04411112591 -441111,2591
[6,] -0,04400230975 -440023,0975
[7,] -0,04389738898 -438973,8898
[8,] -0,04379623231 -437962,3231
[9,] -0,04369871233 -436987,1233
[10,] -0,04360470544 -436047,0544
In [39]:
#==============================================================================
# Comparando os resultados do ajuste normal e t-Student
#==============================================================================
print(previsao)
print("==========================================",quote=F)
print(previsaost)
print("==========================================",quote=F)
# Quantil de 5%
print("Quantil de 5%")
print(cbind(VAR1a.ret,VAR1a.din,VAR2a.ret,VAR2a.din))
print("==========================================",quote=F)
# Quantil de 2,5%
print("Quantil de 2,5%")
print(cbind(VAR1b.ret,VAR1b.din,VAR2b.ret,VAR2b.din))
print("==========================================",quote=F)
# Quantil de 1%
print("Quantil de 1%")
print(cbind(VAR1c.ret,VAR1c.din,VAR2c.ret,VAR2c.din))
meanForecast meanError standardDeviation lowerInterval upperInterval
1 0,0004490575253 0,01634808318 0,01634808318 -0,03159259672 0,03249071177
2 0,0004490575253 0,01629193964 0,01629193964 -0,03148255741 0,03238067246
3 0,0004490575253 0,01623817131 0,01623817131 -0,03137717341 0,03227528846
4 0,0004490575253 0,01618668450 0,01618668450 -0,03127626112 0,03217437617
5 0,0004490575253 0,01613738870 0,01613738870 -0,03117964314 0,03207775819
6 0,0004490575253 0,01609019652 0,01609019652 -0,03108714816 0,03198526321
7 0,0004490575253 0,01604502357 0,01604502357 -0,03099861080 0,03189672585
8 0,0004490575253 0,01600178843 0,01600178843 -0,03091387149 0,03181198654
9 0,0004490575253 0,01596041259 0,01596041259 -0,03083277633 0,03173089138
10 0,0004490575253 0,01592082032 0,01592082032 -0,03075517691 0,03165329196
[1] ==========================================
meanForecast meanError standardDeviation lowerInterval upperInterval
1 0,0004967693489 0,01598842520 0,01598842520 -0,03139986779 0,03239340649
2 0,0004967693489 0,01594381595 0,01594381595 -0,03131087309 0,03230441179
3 0,0004967693489 0,01590078871 0,01590078871 -0,03122503447 0,03221857317
4 0,0004967693489 0,01585929111 0,01585929111 -0,03114224748 0,03213578618
5 0,0004967693489 0,01581927231 0,01581927231 -0,03106241065 0,03205594935
6 0,0004967693489 0,01578068290 0,01578068290 -0,03098542543 0,03197896412
7 0,0004967693489 0,01574347491 0,01574347491 -0,03091119611 0,03190473481
8 0,0004967693489 0,01570760178 0,01570760178 -0,03083962982 0,03183316852
9 0,0004967693489 0,01567301832 0,01567301832 -0,03077063642 0,03176417512
10 0,0004967693489 0,01563968071 0,01563968071 -0,03070412845 0,03169766715
[1] ==========================================
[1] "Quantil de 5%"
VAR1a.ret VAR1a.din VAR2a.ret VAR2a.din
[1,] -0,02644114639 -264411,4639 -0,02880276191 -288027,6191
[2,] -0,02634879848 -263487,9848 -0,02872101339 -287210,1339
[3,] -0,02626035744 -262603,5744 -0,02864216397 -286421,6397
[4,] -0,02617566918 -261756,6918 -0,02856611771 -285661,1771
[5,] -0,02609458482 -260945,8482 -0,02849278140 -284927,8140
[6,] -0,02601696058 -260169,6058 -0,02842206451 -284220,6451
[7,] -0,02594265769 -259426,5769 -0,02835387914 -283538,7914
[8,] -0,02587154222 -258715,4222 -0,02828813996 -282881,3996
[9,] -0,02580348501 -258034,8501 -0,02822476417 -282247,6417
[10,] -0,02573836153 -257383,6153 -0,02816367145 -281636,7145
[1] ==========================================
[1] "Quantil de 2,5%"
VAR1b.ret VAR1b.din VAR2b.ret VAR2b.din
[1,] -0,03159259672 -315925,9672 -0,03565659668 -356565,9668
[2,] -0,03148255741 -314825,5741 -0,03555572530 -355557,2530
[3,] -0,03137717341 -313771,7341 -0,03545843118 -354584,3118
[4,] -0,03127626112 -312762,6112 -0,03536459595 -353645,9595
[5,] -0,03117964314 -311796,4314 -0,03527410459 -352741,0459
[6,] -0,03108714816 -310871,4816 -0,03518684539 -351868,4539
[7,] -0,03099861080 -309986,1080 -0,03510270989 -351027,0989
[8,] -0,03091387149 -309138,7149 -0,03502159280 -350215,9280
[9,] -0,03083277633 -308327,7633 -0,03494339196 -349433,9196
[10,] -0,03075517691 -307551,7691 -0,03486800825 -348680,0825
[1] ==========================================
[1] "Quantil de 1%"
VAR1c.ret VAR1c.din VAR2c.ret VAR2c.din
[1,] -0,03758227102 -375822,7102 -0,04458811082 -445881,1082
[2,] -0,03745166162 -374516,6162 -0,04446231966 -444623,1966
[3,] -0,03732657777 -373265,7777 -0,04434098950 -443409,8950
[4,] -0,03720680154 -372068,0154 -0,04422397273 -442239,7273
[5,] -0,03709212238 -370921,2238 -0,04411112591 -441111,2591
[6,] -0,03698233694 -369823,3694 -0,04400230975 -440023,0975
[7,] -0,03687724894 -368772,4894 -0,04389738898 -438973,8898
[8,] -0,03677666898 -367766,6898 -0,04379623231 -437962,3231
[9,] -0,03668041437 -366804,1437 -0,04369871233 -436987,1233
[10,] -0,03658830898 -365883,0898 -0,04360470544 -436047,0544