Estimation of Value at Risk (VaR) and Conditional Value at Risk (CVaR)

Normal distribution

Normal.VaR.CVaR=function(sdX=1, n=10^3, B=10^3, pu=0.9, pk=0.95)
{
  u <- qnorm(pu,sd=sdX)
  X <- rnorm(n,0,sdX)
  Xu <- X[X>u]-u
  MXu=max(Xu)
  m <- length(Xu)
  
  # Real value
  VaRT <- sdX*qnorm(pk)
  CVaRT <- sdX*dnorm(qnorm(pk))/(1-pk)
  
  ## Empiric
  VaRE <- quantile(X,pk)[[1]]
  CVaRE <- mean(X[X>=VaRE])
  
  ## Parametric
  VaRP <- mean(X)+qnorm(pk)*sd(X)
  CVaRP <- -mean(X)+sd(X)*dnorm(qnorm(1-pk))/(1-pk)
  
  ## Bootstrap
  MR <- function(x){
    rX <- sample(X,n,rep=T)
    VaRB <- quantile(rX,pk)[[1]]
    CVaRB <- mean(rX[rX>=VaRB])
    c(VaRB,CVaRB)
  }
  EstR <- apply(replicate(B,MR(sample(X,n,rep=T))),1,mean)
  VaRB <- 2*VaRE - EstR[1]
  CVaRB <- 2*CVaRE - EstR[2]
  
  ## MCMC
  pcola <- 1-(1-pk)/(1-pu)
  nburn <- 3000
  nthin <- 50
  ndraw <- 10000
  
  kMHi <- 0.15 
  deltaMHi <- sdX/kMHi
  a1 <- 0.1; b1 <- 0.1; sdk <- 0.1
  a2 <- 0.1; b2 <- 0.1; sddelta <- 0.1
  a3 <- 1; b3 <- 1
  mu_xiIBD <- -0.7+0.61*pu
  sd_xiIBD <- 0.03
  mu_sigmaIBD <- (0.34+3.18*(1-pu)-12.4*(1-pu)^2)
  sd_sigmaIBD <- exp(-41.58+83.55*pu-46.24*pu^2)
  xiIBDi <- rnorm(1,mu_xiIBD,sd_xiIBD); sdxi <- 0.03
  sigmaIBDi <- rnorm(1,mu_sigmaIBD*sd(X),sd_sigmaIBD)
  while(MXu>(-sigmaIBDi/xiIBDi)||sigmaIBDi<0||xiIBDi>0)
    {
    xiIBDi <- rnorm(1,mu_xiIBD,sd_xiIBD)
    sigmaIBDi <- rnorm(1,mu_sigmaIBD*sd(X),sd_sigmaIBD)
    }
  sdsigma <- 1
  for(i in 1:nburn)
    {
    #MH
    kMHcan <- rnorm(1,kMHi,sdk)
    while(kMHcan<0) kMHcan <- rnorm(1,kMHi,sdk)
    deltaMHcan <- rnorm(1,deltaMHi,sddelta)
    while(deltaMHcan<MXu) deltaMHcan <- rnorm(1,deltaMHi,sddelta)
    pkMH <- (a1-m-1)*log(kMHcan/kMHi)+b1*(kMHi-kMHcan) + 
      (1/kMHcan-1/kMHi)*sum(log(1-Xu/deltaMHi))
    pdeltaMH <- (a2-m-1)*log(deltaMHcan/deltaMHi) +           
      b2*(deltaMHi-deltaMHcan) + 
      (1/kMHi-1)*sum(log(1-Xu/deltaMHcan)) - 
      (1/kMHi-1)*sum(log(1-Xu/deltaMHi))
    v <- log(runif(2))
    if(v[1]<pkMH) kMHi <- kMHcan
    if(v[2]<pdeltaMH) deltaMHi <- deltaMHcan
    xiMHi <- -kMHi
    sigmaMHi <- kMHi*deltaMHi
    VaRMH <- sigmaMHi/xiMHi*((1-pcola)^(-xiMHi)-1) + u
    CVaRMH <- sigmaMHi/xiMHi*(((1-pcola)^(-xiMHi))/(1-xiMHi)-1) + u
    
    #BD
    sdBDi <- sqrt(1/rgamma(1,a3+n/2,b2+sum(X^2)/2))
    
    VaRBD <- qnorm(pk)*sdBDi
    CVaRBD <- sdBDi*dnorm(qnorm(pk))/(1-pk)
    
    #IBD
    xiIBDcan <- rnorm(1,xiIBDi,sdxi)
    sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
    while(MXu>(-sigmaIBDcan/xiIBDcan)||sigmaIBDcan<0||xiIBDcan>0)
    {
      xiIBDcan <- rnorm(1,xiIBDi,sdxi)
      sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
    }
    mu_sigma <- mu_sigmaIBD*sdBDi
    pIBD <- m*log(sigmaIBDi/sigmaIBDcan) -
      0.5/sd_xiIBD^2 * ((xiIBDcan-mu_xiIBD)^2 - (xiIBDi-mu_xiIBD)^2) -
      0.5/sd_sigmaIBD^2 * ((sigmaIBDcan-mu_sigma)^2 - (sigmaIBDi - mu_sigma)^2) -
      (1+1/xiIBDcan)*sum(log(1+xiIBDcan*Xu/sigmaIBDcan)) +
      (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
    if(log(runif(1))<pIBD)
      {
      xiIBDi <- xiIBDcan
      sigmaIBDi <- sigmaIBDcan
      }
    VaRIBD <- sigmaIBDi/xiIBDi*((1-pcola)^(-xiIBDi)-1) + u
    CVaRIBD <- sigmaIBDi/xiIBDi*(((1-pcola)^(-xiIBDi))/(1-xiIBDi)-1) + u
  }
  draw <- c(xiMHi, sigmaMHi, VaRMH, CVaRMH, 
            sdBDi, VaRBD, CVaRBD,
            xiIBDi, sigmaIBDi, VaRIBD, CVaRIBD)
  
  for(i in 1:ndraw){
    for(j in 1:nthin)
      {
      #MH
      kMHcan <- rnorm(1,kMHi,sdk)
      while(kMHcan<0) kMHcan <- rnorm(1,kMHi,sdk)
      deltaMHcan <- rnorm(1,deltaMHi,sddelta)
      while(deltaMHcan<MXu) deltaMHcan <- rnorm(1,deltaMHi,sddelta)
      pkMH <- (a1-m-1)*log(kMHcan/kMHi)+b1*(kMHi-kMHcan) + 
        (1/kMHcan-1/kMHi)*sum(log(1-Xu/deltaMHi))
      pdeltaMH <- (a2-m-1)*log(deltaMHcan/deltaMHi)+b2*(deltaMHi-deltaMHcan) + 
        (1/kMHi-1)*sum(log(1-Xu/deltaMHcan))-(1/kMHi-1)*sum(log(1-Xu/deltaMHi))
      v <- log(runif(2))
      if(v[1]<pkMH) kMHi <- kMHcan
      if(v[2]<pdeltaMH) deltaMHi <- deltaMHcan
      xiMHi <- -kMHi; sigmaMHi <- kMHi*deltaMHi
      VaRMH <- sigmaMHi/xiMHi*((1-pcola)^(-xiMHi)-1) + u
      CVaRMH <- sigmaMHi/xiMHi*(((1-pcola)^(-xiMHi))/(1-xiMHi)-1) + u
      
      #BD
      sdBDi <- sqrt(1/rgamma(1,a3+n/2,b2+sum(X^2)/2))
      VaRBD <- qnorm(pk)*sdBDi
      CVaRBD <- sdBDi*dnorm(qnorm(pk))/(1-pk)
      
      #IBD
      xiIBDcan <- rnorm(1,xiIBDi,sdxi)
      sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
      while(MXu>=(-sigmaIBDcan/xiIBDcan)||sigmaIBDcan<=0||xiIBDcan>=0){
        xiIBDcan <- rnorm(1,xiIBDi,sdxi)
        sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
      }
      mu_sigma <- mu_sigmaIBD*sdBDi
      pIBD <- m*log(sigmaIBDi/sigmaIBDcan) -
        0.5/sd_xiIBD^2 * ((xiIBDcan-mu_xiIBD)^2 - (xiIBDi-mu_xiIBD)^2) -
        0.5/sd_sigmaIBD^2 * ((sigmaIBDcan-mu_sigma)^2 - (sigmaIBDi - mu_sigma)^2) -
        (1+1/xiIBDcan)*sum(log(1+xiIBDcan*Xu/sigmaIBDcan)) +
        (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
      if(log(runif(1))<pIBD){
        xiIBDi <- xiIBDcan
        sigmaIBDi <- sigmaIBDcan
      }
      VaRIBD <- sigmaIBDi/xiIBDi*(-1+(1-pcola)^(-xiIBDi)) + u 
      CVaRIBD <- sigmaIBDi/xiIBDi*(((1-pcola)^(-xiIBDi))/(1-xiIBDi)-1) + u  
    }
    draw <- c(draw,
              c(xiMHi, sigmaMHi, VaRMH, CVaRMH,
                sdBDi, VaRBD, CVaRBD,
                xiIBDi, sigmaIBDi, VaRIBD, CVaRIBD))
  }
  
  draw <- matrix(draw, ncol=11, byrow = TRUE)
  list(standar.deviation=sdX, length=n, threshold=pu, p=pk, 
    real.value.VaR=VaRT, real.value.CVaR=CVaRT,
    empirical.estimation.VaR=VaRE, empirical.estimation.CVaR=CVaRE,
    bootstrap.VaR=VaRB, bootstrap.CVaR=CVaRB,
    parametric.VaR=VaRP, parametric.CVaR=CVaRP,
    EVB.VaR=mean(draw[,3]), EVB.CVaR=mean(draw[,4]),
    parametric.Bayesian.VaR=mean(draw[,6]), parametric.Bayesian.CVaR=mean(draw[,7]),
    IPB.VaR=mean(draw[,10]), IPB.CVaR=mean(draw[,11]))
}

For example, to estimate the \(VaR(0.95)\) and \(CVaR(0.95)\), from a sample of size \(100\) with a standard normal distribution, one can execute the function:

Normal.VaR.CVaR(n=10^2)

## $standar.deviation
## [1] 1
## 
## $length
## [1] 100
## 
## $threshold
## [1] 0.9
## 
## $p
## [1] 0.95
## 
## $real.value.VaR
## [1] 1.644854
## 
## $real.value.CVaR
## [1] 2.062713
## 
## $empirical.estimation.VaR
## [1] 1.472978
## 
## $empirical.estimation.CVaR
## [1] 1.89999
## 
## $bootstrap.VaR
## [1] 1.455893
## 
## $bootstrap.CVaR
## [1] 1.965748
## 
## $parametric.VaR
## [1] 1.700349
## 
## $parametric.CVaR
## [1] 2.079118
## 
## $EVB.VaR
## [1] 1.630196
## 
## $EVB.CVaR
## [1] 2.013311
## 
## $parametric.Bayesian.VaR
## [1] 1.666037
## 
## $parametric.Bayesian.CVaR
## [1] 2.089277
## 
## $IPB.VaR
## [1] 1.63755
## 
## $IPB.CVaR
## [1] 2.062008

Exponential distribution

Exp.VaR.CVaR=function(lambdaX=1, n=10^3, B=10^3, pu=0.9, pk=0.95)
{
  u <- qexp(pu,lambdaX)
  X <- rexp(n,lambdaX)
  Xu <- X[X>u]-u
  MXu=max(Xu)
  m <- length(Xu)
  
  # Real value
  VaRT <- -log(1-pk)/lambdaX
  CVaRT <- (1-log(1-pk))/lambdaX
  
  ## Empiric
  VaRE <- quantile(X,pk)[[1]]
  CVaRE <- mean(X[X>=VaRE])
  
  ## Parametric
  VaRP <- -log(1-pk)*mean(X)
  CVaRP <- mean(X)*(1-log(1-pk))
  
  ## Bootstrap
  MR <- function(x){
    rX <- sample(X,n,rep=T)
    VaRB <- quantile(rX,pk)[[1]]
    CVaRB <- mean(rX[rX>=VaRB])
    c(VaRB,CVaRB)
  }
  EstR <- apply(replicate(B,MR(sample(X,n,rep=T))),1,mean)
  VaRB <- 2*VaRE - EstR[1]
  CVaRB <- 2*CVaRE - EstR[2]
  
  ## MCMC
  pcola <- 1-(1-pk)/(1-pu)
  nburn <- 3000
  nthin <- 50
  ndraw <- 10000
  
  a1 <- 1; b1 <- 1 
  a2 <- 1; b2 <- 1
  sigmaIBDi <- mean(X)
  for(i in 1:nburn)
    {
    #MH
    sigmaMHi <- 1/rgamma(1,a1+m,b1+sum(Xu))
    VaRMH <- -log(1-pcola)*sigmaMHi + u 
    CVaRMH <- (1-log(1-pcola))*sigmaMHi + u 
                      
    #BD
    lambdaBDi <- rgamma(1,a2+n,b2+sum(X))
    VaRBD <- -log(1-pk)/lambdaBDi
    CVaRBD <- (1-log(1-pk))/lambdaBDi
                      
    #IBD
    sigmaIBDcan <- 1/lambdaBDi
    pIBD <- m*log(sigmaIBDi/sigmaIBDcan) + sum(Xu)*(1/sigmaIBDi-1/sigmaIBDcan)
    if(log(runif(1))<pIBD) sigmaIBDi <- sigmaIBDcan
    VaRIBD <- -log(1-pcola)*sigmaIBDi + u
    CVaRIBD <- (1-log(1-pcola))*sigmaIBDi + u 
    }
  draw <- c(sigmaMHi, VaRMH, CVaRMH,
            lambdaBDi, VaRBD, CVaRBD,
            sigmaIBDi, VaRIBD, CVaRIBD)
                    
  for(i in 1:ndraw)
    {
    for(j in 1:nthin)
      {
      #MH
      sigmaMHi <- 1/rgamma(1,a1+m,b1+sum(Xu))
      VaRMH <- -log(1-pcola)*sigmaMHi + u
      CVaRMH <- (1-log(1-pcola))*sigmaMHi + u 
      #BD
      lambdaBDi <- rgamma(1,a2+n,b2+sum(X))
      VaRBD <- -log(1-pk)/lambdaBDi
      CVaRBD <- (1-log(1-pk))/lambdaBDi
      
      #IBD
      sigmaIBDcan <- 1/lambdaBDi
      pIBD <- m*log(sigmaIBDi/sigmaIBDcan) +sum(Xu)*(1/sigmaIBDi-1/sigmaIBDcan)
      if(log(runif(1))<pIBD) sigmaIBDi <- sigmaIBDcan
      VaRIBD <- -log(1-pcola)*sigmaIBDi + u
      CVaRIBD <- (1-log(1-pcola))*sigmaIBDi + u 
      }
    draw <- c(draw,c(sigmaMHi, VaRMH, CVaRMH,
                     lambdaBDi, VaRBD, CVaRBD,
                     sigmaIBDi, VaRIBD, CVaRIBD))
                    }
  draw <- matrix(draw,ncol=9,byrow = TRUE)
  
  list(lambda=lambdaX, length=n, threshold=pu, p=pk, 
    real.value.VaR=VaRT, real.value.CVaR=CVaRT,
    empirical.estimation.VaR=VaRE, empirical.estimation.CVaR=CVaRE,
    bootstrap.VaR=VaRB, bootstrap.CVaR=CVaRB,
    parametric.VaR=VaRP, parametric.CVaR=CVaRP,
    EVB.VaR=mean(draw[,2]), EVB.CVaR=mean(draw[,3]),
    parametric.Bayesian.VaR=mean(draw[,5]), parametric.Bayesian.CVaR=mean(draw[,6]),
    IPB.VaR=mean(draw[,8]), IPB.CVaR=mean(draw[,9]))
}

Exp.VaR.CVaR()

## $lambda
## [1] 1
## 
## $length
## [1] 1000
## 
## $threshold
## [1] 0.9
## 
## $p
## [1] 0.95
## 
## $real.value.VaR
## [1] 2.995732
## 
## $real.value.CVaR
## [1] 3.995732
## 
## $empirical.estimation.VaR
## [1] 3.007078
## 
## $empirical.estimation.CVaR
## [1] 3.880637
## 
## $bootstrap.VaR
## [1] 3.013336
## 
## $bootstrap.CVaR
## [1] 3.896095
## 
## $parametric.VaR
## [1] 2.92424
## 
## $parametric.CVaR
## [1] 3.900376
## 
## $EVB.VaR
## [1] 2.970319
## 
## $EVB.CVaR
## [1] 3.933656
## 
## $parametric.Bayesian.VaR
## [1] 2.928706
## 
## $parametric.Bayesian.CVaR
## [1] 3.906332
## 
## $IPB.VaR
## [1] 2.978757
## 
## $IPB.CVaR
## [1] 3.954268

Cauchy Distribution

Cauchy.VaR.CVaR=function(deltaX=1, n=10^3, B=10^3, pu=0.9, pk=0.95)
{
  u <- qcauchy(pu,0,deltaX)
  X <- rcauchy(n,0,deltaX)
  Xu <- X[X>u]-u
  MXu=max(Xu)
  m <- length(Xu)
  
  # Real value
  VaRT <- deltaX*tan(pi*(pk-0.5))
  
  ## Empiric
  VaRE <- quantile(X,pk)[[1]]
  
  ## Parametric
  VaRP <- IQR(X)*tan(pi*(pk-0.5))/2
  
  ## Bootstrap
  B <- 1000
  MR <- function(x)
    {
    rX <- sample(X,n,rep=T)
    quantile(rX,pk)[[1]]
    }
  EstR <- mean(replicate(B,MR(sample(X,n,rep=T))))
  VaRB <- 2*VaRE - EstR
  
  ## MCMC
  pcola <- 1-(1-pk)/(1-pu)
  nburn <- 3000
  nthin <- 50
  ndraw <- 10000
  
  xiMHi <- 1
  a1 <- 10^(-6); b1 <- 10^(-3)
  sd_xiMH <- 0.01
  a2 <- 1; b2 <- 1
  sigmaMHi <- 1/rgamma(1,a2,b2); sd_sigmaMH <- 1
  deltaBDi <- as.numeric((quantile(X,3/4)-quantile(X,1/4))/2)
  a3 <- 1; b3 <- 1; sd_deltaBD <- 1
  mu_xiIBD <- 1; sd_xiIBD <- 0.065
  mu_sigmaIBD <- 1/(pi*(1-pu))
  sd_sigmaIBD <- exp(266.13-588.51*pu+323.57*pu^2)
  xiIBDi <- rnorm(1,mu_xiIBD,sd_xiIBD); 
  sdxi <- 0.1
  while(xiIBDi<0) xiIBDi <- rnorm(1,mu_xiIBD,sd_xiIBD)
  sigmaIBDi <- rnorm(1,mu_sigmaIBD*deltaBDi,sd_sigmaIBD)
  while(sigmaIBDi<0) sigmaIBDi <- rnorm(1,mu_sigmaIBD*deltaBDi,sd_sigmaIBD)
  sdsigma <- 0.1
  for(i in 1:nburn)
    {
    #MH
    xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
    while(xiMHcan<b1) xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
    sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH)
    while(sigmaMHcan<0) sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH) 
    pxiMH <- (a1+1)*log(xiMHi/xiMHcan) + (1+1/xiMHi)*sum(log(1+xiMHi*Xu/sigmaMHi)) - 
      (1+1/xiMHcan)*sum(log(1+xiMHcan*Xu/sigmaMHi))
    psigmaMH <- (m+a2+1)*log(sigmaMHi/sigmaMHcan) + 
      b2*(1/sigmaMHi - 1/sigmaMHcan) + 
      (1+1/xiMHi)*sum(log(1+xiMHi*Xu/sigmaMHi)) -
      (1+1/xiMHi)*sum(log(1+xiMHi*Xu/sigmaMHcan))
    v <- log(runif(2))
    if(v[1]<pxiMH) xiMHi <- xiMHcan
    if(v[2]<psigmaMH) sigmaMHi <- sigmaMHcan
    VaRMH <- sigmaMHi/xiMHi*((1-pcola)^(-xiMHi)-1) + u
    
    #BD
    deltaBDcan <- rnorm(1,deltaBDi,sd_deltaBD)
    while(deltaBDcan<0) deltaBDcan <- rnorm(1,deltaBDi,sd_deltaBD)
    pBD <- (a3-n-1)*log(deltaBDcan/deltaBDi) - b3*(deltaBDcan-deltaBDi) - 
      sum(log(1+X^2/deltaBDcan^2)) + sum(log(1+X^2/deltaBDi^2))
    if(log(runif(1))<pBD) deltaBDi <- deltaBDcan
    VaRBD <- deltaBDi*tan(pi*(pk-0.5))
    
    #IBD
    xiIBDcan <- rnorm(1,xiIBDi,sdxi)
    while(xiIBDcan<0) xiIBDcan <- rnorm(1,xiIBDi,sdxi)
    sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
    while(sigmaIBDcan<0) sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
    mu_sigma <- mu_sigmaIBD*deltaBDi
    pxiIBD <- - 0.5/sd_xiIBD^2 * ((xiIBDcan-mu_xiIBD)^2 - (xiIBDi-mu_xiIBD)^2) -
      (1+1/xiIBDcan)*sum(log(1+xiIBDcan*Xu/sigmaIBDi)) + 
      (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
    psigmaIBD <- m*log(sigmaIBDi/sigmaIBDcan) - 
      0.5/sd_sigmaIBD^2 * ((sigmaIBDcan-mu_sigma)^2 - (sigmaIBDi-mu_sigma)^2) -
      (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDcan)) +
      (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
    v <- log(runif(2))
    if(v[1]<pxiIBD) xiIBDi <- xiIBDcan
    if(v[2]<psigmaIBD) sigmaIBDi <- sigmaIBDcan
    VaRIBD <- sigmaIBDi/xiIBDi*((1-pcola)^(-xiIBDi)-1) + u
    }
  draw <- c(xiMHi, sigmaMHi, VaRMH,
            deltaBDi, VaRBD,
            xiIBDi, sigmaIBDi, VaRIBD)
  for(i in 1:ndraw)
    {
    for(j in 1:nthin)
      {
      #MH
      xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
      while(xiMHcan<b1) xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
      sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH)
      while(sigmaMHcan<0) sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH)
      pxiMH <- (a1+1)*log(xiMHi/xiMHcan) + 
        (1+1/xiMHi)*sum(log(1+xiMHi*Xu/sigmaMHi)) - 
        (1+1/xiMHcan)*sum(log(1+xiMHcan*Xu/sigmaMHi))
      psigmaMH <- (m+a2+1)*log(sigmaMHi/sigmaMHcan) + 
        b2*(1/sigmaMHi - 1/sigmaMHcan) + 
        (1+1/xiMHi)*sum(log(1+xiMHi*Xu/sigmaMHi)) -
        (1+1/xiMHi)*sum(log(1+xiMHi*Xu/sigmaMHcan))
      v <- log(runif(2))
      if(v[1]<pxiMH) xiMHi <- xiMHcan
      if(v[2]<psigmaMH) sigmaMHi <- sigmaMHcan
      VaRMH <- sigmaMHi/xiMHi*((1-pcola)^(-xiMHi)-1) + u
      #BD
      deltaBDcan <- rnorm(1,deltaBDi,sd_deltaBD)
      while(deltaBDcan<0) deltaBDcan <- rnorm(1,deltaBDi,sd_deltaBD)
      pBD <- (a3-n-1)*log(deltaBDcan/deltaBDi) -
        b3*(deltaBDcan-deltaBDi) - 
        sum(log(1+X^2/deltaBDcan^2)) + 
        sum(log(1+X^2/deltaBDi^2))
      if(log(runif(1))<pBD) deltaBDi <- deltaBDcan
      VaRBD <- deltaBDi*tan(pi*(pk-0.5))
      
      #IBD
      xiIBDcan <- rnorm(1,xiIBDi,sdxi)
      while(xiIBDcan<0) xiIBDcan <- rnorm(1,xiIBDi,sdxi)
      sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
      while(sigmaIBDcan<0) sigmaIBDcan <- rnorm(1,sigmaIBDi,sdsigma)
      mu_sigma <- mu_sigmaIBD*deltaBDi
      pxiIBD <- - 0.5/sd_xiIBD^2 * ((xiIBDcan-mu_xiIBD)^2 - (xiIBDi-mu_xiIBD)^2) -
        (1+1/xiIBDcan)*sum(log(1+xiIBDcan*Xu/sigmaIBDi)) +
        (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
      psigmaIBD <- m*log(sigmaIBDi/sigmaIBDcan) - 
        0.5/sd_sigmaIBD^2 * ((sigmaIBDcan-mu_sigma)^2 - (sigmaIBDi-mu_sigma)^2) -
        (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDcan)) +
        (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
      v <- log(runif(2))
      if(v[1]<pxiIBD) xiIBDi <- xiIBDcan
      if(v[2]<psigmaIBD) sigmaIBDi <- sigmaIBDcan
      VaRIBD <- sigmaIBDi/xiIBDi*((1-pcola)^(-xiIBDi)-1) + u
      }
    draw <- c(draw, c(xiMHi, sigmaMHi, VaRMH,
                      deltaBDi, VaRBD,
                      xiIBDi, sigmaIBDi, VaRIBD))
    }
  draw <- matrix(draw,ncol=8,byrow = TRUE)
                    
  list(delta=deltaX, length=n, threshold=pu, p=pk, 
    real.value.VaR=VaRT, 
    empirical.estimation.VaR=VaRE, 
    bootstrap.VaR=VaRB, 
    parametric.VaR=VaRP, 
    EVB.VaR=mean(draw[,3]), 
    parametric.Bayesian.VaR=mean(draw[,5]), 
    IPB.VaR=mean(draw[,8]))
}

