Sample code for computing the figures in Table 4 in:
Embrechts, P., Puccetti, G. and L. Rüschendorf (2013). Model uncertainty and VaR aggregation. J. Bank. Financ., 37(8), 2750-2764. paper - post-print (updated) version
Disclaimer
The sample code below uses a randomized starting matrix. This is not really necessary (can be dropped), time consuming and will produce slightly different results at each run of the algorithm); 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.
#Algorithm in R to compute the WORST possible VaR range for the sum of d random losses identically distributed like F; see Table 4 in [EPR13] for F=Pareto.
# INPUT PARAMETERS
# 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=10^5
# d = dimensionality of the risk portfolio
d=8
# Q = quantile function of the marginal distribution F (this can be changed according to the end-user needs).
#PARETO DISTRIBUTION
omega=2;
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)}
# CREATION OF THE MATRICES LX AND UX
# X=auxiliary matrix
X=matrix(rep(Q(seq(alpha,1,length=N+1)),d),nrow=N+1)
#LX = matrix representing an approximation of the marginal distributions from below
LX=X[1:N,]
#UX = matrix representing an approximation of the marginal distributions from above
UX=X[2:(N+1),]
#randomization of LX and UX (this is not really necessary; see the Disclaimer above)
for(j in 1:d){
LX[,j]<-LX[rank(runif(N)),j]
UX[,j]<-UX[rank(runif(N)),j]}
#REARRANGEMENT PROCEDURE FOR LX AND UX
# epsilon = desired accuracy
epsilon=0.001
# VaR = recursive value for the worst possible VaR
LVaR=Inf
# delta = difference between two consecutive estimates
delta=Inf
# iteratively rearrange each column of the matrix LX until the difference between
# two consecutive estimates of the worst VaR is less than epsilon
while(delta >epsilon ){
#LVaRtemp = auxiliary variable
LVaRtemp=LVaR
for(j in 1:d) {
LX[,j]=sort(LX[,j],decreasing=TRUE)[rank(rowSums(LX[,c(1:d)[-j]]))] }
LVaR= min(rowSums(LX))
delta=abs(LVaR-LVaRtemp)
}
# iteratively rearrange each column of the matrix UX until the difference between
# two consecutive estimates of the worst VaR is less than epsilon
UVaR=Inf
delta2=Inf
while(delta2 >epsilon ){
UVaRtemp=UVaR
for(j in 1:d) {
UX[,j]=sort(UX[,j],decreasing=TRUE)[rank(rowSums(UX[,c(1:d)[-j]]))] }
UVaR= min(rowSums(UX))
delta2=abs(UVaR-UVaRtemp)
}
# OUTPUT
#WORST VAR is the interval [LVaR,UVaR].
print(LVaR)
print(UVaR)
#Algorithm in R to compute the BEST possible VaR range for the sum of d random losses identically distributed like F; see Table 4 in [EPR13] for the F=Pareto.
# INPUT PARAMETERS
# 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=10^5
# d = dimensionality of the risk portfolio
d=8
# Q = quantile function of the marginal distribution F (this can be changed according to the end-user needs).
#PARETO DISTRIBUTION
omega=2;
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)}
# CREATION OF THE MATRICES LX AND UX
# X=auxiliary matrix
X=matrix(rep(Q(seq(0,alpha,length=N+1)),d),nrow=N+1)
#LX = matrix representing an approximation of the marginal distributions from below
LX=X[1:N,]
#UX = matrix representing an approximation of the marginal distributions from above
UX=X[2:(N+1),]
#randomization of LX and UX (this is not really necessary; see the Disclaimer above)
for(j in 1:d){
LX[,j]<-LX[rank(runif(N)),j]
UX[,j]<-UX[rank(runif(N)),j]}
#REARRANGEMENT PROCEDURE FOR LX AND UX
# epsilon = desired accuracy
epsilon=0.001
# VaR = recursive value for the best possible VaR
LVaR=Inf
# delta = difference between two consecutive estimates
delta=Inf
# iteratively rearrange each column of the matrix LX until the difference between
# two consecutive estimates of the best VaR is less than epsilon
while(delta >epsilon ){
#LVaRtemp = auxiliary variable
LVaRtemp=LVaR
for(j in 1:d) {
LX[,j]=sort(LX[,j],decreasing=TRUE)[rank(rowSums(LX[,c(1:d)[-j]]))] }
LVaR= max(rowSums(LX))
delta=abs(LVaR-LVaRtemp)
}
# iteratively rearrange each column of the matrix UX until the difference between
# two consecutive estimates of the best VaR is less than epsilon
UVaR=Inf
delta2=Inf
while(delta2 >epsilon ){
UVaRtemp=UVaR
for(j in 1:d) {
UX[,j]=sort(UX[,j],decreasing=TRUE)[rank(rowSums(UX[,c(1:d)[-j]]))] }
UVaR= max(rowSums(UX))
delta2=abs(UVaR-UVaRtemp)
}
# OUTPUT
#BEST VAR is the interval [LVaR,UVaR].
print(LVaR)
print(UVaR)