A seguinte função pode ser facilmente programada no R para o cálculo do VaR:
RMVaR<-function(rt, lambda = 0.94, k = 1, long = TRUE, alpha = 0.95,range.optim=seq(0.9,1,0.001),plot=FALSE){
#### rt é a série de log-retornos
#### se long == TRUE, entao estamos no caso de uma posicao comprada. Se long == FALSE, entao a posicao será considerada vendida
#### lambda é o parâmetro de suavização
#### 'alpha' é a probabilidade fixada da distribuição Normal padrao - Lembrete: o RiskMetrics supõe erros normais.
#### k é o horizonte para o qual se deseja calcular o VaR
#### Se lambda == NULL, entao o programa executa a rotina para todos os valores em range.optim=seq(0.9,1,0.001)
#### Se lambda == NULL e plot == TRUE, entao o programa imprime a soma de quadrado dos resíduos para os lambda em range.optim
# Se long=T, entao estou em uma posicao comprada e o VaR é definido na cauda inferior, i.e., lower.tail=T.
# Se long=F, entao o VaR é calculado na cauda superior e lower.tail = F.
Z=qnorm(0.05,lower.tail=long)
r2t<-rt^2
T=length(r2t)
SQR<-numeric()
VaR<-numeric()
sigma2hat<-numeric()
sigma2hat[1]<-r2t[1]
if(!is.null(lambda)){ # Entao nao ira para a otimizacao de lambda
for(t in 2:T){
sigma2hat[t]<-sigma2hat[t-1]*lambda+(1-lambda)*r2t[t-1]
#print(paste("hello 24-",t))
}
lambda_ = lambda
S2t1<-sigma2hat[T]*lambda+(1-lambda)*r2t[T]
VaR<-sqrt(S2t1*k)*Z
SQR<-sum((sigma2hat-r2t)^2,na.rm = TRUE)
}
else {
if(!is.null(range.optim)){
for( i in 1: length(range.optim)){
tmp<-RMVaR(rt, lambda = range.optim[i], k = k, long = long, alpha =alpha)
VaR[i]<-tmp[[1]]
SQR[i]<-tmp[[2]]
}
if(plot==TRUE){
plot(range.optim,SQR, type="l", col="darkblue",lwd=2, main="RiskMetrics\n Sum of Squared Residuals",xlab="lambda",ylab="SQR")
abline(h=min(SQR),v=range.optim[SQR==min(SQR)],lty=2,col="gray")
}
}
else {
return("Error! You must specify 'lambda' or 'range.optim' ")
}
lambda=range.optim[SQR==min(SQR)]
lambda_=range.optim
print(paste("RiskMetrics - The Minimum SQR is: ",round(min(SQR),4)," at lambda =",lambda))
}
return(list(VaR=VaR,SQR=SQR,lambda=lambda_))
}
Exemplo de Utilização utilizando dados e log-retornos do Ibovespa:
load("VaR.RData")
RMVaR(r_ibov,long=F,k=30) #Cálculo do VaR para periodo de 30 dias com o lambda padrão, sugerido pelo RiskMetrics, i.e., lambda=0.94
## $VaR
## [1] 0.19625
##
## $SQR
## [1] 0.014213
##
## $lambda
## [1] 0.94
RMVaR(r_ibov,long=F,k=30,lambda=0.88) #Cálculo do VaR para periodo de 30 dias com lambda = 0.88
## $VaR
## [1] 0.20049
##
## $SQR
## [1] 0.014179
##
## $lambda
## [1] 0.88
Optim<-RMVaR(r_ibov,long=F,k=30,lambda=NULL,range.optim=seq(0.8,1,0.001),plot=T) # Cálculo do VaR para periodo de 30 dias com lambda selecionado iterativamente de acordo com a menor soma de quadrados dos residuos
## [1] "RiskMetrics - The Minimum SQR is: 0.0142 at lambda = 0.901"
qtEmpirico<-function(x,p,long=TRUE,plot=T,nbreaks=30,print=TRUE,title="Empirical VaR"){
dnst<-density(x)
qt<-quantile(x,p*long+(1-p)*!long)
if(plot==T){
hist(x,breaks=nbreaks,col="lightgrey",border="white",freq=F,
ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),main=title)
par(new=T)
if(long != TRUE){ ### Posicao Vendida
plot(dnst$x[dnst$x>qt],
dnst$y[dnst$x>qt],type="h",
ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),
xaxt="n",yaxt="n",lwd=1,col="blue",ylab="",xlab="",main="")
}
else{ ### Posicao comprada
plot(dnst$x[dnst$x<qt],
dnst$y[dnst$x<qt],type="h",
ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),
xaxt="n",yaxt="n",lwd=0,col="blue",ylab="",xlab="",main="")
}
par(new=T)
plot(dnst,xaxt="n",yaxt="n",
ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),
lwd=1,col="darkblue",ylab="",xlab="",main="")
}
if(print==TRUE){
print(paste0("Probability: ",p))
print(paste0("Position long?: ",long))
print(paste0("Empirical Quantile: ",round(qt,5)))
}
return(qt)
}
qtEmpirico(rnorm(500),0.05,long=T)
## [1] "Probability: 0.05"
## [1] "Position long?: TRUE"
## [1] "Empirical Quantile: -1.75689"
## 5%
## -1.7569
O Pacote evir é análogo do R ao pacote evis do S-plus Assim, para estimar os parâmetros da distribuição de Valores extremos:
\[ \begin{align} VaR & = \mu_n - \frac{\sigma_n}{\xi_n}\left\{1-\left(-n log(1-p)\right)^{\xi_n}\right\} , \xi_n\neq 0 \\ & = \mu_n - \sigma_n log\left(-n log(1-p)\right) , \xi_n = 0 \\ VaR(k) & = k^{\xi}VaR \end{align} \]
utilizamos o comando gev:
gev(data,block)
data é o objeto contendo os dados que se deseja estimar a distribuição;block é o tamanho da janela utilizada para se calcular os máximos locais.Pode-se ainda utilizar de forma adaptada a função qtEmpirico(.) apresentada anteriormente para se calcular o VaR via a distribuição de valores extremos estimada aqui.