Cauchy.VaR.CVaR()

## $delta
## [1] 1
## 
## $length
## [1] 1000
## 
## $threshold
## [1] 0.9
## 
## $p
## [1] 0.95
## 
## $real.value.VaR
## [1] 6.313752
## 
## $empirical.estimation.VaR
## [1] 7.491547
## 
## $bootstrap.VaR
## [1] 7.486486
## 
## $parametric.VaR
## [1] 6.549596
## 
## $EVB.VaR
## [1] 6.253361
## 
## $parametric.Bayesian.VaR
## [1] 6.544996
## 
## $IPB.VaR
## [1] 6.363789

Gamma Distribution

Gamma.VaR.CVaR=function(alphaX=2, betaX=2, n=10^3, B=10^3, pu=0.9, pk=0.95)
{
  u <- qgamma(pu,alphaX,betaX)
  X <- rgamma(n,alphaX,betaX)
  Xu <- X[X>u]-u
  MXu=max(Xu)
  m <- length(Xu)
  
  # Real value
  VaRT <- qgamma(pk,alphaX,betaX)
  CVaRT <- 1/((1-pk)*gamma(alphaX)*betaX)*gamma_inc(alphaX+1,betaX*VaRT)
  
  ## Empiric
  VaRE <- quantile(X,pk)[[1]]
  CVaRE <- mean(X[X>=VaRE])
  
  ## Parametric
  library(gsl)
  GammaNewton <- function(x,alpha0)
    {
    m <- mean(x)
    mlog <- mean(log(x))
    alpha0*(1+(digamma(alpha0)+log(m/alpha0)-mlog)/(1-alpha0*trigamma(alpha0)))
    }
  emvalphaGamma <- function(x)
    {
    m <- mean(x)
    v <- var(x)
    alpha0 <- m^2/v
    error <- 1
    i <- 0
    while(error>1E-4)
      {
      i <- i+1
      alpha1 <- GammaNewton(x,alpha0)
      error <- abs(alpha1-alpha0); alpha0 <- alpha1
      }
    alpha1
    }
  alphaX1 <- emvalphaGamma(X); betaX1 <- alphaX1/mean(X)
  VaRP <- qgamma(pk,alphaX1,betaX1)
  CVaRP <- 1/((1-pk)*gamma(alphaX1)*betaX1)*gamma_inc(alphaX1+1,betaX1*VaRP)
  
  ## Bootstrap
  MR <- function(x){
    rX <- sample(X,n,rep=T)
    VaRB <- quantile(rX,pk)[[1]]
    CVaRB <- mean(rX[rX>=VaRB])
    c(VaRB,CVaRB)
  }
  EstR <- apply(replicate(B,MR(sample(X,n,rep=T))),1,mean)
  VaRB <- 2*VaRE - EstR[1]
  CVaRB <- 2*CVaRE - EstR[2]
  
  ## MCMC
  pcola <- 1-(1-pk)/(1-pu)
  nburn <- 3000
  nthin <- 50
  ndraw <- 10000
  
  xiMHi <- 0.1; sd_xiMH <- 0.01
  sigmaMHi <- ifelse(m>=2,var(Xu)/mean(Xu),1/betaX); sd_sigmaMH <- 1
  alphaBDi <- mean(X)^2/var(X); sd_alphaBD <- 0.5
  betaBDi <- mean(X)/var(X); sd_betaBD <- 0.5
  xiIBDi <- -0.032+0.014/alphaBDi
  sigmaIBDi <- 1/betaBDi*(0.5+0.5*sqrt(alphaBDi))
  for(i in 1:nburn){
    #MH
    xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
    sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH)
    mXu <- min(1+xiMHcan*Xu/sigmaMHcan)
    while(mXu<0||sigmaMHcan<0||xiMHcan<(-0.5))
      {
      xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
      sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH)
      mXu <- min(1+xiMHcan*Xu/sigmaMHcan)
      }
    pMH <- -(m+1)*log(sigmaMHcan) -
      (1+xiMHcan)/xiMHcan*sum(log(1+xiMHcan*Xu/sigmaMHcan)) - 
      log(1+xiMHcan) - 0.5*log(1+2*xiMHcan) + 
      (m+1)*log(sigmaMHi) + 
      (1+xiMHi)/xiMHi*sum(log(1+xiMHi*Xu/sigmaMHi)) + 
      log(1+xiMHi) + 0.5*log(1+2*xiMHi)
    if(log(runif(1))<pMH)
      {
      xiMHi <- xiMHcan
      sigmaMHi <- sigmaMHcan
      }
    VaRMH <- sigmaMHi/xiMHi*((1-pcola)^(-xiMHi)-1) + u
    CVaRMH <- sigmaMHi/xiMHi*(((1-pcola)^(-xiMHi))/(1-xiMHi)-1) + u
    
    #BD
    alphaBDcan <- rnorm(1,alphaBDi,sd_alphaBD)
    betaBDcan <- rnorm(1,betaBDi,sd_betaBD)
    while(alphaBDcan<0||betaBDcan<0)
      {
      alphaBDcan <- rnorm(1,alphaBDi,sd_alphaBD)
      betaBDcan <- rnorm(1,betaBDi,sd_betaBD)
      }
    pBD <- -n*log(gamma(alphaBDcan)) + 
      (n*alphaBDcan-1)*log(betaBDcan) + 
      (alphaBDcan-1)*sum(log(X)) - betaBDcan*sum(X) + 
      0.5*log(alphaBDcan*trigamma(alphaBDcan)-1) +
      n*log(gamma(alphaBDi)) - 
      (n*alphaBDi-1)*log(betaBDi) -
      (alphaBDi-1)*sum(log(X)) + betaBDi*sum(X) - 
      0.5*log(alphaBDi*trigamma(alphaBDi)-1) 
    if(log(runif(1))<pBD)
      {
      alphaBDi <- alphaBDcan 
      betaBDi <- betaBDcan
      }
    VaRBD <- qgamma(pk,alphaBDi,betaBDi)
    CVaRBD <- 1/((1-pk)*gamma(alphaBDi)*betaBDi)*gamma_inc(alphaBDi+1,betaBDi*VaRBD)
    
    #IBD
    xiIBDcan <- -0.032+0.014/alphaBDi
    sigmaIBDcan <- 1/betaBDi * (0.5+0.5*sqrt(alphaBDi))
    pIBD <- -m*log(sigmaIBDcan) -
      (1+1/xiIBDcan)*sum(log(1+xiIBDcan*Xu/sigmaIBDcan)) +
      m*log(sigmaIBDi) + 
      (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
    if(log(runif(1))<pIBD)
      {
      xiIBDi <- xiIBDcan
      sigmaIBDi <- sigmaIBDcan
      }
    VaRIBD <- sigmaIBDi/xiIBDi*((1-pcola)^(-xiIBDi)-1) + u
    CVaRIBD <- sigmaIBDi/xiIBDi*(((1-pcola)^(-xiIBDi))/(1-xiIBDi)-1) + u
    }
  draw <- c(xiMHi, sigmaMHi, VaRMH, CVaRMH,
            alphaBDi, betaBDi, VaRBD, CVaRBD, 
            xiIBDi, sigmaIBDi, VaRIBD, CVaRIBD)
  for(i in 1:ndraw)
    {
    for(j in 1:nthin)
      {
      #MH
      xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
      sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH)
      mXu <- min(1+xiMHcan*Xu/sigmaMHcan)
      while(mXu<0||sigmaMHcan<0||xiMHcan<(-0.5))
        {
        xiMHcan <- rnorm(1,xiMHi,sd_xiMH)
        sigmaMHcan <- rnorm(1,sigmaMHi,sd_sigmaMH)
        mXu <- min(1+xiMHcan*Xu/sigmaMHcan)
        }
      pMH <- -(m+1)*log(sigmaMHcan) -
        (1+xiMHcan)/xiMHcan*sum(log(1+xiMHcan*Xu/sigmaMHcan)) - 
        log(1+xiMHcan) - 0.5*log(1+2*xiMHcan) + 
        (m+1)*log(sigmaMHi) + 
        (1+xiMHi)/xiMHi*sum(log(1+xiMHi*Xu/sigmaMHi)) + 
        log(1+xiMHi) + 0.5*log(1+2*xiMHi)
      if(log(runif(1))<pMH)
        {
        xiMHi <- xiMHcan
        sigmaMHi <- sigmaMHcan
        }
      VaRMH <- sigmaMHi/xiMHi*((1-pcola)^(-xiMHi)-1) + u
      CVaRMH <- sigmaMHi/xiMHi*(((1-pcola)^(-xiMHi))/(1-xiMHi)-1) + u
      
      #BD
      alphaBDcan <- rnorm(1,alphaBDi,sd_alphaBD)
      betaBDcan <- rnorm(1,betaBDi,sd_betaBD)
      while(alphaBDcan<0||betaBDcan<0)
        {
        alphaBDcan <- rnorm(1,alphaBDi,sd_alphaBD)
        betaBDcan <- rnorm(1,betaBDi,sd_betaBD)
        }
      pBD <- -n*log(gamma(alphaBDcan)) + 
        (n*alphaBDcan-1)*log(betaBDcan) + 
        (alphaBDcan-1)*sum(log(X)) - betaBDcan*sum(X) + 
        0.5*log(alphaBDcan*trigamma(alphaBDcan)-1) +
        n*log(gamma(alphaBDi)) - 
        (n*alphaBDi-1)*log(betaBDi) -
        (alphaBDi-1)*sum(log(X)) + betaBDi*sum(X) - 
        0.5*log(alphaBDi*trigamma(alphaBDi)-1) 
      if(log(runif(1))<pBD)
        {
        alphaBDi <- alphaBDcan
        betaBDi <- betaBDcan
        }
      VaRBD <- qgamma(pk,alphaBDi,betaBDi)
      CVaRBD <- 1/((1-pk)*gamma(alphaBDi)*betaBDi)*gamma_inc(alphaBDi+1,betaBDi*VaRBD)
      
      #IBD
      xiIBDcan <- -0.032+0.014/alphaBDi
      sigmaIBDcan <- 1/betaBDi * (0.5+0.5*sqrt(alphaBDi))
      pIBD <- -m*log(sigmaIBDcan) -
        (1+1/xiIBDcan)*sum(log(1+xiIBDcan*Xu/sigmaIBDcan)) +
        m*log(sigmaIBDi) +
        (1+1/xiIBDi)*sum(log(1+xiIBDi*Xu/sigmaIBDi))
      if(log(runif(1))<pIBD)
        {
        xiIBDi <- xiIBDcan
        sigmaIBDi <- sigmaIBDcan
        }
      VaRIBD <- sigmaIBDi/xiIBDi*((1-pcola)^(-xiIBDi)-1) + u
      CVaRIBD <- sigmaIBDi/xiIBDi*(((1-pcola)^(-xiIBDi))/(1-xiIBDi)-1) + u
      }
    draw <- c(draw, c(xiMHi, sigmaMHi, VaRMH, CVaRMH,
                      alphaBDi, betaBDi, VaRBD, CVaRBD,
                      xiIBDi, sigmaIBDi, VaRIBD, CVaRIBD))
    }
  draw <- matrix(draw,ncol=12,byrow = TRUE)
  list(alpha=alphaX, beta=betaX, length=n, threshold=pu, p=pk, 
    real.value.VaR=VaRT, real.value.CVaR=CVaRT,
    empirical.estimation.VaR=VaRE, empirical.estimation.CVaR=CVaRE,
    bootstrap.VaR=VaRB, bootstrap.CVaR=CVaRB,
    parametric.VaR=VaRP, parametric.CVaR=CVaRP,
    EVB.VaR=mean(draw[,3]), EVB.CVaR=mean(draw[,4]),
    parametric.Bayesian.VaR=mean(draw[,7]), parametric.Bayesian.CVaR=mean(draw[,8]),
    IPB.VaR=mean(draw[,11]), IPB.CVaR=mean(draw[,12]))
}

library(gsl)
Gamma.VaR.CVaR()

## $alpha
## [1] 2
## 
## $beta
## [1] 2
## 
## $length
## [1] 1000
## 
## $threshold
## [1] 0.9
## 
## $p
## [1] 0.95
## 
## $real.value.VaR
## [1] 2.371932
## 
## $real.value.CVaR
## [1] 2.958982
## 
## $empirical.estimation.VaR
## [1] 2.344989
## 
## $empirical.estimation.CVaR
## [1] 2.805226
## 
## $bootstrap.VaR
## [1] 2.347253
## 
## $bootstrap.CVaR
## [1] 2.81593
## 
## $parametric.VaR
## [1] 2.391305
## 
## $parametric.CVaR
## [1] 2.986132
## 
## $EVB.VaR
## [1] 2.316178
## 
## $EVB.CVaR
## [1] 2.790418
## 
## $parametric.Bayesian.VaR
## [1] 2.394281
## 
## $parametric.Bayesian.CVaR
## [1] 2.990394
## 
## $IPB.VaR
## [1] 2.364712
## 
## $IPB.CVaR
## [1] 2.95062