sampleES

Algorithm in R to compute the BEST possible ES for the sum of d random losses identically distributed as F as described in Table 1 in

Puccetti, G. (2013). Sharp bounds on the expected shortfall for a sum of dependent random variables. Stat. Probabil. Lett., 83(4), 1227-1232 paper - preprint

Disclaimer

The sample code below does not use a randomized starting matrix and will produce slightly different than those in the reference above; see the remark "Randomise or not to randomise" in: 

Aas, K. and G. Puccetti (2014). Bounds for total economic capital: the DNB case study. Extremes, 17(4), 693-715 paper  SSRN.

# STEP 1.  CREATION OF THE MATRIX X

# epsilon = desired accuracy

epsilon=0.0001

# alpha = confidence level at which VaR is computed

alpha=0.999

# N = cardinality of the support of the approximation of EACH of the marginal distribution (the larger N, the more accurate, the more time expensive the algorithm)

N=2*10^6

# d = dimensionality of the risk portfolio 

d=3

# Q = quantile function of the marginal distribution F (this can be changed according to the end-user needs). 

#PARETO DISTRIBUTION

omega=4;

Q=function(x){(1-x)^(-1/omega)-1};

#LOGNORMAL DISTRIBUTION

#Q=function(x){qlnorm(x, meanlog = 6.474105, sdlog = 0.7213475, lower.tail = TRUE, log.p = FALSE)}

#EXPONENTIAL DISTRIBUTION

#Q=function(x){qexp(x,rate=expmu,log=FALSE,lower.tail=TRUE)}

# X=auxiliary matrix

X=matrix(rep(Q(seq(0,1,length=N+1)[1:N]),d),nrow=N)

# STEP 2-3.  ITERATIVE REARRANGEMENT

# ES = recursive value for the best possible ES

ES=Inf

# delta = difference between two consecutive estimates 

delta=Inf

# iteratively rearrange each column of the matrix until the difference between

# two consecutive estimates of the best ES is less than epsilon

while(delta >epsilon ){

#EStemp = auxiliary variable

EStemp=ES

for(j in 1:d) {

# rearrangement procedure

X[,j]=sort(X[,j],decreasing=TRUE)[rank(rowSums(X[,c(1:d)[-j]]))] }

# computation of the best ES in one iteration

#Y auxiliary variable representing the vector of the row-wise sums of X in increasing order

Y=sort(rowSums(X))

ES =sum(Y[(floor(N*alpha)+1):N])/N/(1-alpha)

delta=abs(ES-EStemp) 

}

# OUTPUT

# Output=worst possible ES at the fixed level of accuracy

print(ES)