Exercícios de TVE e Cópulas 
Profa Cira
20/02/2021
Os dados para esses exercícios foram obtidos em https://finance.yahoo.com/quote/PBR/history?p=PBR . (Petróleo Brasileiro S.A. - Petrobras (PBR) Feb 09, 2020 - Feb 09, 2021).
1 Exercício 1
Estimação da estável - Método dos quantis
1.1 Pacotes Utilizados
## Instalando e Carregando os pacotes (bibliotecas)
lbs<-c('tidyverse','lubridate','copula','fBasics','StableEstim','stabledist','DT','kableExtra','PerformanceAnalytics')
req <- substitute(require(x, character.only = TRUE))
sapply(lbs, function(x) eval(req) || {install.packages(x,dependencies = TRUE); eval(req)})## tidyverse lubridate copula
## TRUE TRUE TRUE
## fBasics StableEstim stabledist
## TRUE TRUE TRUE
## DT kableExtra PerformanceAnalytics
## TRUE TRUE TRUE1.2 Leitura e Apresentação dos Dados
# Leitura dos dados
PBR1<-read.csv(file.choose())
PBR1$retorno<-log(PBR1$High/dplyr::lag(PBR1$High))
head(PBR1$retorno)## [1] NA 0.025887004 0.010036886 -0.007350518 -0.004032264
## [6] -0.010832875PBR1<-PBR1[is.finite(PBR1$retorno), ]
kable(as.data.frame(PBR1[1:10,])) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Date | Open | High | Low | Close | Adj.Close | Volume | retorno | |
|---|---|---|---|---|---|---|---|---|
| 2 | 2020-02-11 | 14.73 | 14.87 | 14.63 | 14.82 | 14.68136 | 26408400 | 0.0258870 |
| 3 | 2020-02-12 | 14.81 | 15.02 | 14.81 | 14.92 | 14.78043 | 24175400 | 0.0100369 |
| 4 | 2020-02-13 | 14.91 | 14.91 | 14.56 | 14.63 | 14.49314 | 21952600 | -0.0073505 |
| 5 | 2020-02-14 | 14.83 | 14.85 | 14.58 | 14.64 | 14.50305 | 14207900 | -0.0040323 |
| 6 | 2020-02-18 | 14.32 | 14.69 | 14.31 | 14.56 | 14.42379 | 20754800 | -0.0108329 |
| 7 | 2020-02-19 | 14.77 | 14.94 | 14.73 | 14.89 | 14.75071 | 16214300 | 0.0168752 |
| 8 | 2020-02-20 | 14.77 | 14.93 | 14.37 | 14.40 | 14.26529 | 21453600 | -0.0006696 |
| 9 | 2020-02-21 | 14.17 | 14.18 | 13.89 | 14.03 | 13.89875 | 22637200 | -0.0515401 |
| 10 | 2020-02-24 | 13.14 | 13.28 | 12.71 | 13.08 | 12.95764 | 38318000 | -0.0655734 |
| 11 | 2020-02-25 | 13.20 | 13.30 | 12.62 | 12.82 | 12.70007 | 26389700 | 0.0015049 |
1.3 Estatísticas do retorno
ret<-PBR1$retorno
summary(ret)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.277014 -0.019388 0.000000 -0.001083 0.021181 0.149230media<-mean(ret)
varianc<-var(ret)
dp<-sqrt(varianc)
setNames(c(media,varianc,dp),c("Média "," Variância "," Desvio Padrão"))## Média Variância Desvio Padrão
## -0.001083401 0.002316954 0.0481347481.4 Gráficos do retorno
par(mfrow=c(1,2))
#==========================================================================#
plot(ret,type="l",ylab="retorno")
#============================ HISTOGRAMA =================================#
hist(ret, n = 50, probability = TRUE, border = "white",
col = "steelblue", main="", ylim=c(0,28), xlab="log Retorno diário PBR",
ylab="Densidade")
1.5 Ajustes dos Parâmetros - Método Quantis
(st1<-stableFit(ret, "q",doplot = TRUE))
##
## Title:
## Stable Parameter Estimation
##
## Call:
## .qStableFit(x = x, doplot = doplot, title = title, description = description)
##
## Model:
## Stable Distribution
##
## Estimated Parameter(s):
## alpha beta gamma delta
## 1.4860000000 0.0490000000 0.0209139098 -0.0002736295
##
## Description:
## Tue Mar 09 14:10:22 2021 by user: adolfoamds1.6 Extração das estimativas dos parâmetros
Extrai-se com o comando result@fit$estimate
(par.est<-st1@fit$estimate)## alpha beta gamma delta
## 1.4860000000 0.0490000000 0.0209139098 -0.0002736295alpha.est<-par.est[1] ; beta.est<-par.est[2]
gamma.est<-par.est[3] ; delta.est<-par.est[4]1.7 Cálculo do Value at Risk (VaR)*
O cálculo do VaR foi obtido através dos pacotes PerformanceAnalytics e stabledist.
Com o auxílio da função sapply do próprio R, apresentamos os valores do VaR para diversos níveis de confiança.
## Nível de Confianca
p<-c(0.90,0.95,0.975,0.99,0.999)
## VaR historico
V.hist<-sapply(p, VaR, R=ret, method="historical")
## VaR normal
V.normal<-sapply(p, VaR, R=ret, method="gaussian")
## VaR alpha-estavel
V.alpha<-sapply(p, qstable,alpha= alpha.est,beta= beta.est,gamma= gamma.est,delta= delta.est,pm=0)
## Apresentação dos VaRs obtidos
df_info<-data.frame(p,V.hist,V.normal,V.alpha)
names(df_info)<-c("Nível de Confiança","VaR Histórico","VaR Normal","VaR Alpha-Estável")
kable(df_info) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Nível de Confiança | VaR Histórico | VaR Normal | VaR Alpha-Estável |
|---|---|---|---|
| 0.900 | -0.0437765 | -0.0626476 | 0.0441451 |
| 0.950 | -0.0637639 | -0.0801001 | 0.0662240 |
| 0.975 | -0.1017346 | -0.0952377 | 0.0985055 |
| 0.990 | -0.1489494 | -0.1128383 | 0.1721705 |
| 0.999 | -0.2764062 | -0.1495343 | 0.7772120 |
1.8 Densidade alpha-estável
x<-ret ; hist(ret[-1],n=50, prob=T,ylim=c(0,28),main="Histograma do Retorno",xlab="retorno")
lines(seq(min(x,na.rm=T),max(x,na.rm=T),length=1000),dnorm(seq(min(x,na.rm=T),
max(x,na.rm=T),length=1000),mean(x,na.rm=T),sd(x,na.rm=T)),lwd=2)
curve(dstable(x, alpha=alpha.est, beta = beta.est,
gamma = gamma.est, delta = delta.est), -1, 1,add=TRUE)
2 Exercício 2
2.1 Pacotes Utilizados
## Instalando e Carregando os pacotes (bibliotecas)
lbs<-c('extRemes','ismev','evmix','copula','extremis','DT','kableExtra','VGAM','evd','fExtremes')
req <- substitute(require(x, character.only = TRUE))
sapply(lbs, function(x) eval(req) || {install.packages(x,dependencies = TRUE); eval(req)})## extRemes ismev evmix copula extremis DT kableExtra
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## VGAM evd fExtremes
## TRUE TRUE TRUE2.2 Leitura dos dados
# Leitura dos dados
PBR1<-read.csv(file.choose())
PBR1$retorno<-log(PBR1$High/dplyr::lag(PBR1$High)) ; ret<-PBR1$retorno ################## Bloco Máximo ###########
plot.new() ; par(mfrow=c(1,2))
N<-length(ret) ; n<-10
tau<-floor(N/n)
m<-numeric(tau) ; j<-1
for (i in 1:tau){
m[i]<-max(ret[j:(j+n-1)])
j<-j+n }
m<-m[-1]
{hist(m, prob=T,ylim=c(0,28))
lines(density(m))}
plot(m, type="l",main="Some Extremes") #Some extremes
p<-c((1:tau)/(tau+1)) ; ginv<- -log(-log(p))
qqplot(m,ginv,xlab="Quantis empíricos",ylab="Quantis da Gumbel",main="qqplot") ; grid()
2.3 Estimação dos Parâmetros
#################### Estimando os parametros ############
## R MLE - GEV
fitmv = gevFit(m, type ="mle")
fitmv##
## Title:
## GEV Parameter Estimation
##
## Call:
## gevFit(x = m, type = "mle")
##
## Estimation Type:
## gev mle
##
## Estimated Parameters:
## xi mu beta
## 0.24399924 0.04047920 0.02300403
##
## Description
## Tue Mar 09 14:10:32 2021 par(mfrow = c(2, 2))
summary(fitmv)##
## Title:
## GEV Parameter Estimation
##
## Call:
## gevFit(x = m, type = "mle")
##
## Estimation Type:
## gev mle
##
## Estimated Parameters:
## xi mu beta
## 0.24399924 0.04047920 0.02300403
##
## Standard Deviations:
## xi mu beta
## 0.242155030 0.005748520 0.004717997
##
## Log-Likelihood Value:
## -49.30563
##
## Type of Convergence:
## 0
##
## Description
## Tue Mar 09 14:10:32 2021
## R PWM - GEV
fitpwm = gevFit(m, type ="pwm")
fitpwm##
## Title:
## GEV Parameter Estimation
##
## Call:
## gevFit(x = m, type = "pwm")
##
## Estimation Type:
## gev pwm
##
## Estimated Parameters:
## xi mu beta
## 0.12583692 0.04119796 0.02594432
##
## Description
## Tue Mar 09 14:10:32 2021 par(mfrow = c(2, 2))
summary(fitpwm)##
## Title:
## GEV Parameter Estimation
##
## Call:
## gevFit(x = m, type = "pwm")
##
## Estimation Type:
## gev pwm
##
## Estimated Parameters:
## xi mu beta
## 0.12583692 0.04119796 0.02594432
##
## Description
## Tue Mar 09 14:10:32 2021
2.4 Estimação semi-paramétrica do índice caudal
2.4.1 Ajuste da GEV
par(mfrow=c(1,1))
fit <- fevd(m,type="GEV") ; fit #positive shape estimate but fairly large standard error##
## fevd(x = m, type = "GEV")
##
## [1] "Estimation Method used: MLE"
##
##
## Negative Log-Likelihood Value: -49.30563
##
##
## Estimated parameters:
## location scale shape
## 0.04047789 0.02300109 0.24408832
##
## Standard Error Estimates:
## location scale shape
## 0.005747584 0.004745185 0.242146504
##
## Estimated parameter covariance matrix.
## location scale shape
## location 3.303473e-05 1.799104e-05 -0.0006339456
## scale 1.799104e-05 2.251678e-05 -0.0003433678
## shape -6.339456e-04 -3.433678e-04 0.0586349296
##
## AIC = -92.61127
##
## BIC = -89.0771plot(fit) #The fit looks reasonable
ci(fit,type="parameter") #As expected the 95% confidence interval includes negative values## fevd(x = m, type = "GEV")
##
## [1] "Normal Approx."
##
## 95% lower CI Estimate 95% upper CI
## location 0.02921283 0.04047789 0.05174295
## scale 0.01370069 0.02300109 0.03230148
## shape -0.23051011 0.24408832 0.71868674return.level(fit,do.ci=T)## fevd(x = m, type = "GEV")
##
## [1] "Normal Approx."
##
## 95% lower CI Estimate 95% upper CI
## 2-year return level 0.035747834 0.04929667 0.06284551
## 20-year return level 0.063909413 0.14080836 0.21770730
## 100-year return level -0.002328087 0.23587671 0.474081502.4.2 Estimador de Pickands
amost<-PBR1$retorno
xord<-sort(amost)
n <- length(amost)
gamap<-numeric(n/4)
for(k in 1:(n/4)){
gamap[k]<-(1/log(2))*log( (xord[n-k+1]-xord[n-2*k+1]) /
(xord[n-2*k+1]-xord[n-4*k+1]) )}
plot(gamap,main='Estimador de Pickands',cex.lab=1.5,xlab="k",ylab='gamap', type="l")
2.4.3 Estimadores de Hill
par(mfrow=c(1,2))
## Estimador de hill construido manualmente
n<-length(xord) ; alphah<-R<-numeric(n-1)
(n-1)/2## [1] 125 for (k in 1:(n-1)/2){
for (i in 1:k){
R[i] <-xord[n-i+1]/xord[n-k]
}
R<-R[!is.na(R) & !is.infinite(R)]
alphah[k]<-(1/k)*sum( log(R[R>0]) )
}
plot(alphah,xlim=c(15,120),ylim=c(0,3.5),main='Estimador de Hill',cex.lab=1.5,ylab='Gama',xlab='k',type="l")
## Estimador de Hill gerado no R
hillPlot(xord,plottype = "xi")
2.4.4 Estimador de Deckers e Haan
gamad<-numeric(n/2)
n <- length(amost) # definido anteriormente
for(k in 1:(n/2)){
h1<- (1/k)*sum(log(xord[(n-1):(n-(k+1)+1)]/xord[n-(k+1)]))
h2<- (1/k)*sum(log((xord[(n-1):(n-(k+1)+1)]/xord[n-(k+1)])^2))
gamad[k]<-1+h1+(1/2)*(h1/h2-1)}
plot(gamad,main='Estimador de Deckers e Haan',type="l")
2.4.5 Gráficos dos estimadores semiparamétricos
#par(mfcol = c(3,2)) # shaparmPlot(xord) ---- Verificar!
par(mfrow=c(1,3))
plot(gamap,main='Pickands',ylab='',xlab='k',type="l")
plot(alphah,main='Hill',xlim=c(15,120),ylab='',xlab='k',type="l")
plot(gamad,main='Deckers e Haan',ylab='',xlab='k',type="l")