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)