1 Homogeneous case
We start by considering the following setup in the homogeneous case, that is, when all marginal distributions are equal.
qF2 <- function(p, th = 2) qPar(p, shape = th) # Par(2) quantile function
pF2 <- function(q, th = 2) pPar(q, shape = th) # Par(2) distribution function
dim <- 8 # variable dimension (we use 8 or 100 here)1.1 Checks for method = “dual”
Investigate the helper function \(h(s,t)\) (the function for the inner
root-finding to compute \(D(s)\); see
dual_bound()).
d <- 8 # dimension
s <- c(1, 5, 10, 100, 500, 1000)
t <- sapply(seq_along(s), function(i) {
res <- exp(seq(log(1e-3), log(s[i]/d), length.out = 257))
res[length(res)] <- s[i]/d # to avoid numerical issues (t > s/d)
res
})
f <- sapply(seq_along(s), function(i)
sapply(t[,i], function(t.)
qrmtools:::dual_bound_2(s[i], t = t., d = d, pF = pF2) -
qrmtools:::dual_bound_2_deriv_term(s[i], t = t., d = d, pF = pF2)))
palette <- colorRampPalette(c("maroon3", "darkorange2", "royalblue3"), space = "Lab")
cols <- palette(6)
if(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_hom_dual_h_Par=2_d=",d,".pdf")),
width = 6, height = 6)
plot(t[,1], f[,1], type = "l", log = "x", xlim = range(t), ylim = range(f), col = cols[1],
xlab = "t", ylab = expression("h(s,t) for d = 8 and F being Par(2)"))
lines(t[,2], f[,2], col = cols[2])
lines(t[,3], f[,3], col = cols[3])
lines(t[,4], f[,4], col = cols[4])
lines(t[,5], f[,5], col = cols[5])
lines(t[,6], f[,6], col = cols[6])
abline(h = 0, lty = 2)
legend("topright", lty = rep(1,6), col = cols,
bty = "n", legend = as.expression(lapply(1:6,
function(i) substitute(s==s., list(s. = s[i])))))
if(doPDF) dev.off()
As we know, \(h(s,s/d) = 0\). We also
see that \(s\) has to be sufficiently
large in order to find a root \(h(s,t) =
0\) for \(t < s/d\).
Now let’s plot the dual bound \(D(s)\) for various \(\theta\) (this checks the outer root-finding).
theta <- c(0.5, 1, 2, 4) # theta values
s <- seq(48, 2000, length.out = 257) # s values
D <- sapply(theta, function(th)
sapply(s, function(s.)
dual_bound(s., d = d, pF = function(q) pPar(q, shape = th)))) # (s, theta) matrix
if(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_hom_dual_D_s_Par=",
paste0(theta, collapse="_"),"_d=",d,".pdf")),
width = 6, height = 6)
plot(s, D[,1], type = "l", ylim = range(D), col = "maroon3",
ylab = substitute("Dual bound D(s) for d ="~d.~"Par("*theta*") margins", list(d. = d)))
lines(s, D[,2], col = "darkorange2")
lines(s, D[,3], col = "royalblue3")
lines(s, D[,4], col = "black")
legend("topright", lty = rep(1,4),
col = c("maroon3", "darkorange2", "royalblue3", "black"),
bty = "n", legend = as.expression(lapply(1:4,
function(j) substitute(theta==j, list(j = theta[j])))))
if(doPDF) dev.off()
1.2 Checks for method = “Wang”/“Wang.Par”
1.2.1 Check of auxiliary functions with numerical integration (for \(\theta = 2\))
Check Wang_h_aux().
d <- 8 # dimension
alpha <- 0.99 # confidence level
c <- seq(0, (1-alpha)/d, length.out = 129) # domain of h
h.aux <- qrmtools:::Wang_h_aux(c, level = alpha, d = d, qF = qF2)
par(mar = c(5, 4+1, 4, 2) + 0.1) # increase space (for y axis label)
plot(c, h.aux, type = "l", xlab = "c (in initial interval)",
ylab = expression(frac(d-1,d)~{F^{-1}}(a[c])+frac(1,d)~{F^{-1}}(b[c])))
Check the objective function \(h(c)\) (that is, Wang_h() with
numerical integration) on its domain. Note that \(h\) is used to determine the root in the
open interval \((0,(1-\alpha)/d)\) for
computing worst value-at-risk.
h <- sapply(c, function(c.) qrmtools:::Wang_h(c., level = alpha, d = d, qF = qF2))
if(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_",alpha,"_hom_Wang_h_Par=2_d=",d,"_num.pdf")),
width = 6, height = 6)
plot(c, h, type = "l", xlab = "c (in initial interval)",
ylab = substitute("h(c) for"~~alpha~"= 0.99 and d ="~d.~"Par(2) margins",
list(d. = d)))
abline(h = 0, lty = 2)
if(doPDF) dev.off()
Check the objective function \(h\) at the endpoints of its domain.
sapply(c(0, (1-alpha)/d), function(c.)
qrmtools:::Wang_h(c., level = alpha, d = d, qF = qF2)) # -Inf, 0## [1] -Inf 0\(-\infty\) is not a problem for
root finding (for \(\theta > 1\);
for \(\theta <= 1\) it is
NaN, see below, and thus a problem!), but the \(0\) at the right endpoint is a problem.
method <- "Wang.Par" # this also holds for (the numerical) method = "Wang"
th <- 0.99
qrmtools:::Wang_h(0, level = alpha, d = d, method = method, shape = th) # NaN => uniroot() fails## [1] NaN## Note: Wang_h() is actually already NaN for c <= 1e-17
qrmtools:::Wang_h_aux(0, level = alpha, d = d, method = method, shape = th) # Inf## [1] InfA proper initial interval \([c_l,u_l]\) with \(0<c_l<=c_u<(1-\alpha)/d\) (containing the root) is thus required. We have derived this in Hofert et al. (2017, ``Improved Algorithms for Computing Worst Value-at-Risk’’)).
1.2.2 Check of \(h(c)\) without numerical integration (for a range of \(\theta\))
Check the objective function \(h(c)\) (Wang_h() without
numerical integration)
d <- dim # dimension
alpha <- 0.99 # confidence level
theta <- c(0.1, 0.5, 1, 5, 10, 50) # theta values
palette <- colorRampPalette(c("darkorange2", "maroon3", "royalblue3", "black"), space = "Lab")
cols <- palette(length(theta))
c <- seq(0, (1-alpha)/d, length.out = 2^13+1)
## => They all go further down to 0 if length.out is increased.
## Smaller theta thus corresponds to a larger derivative in the root
## Root-finding thus requires higher precision for smaller theta
h <- matrix(, nrow = length(c), ncol = length(theta))
for(j in 1:length(theta))
h[,j] <- sapply(c, function(c.)
qrmtools:::Wang_h(c., level = alpha, d = d, method = "Wang.Par", shape = theta[j]))
z <- h
z[z <= 0] <- NA # > 0 => makes log-scale possibleif(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_",alpha,"_hom_Wang_h_Par_d=",d,".pdf")),
width = 6, height = 6)
plot(NA, xlim = range(c), ylim = range(z, na.rm = TRUE), log = "y", yaxt = "n", xlab = "c",
ylab = substitute("h(c) for"~~alpha~"= 0.99, d ="~d.~"Par("*theta*") margins",
list(d. = d)))
eaxis(2)
for(k in 1:length(theta))
lines(c, z[,k], col = cols[k])
legend("topleft", bty = "n", lty = rep(1,length(theta)), col = cols,
legend = as.expression(lapply(1:length(theta),
function(k) substitute(theta==k, list(k = theta[k])))))
if(doPDF) dev.off()
1.3 Compute best/worst \(\mathrm{VaR}_\alpha\) (via “Wang.Par”)
After dealing with various numerical issues, we can now look at some example calculations of best/worst value-at-risk. We first plot value-at-risk as a function of \(\alpha\).
d <- dim # dimension
alpha <- 1-2^seq(-0.001, -10, length.out = 128) # confidence levels; concentrated near 1
theta <- c(0.1, 0.5, 1, 5, 10, 50) # theta values
VaR <- simplify2array(sapply(alpha, function(a)
sapply(theta, function(th) VaR_bounds_hom(a, d = d, method = "Wang.Par",
shape = th)), simplify = FALSE))
## => (best/worst VaR, theta, alpha)-matrixif(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_hom_Wang_Par_d=",d,".pdf")),
width = 7, height = 7)
par(mar = c(5, 4+1, 4, 2) + 0.1) # increase space (for y axis label)
plot(NA, xlim = range(alpha), ylim = range(VaR), log = "xy", yaxt = "n",
xlab = expression(1-alpha),
ylab = as.expression(substitute(underline(VaR)[alpha](L^{"+"})~"(dashed) and"~bar(VaR)[alpha](L^{"+"})~
"(solid) for d ="~d.~"and Par("*theta*") margins",
list(d. = d))))
eaxis(2)
for(k in 1:length(theta)) {
lines(1-alpha, VaR[2,k,], col = cols[k]) # worst VaR
lines(1-alpha, VaR[1,k,], col = cols[k], lty = 2) # best VaR
}
legend("topright", bty = "n", lty = rep(1,length(theta)), col = cols,
legend = as.expression(lapply(1:length(theta),
function(k) substitute(theta==k, list(k = theta[k])))))
if(doPDF) dev.off()
Now consider best/worst value-at-risk as a function of \(d\).
d <- seq(2, 1002, by = 20) # dimensions
alpha <- 0.99 # confidence level
theta <- c(0.1, 0.5, 1, 5, 10, 50) # theta values
VaR <- simplify2array(sapply(d, function(d.)
sapply(theta, function(th) VaR_bounds_hom(alpha, d = d., method = "Wang.Par",
shape = th)), simplify = FALSE))
## => (best/worst VaR, theta, d)-matrixif(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_hom_Wang_Par_alpha=",alpha,"_in_d.pdf")),
width = 7, height = 7)
par(mar = c(5, 4+1, 4, 2) + 0.1) # increase space (for y axis label)
plot(NA, xlim = range(d), ylim = range(VaR), log = "xy", yaxt = "n",
xlab = expression(d),
ylab = as.expression(substitute(underline(VaR)[a](L^{"+"})~"(dashed) and"~bar(VaR)[a](L^{"+"})~
"(solid) for Par("*theta*") margins",
list(a = alpha))))
eaxis(2)
for(k in 1:length(theta)) {
lines(d, VaR[2,k,], col = cols[k]) # worst VaR
lines(d, VaR[1,k,], col = cols[k], lty = 2) # best VaR
}
legend("topleft", bty = "n", lty = rep(1,length(theta)), col = cols,
legend = as.expression(lapply(1:length(theta),
function(k) substitute(theta==k, list(k = theta[k])))))
if(doPDF) dev.off()
We can also consider best/worst value-at-risk as a function of \(\theta\). Note that, as before, best value-at-risk is the same for all \(d\) here (depicted in a green dashed line below).
d <- c(2, 10, 20, 100, 200, 1000) # dimensions
theta <- 10^seq(-1, log(50, base = 10), length.out = 50) # theta values
alpha <- 0.99 # confidence level
VaR <- simplify2array(sapply(d, function(d.)
sapply(theta, function(th) VaR_bounds_hom(alpha, d = d., method = "Wang.Par",
shape = th)), simplify = FALSE))
## => (best/worst VaR, theta, d)-matrixif(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_hom_Wang_Par_alpha=",alpha,"_in_theta.pdf")),
width = 6, height = 6)
par(mar = c(5, 4+1, 4, 2) + 0.1) # increase space (for y axis label)
plot(NA, xlim = range(theta), ylim = range(VaR), log = "xy", yaxt = "n",
xlab = expression(theta),
ylab = as.expression(substitute(underline(VaR)[a](L^{"+"})~"(dashed) and"~
bar(VaR)[a](L^{"+"})~"(solid) for Par("*theta*") margins",
list(a = alpha))))
eaxis(2)
bestVaR <- VaR[1,,1]
for(k in 1:length(d)) {
lines(theta, VaR[2,,k], col = cols[k]) # worst VaR
if(k >= 2) stopifnot(all.equal(bestVaR, VaR[1,,k]))
}
lines(theta, bestVaR, col = "darkgreen", lty = 2) # best VaR
legend("topright", bty = "n", lty = rep(1,length(d)), col = cols,
legend = as.expression(lapply(1:length(d),
function(k) substitute(d==k, list(k = d[k])))))
if(doPDF) dev.off()
1.4 Comparison between various methods for computing worst value-at-risk
Now consider a graphical comparison between the various methods. To stress the numerical challenges, let us include the following implementation of Wang’s approach (with various switches to highlight the individual numerical challenges) for the Pareto case.
## Initial interval for the root finding in case of worst VaR
init_interval <- function(alpha, d, shape, trafo = FALSE, adjusted = FALSE)
{
if(trafo) {
low <- if(shape == 1) {
d/2
} else {
(d-1)*(1+shape)/(d-1+shape)
}
up <- if(shape > 1) {
r <- (1+d/(shape-1))^shape
if(adjusted) 2*r else r
} else if(shape == 1) {
e <- exp(1)
(d+1)^(e/(e-1))
} else {
d*shape/(1-shape)+1
}
c(low, up)
} else {
low <- if(shape > 1) {
r <- (1-alpha)/((d/(shape-1)+1)^shape + d-1)
if(adjusted) r/2 else r
} else if(shape == 1) {
e <- exp(1)
(1-alpha)/((d+1)^(e/(e-1))+d-1)
} else {
r <- (1-shape)*(1-alpha)/d
if(adjusted) r/2 else r
}
up <- if(shape == 1) (1-alpha)/(3*d/2-1)
else (1-alpha)*(d-1+shape)/((d-1)*(2*shape+d))
c(low, up)
}
}
## Function to compute the best/worst value-at-risk in the homogeneous case with
## Par(theta) margins
VaR_hom_Par <- function(alpha, d, shape, method = c("worst", "best"),
trafo = FALSE, interval = NULL, adjusted = FALSE,
avoid.cancellation = FALSE, ...)
{
## Pareto quantile function
qF <- function(p) (1 - p)^(-1/shape) - 1
## Compute \bar{I}
Ibar <- function(a, b, alpha, d, shape)
{
if(shape == 1) log((1-a)/(1-b))/(b-a) - 1
else (shape/(1-shape))*((1-b)^(1-1/shape)-(1-a)^(1-1/shape))/(b-a) - 1
}
## Main
method <- match.arg(method)
switch(method,
"worst" = {
## Distinguish according to whether we optimize the auxiliary function
## on a transformed scale
h <- if(trafo) {
## Auxiliary function to find the root of on (1, Inf)
if(shape == 1) {
function(x) x^2 + x*(-d*log(x)+d-2)-(d-1)
} else {
function(x)
(d/(1-shape)-1)*x^(-1/shape + 1) - (d-1)*x^(-1/shape) + x - (d*shape/(1-shape) + 1)
}
} else {
## Auxiliary function to find the root of on (0, (1-alpha)/d)
function(c) {
a <- alpha+(d-1)*c
b <- 1-c
Ib <- if(c == (1-alpha)/d) { # Properly deal with limit c = (1-alpha)/d
((1-alpha)/d)^(-1/shape) - 1
} else {
Ibar(a = a, b = b, alpha = alpha, d = d, shape = shape)
}
Ib - (qF(a)*(d-1)/d + qF(b)/d)
}
}
## Do the optimization
if(is.null(interval)) interval <- init_interval(alpha, d, shape,
trafo = trafo, adjusted = adjusted)
c <- uniroot(h, interval = interval, ...)$root
if(trafo) # convert value back to the right scale (c-scale)
c <- (1-alpha)/(c+d-1)
if(avoid.cancellation) {
t1 <- (1-alpha)/c-(d-1)
d * ((c^(-1/shape)/d) * ((d-1)*t1^(-1/shape) + 1) - 1) # = qF(a)*(d-1) + qF(b)
} else {
a <- alpha+(d-1)*c
b <- 1-c
qF(a)*(d-1) + qF(b)
}
},
"best" = {
max((d-1)*0 + (1-alpha)^(-1/shape)-1, # Note: Typo in Wang, Peng, Yang (2013)
d*Ibar(a = 0, b = alpha, alpha = alpha, d = d, shape))
},
stop("Wrong 'method'"))
}For the comparison, consider the following setup.
alpha <- 0.99 # confidence level
d <- dim # dimension
n.th <- 32 # number of thetas
th <- seq(0.2, 5, length.out = n.th) # thetas
qFs <- lapply(th, function(th.) {th.; function(p) qPar(p, shape = th.)}) # n.th-vector of Pareto quantile functions
pFs <- lapply(th, function(th.) {th.; function(q) pPar(q, shape = th.)}) # n.th-vector of Pareto dfs
N <- 1e4 # number of discretization points for RA(); N = 1e5 does not improve the situationNow compute the values and plot them.
res <- matrix(, nrow = n.th, ncol = 7)
colnames(res) <- c("Wang", "straightforward", "transformed", "Wang.Par",
"dual", "RA.low", "RA.up")
pb <- txtProgressBar(max = n.th, style = if(isatty(stdout())) 3 else 1) # setup progress bar
for(i in seq_len(n.th)) {
## "Wang" (numerical integration with smaller uniroot() tolerance; still
## numerically critical -- we catch "the integral is probably divergent"-errors here)
Wang.num.res <- tryCatch(VaR_bounds_hom(alpha, d = d, qF = qFs[[i]])[2], error = function(e) e)
res[i,"Wang"] <- if(is(Wang.num.res, "simpleError")) NA else Wang.num.res
## Our straightforward implementation
res[i,"straightforward"] <- VaR_hom_Par(alpha, d = d, shape = th[i])
## Our straightforward implementation based on the transformed auxiliary function
res[i,"transformed"] <- VaR_hom_Par(alpha, d = d, shape = th[i], trafo = TRUE)
## "Wang.Par" (using a smaller uniroot() tolerance and adjusted initial interval)
res[i,"Wang.Par"] <- VaR_bounds_hom(alpha, d = d, method = "Wang.Par", shape = th[i])[2]
## "dual" (with uniroot()'s default tolerance)
res[i,"dual"] <- VaR_bounds_hom(alpha, d = d, method = "dual",
interval = crude_VaR_bounds(alpha, qF = qFs[[i]], d = d),
pF = pFs[[i]])[2]
## Rearrangement Algorithm
set.seed(271) # use the same random permutation for each theta
RA. <- RA(alpha, qF = rep(qFs[i], d), N = N)
res[i,"RA.low"] <- RA.$bounds[1]
res[i,"RA.up"] <- RA.$bounds[2]
## Progress
setTxtProgressBar(pb, i) # update progress bar
}
close(pb) # close progress barres. <- res/res[,"dual"] # standardize (by dual method)
ylim <- range(res., na.rm = TRUE)
if(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_",alpha,"_hom_comparison_d=",
d,"_N=",N,".pdf")), width = 7, height = 7)
par(mar = c(5, 4+1, 4, 2) + 0.1) # increase space (for y axis label)
plot(th, res.[,"Wang"], type = "l", ylim = ylim,
xlab = expression(theta), ylab = substitute("Standardized (by dual method)"~
bar(VaR)[0.99]~"for d ="~d.~"Par("*theta*") margins", list(d. = d)),
col = "gray60", lty = 2, lwd = 5.5) # Wang (with numerical integration)
lines(th, res.[,"straightforward"], col = "maroon3", lty = 1, lwd = 1) # still bad (although we have an initial interval)
lines(th, res.[,"transformed"], col = "black", lty = 1, lwd = 1) # okay
lines(th, res.[,"Wang.Par"], col = "royalblue3", lty = 2, lwd = 2.5) # Wang Pareto (wo num. integration)
lines(th, res.[,"RA.low"], col = "black", lty = 3, lwd = 1) # lower RA bound
lines(th, res.[,"RA.up"], col = "black", lty = 2, lwd = 1) # upper RA bound
legend("topright", bty = "n",
col = c("gray60", "maroon3", "black", "royalblue3", "black", "black"),
lty = c(2,1,1,2,3,2), lwd = c(5.5,1,1,2.5,1,1),
legend = c("Wang (num. int.)", "Wang Pareto (straightforward)",
"Wang Pareto (transformed)", "Wang Pareto (wo num. int.)",
"Lower RA bound", "Upper RA bound"))
if(doPDF) dev.off()
So what goes wrong with our straightforward implementation? Let’s
start with the obvious. As we can infer from the x-axis scale in the
plots of \(h\) above,
uniroot()’s default tolerance
tol = .Machine$double.eps^0.25(\(\approx 0.0001220703\) here) is too large
to determine the root accurately. In fact, if we choose a smaller
tol, we have no problem.
tol <- 2.2204e-16
wVaR.tol <- sapply(th, function(th.)
VaR_hom_Par(alpha = alpha, d = d, shape = th., tol = tol))
plot(th, wVaR.tol/res[,"dual"], type = "l", ylim = ylim,
xlab = expression(theta), ylab = "Wang's approach (straightforward) but with smaller tol")
So is this the solution to all the problems? No. Consider the setup as before, where \(d\) is running (we catch the appearing errors here).
alpha <- 0.99 # confidence level
d <- seq(2, 1002, by = 20) # dimensions
theta <- c(0.1, 0.5, 1, 5, 10, 50) # theta values
VaR <- simplify2array(sapply(d, function(d.)
sapply(theta, function(th) {
res <- tryCatch(VaR_hom_Par(alpha, d = d., shape = th, tol = tol),
error = function(e) e)
if(is(res, "simpleError")) { warning(conditionMessage(res)); NA }
else res
}), simplify = FALSE))By using warnings() we see that we cannot even determine
the root in all cases as the function values are numerically
not of opposite sign at the endpoints of the theoretically
correct initial interval containing the root. We can solve this
(non-elegantly) by simply adjusting the lower bound.
VaR <- simplify2array(sapply(d, function(d.)
sapply(theta, function(th) VaR_hom_Par(alpha, d = d., shape = th, tol = tol, adjusted = TRUE)),
simplify = FALSE))So we see that we need a smaller tolerance and to extend the initial interval in order to get reasonable results.
What about using the transformed auxiliary function directly? Is this a numerically more stable approach?
## [1] InfAs we can see, there is a problem as well. To understand it, consider
the following minimized version of VaR_hom_Par() (for
computing worst value-at-risk for \(\theta\neq
1\) only).
h <- function(x)
(d/(1-th)-1)*x^(-1/th + 1) - (d-1)*x^(-1/th) + x - (d*th/(1-th) + 1)
interval <- init_interval(alpha, d, th, trafo = TRUE)
x <- uniroot(h, interval = interval)$root
(c <- (1-alpha)/(x+d-1)) # convert back to c-scale## [1] 1.869497e-31## [1] InfAs we can see, \(c>0\) is so
small that \(b = 1-c\) is numerically
equal to 1 and thus the Pareto quantile evaluated as \(\infty\). This cancellation can be avoided
by simplifying the terms (note that
qPar(1-c, shape = th) = c^(-1/th)-1).
## [1] 162.5923We can apply this to the variable combinations as before. Note that due to the partly extreme values of \(d\) and \(\theta\), we also need to adjust the range here in order for a root to be found.
d <- seq(2, 1002, by = 20) # dimensions
theta <- c(0.1, 0.5, 1, 5, 10, 50) # theta values
VaR. <- simplify2array(sapply(d, function(d.)
sapply(theta, function(th) VaR_hom_Par(alpha, d = d., shape = th,
trafo = TRUE, adjusted = TRUE,
avoid.cancellation = TRUE)),
simplify = FALSE))We can now compare the values (based on the non-transformed/transformed auxiliary function) in a graph.
ylim <- range(VaR, VaR.)
if(doPDF)
pdf(file = (file <- paste0("fig_worst_VaR_",alpha,"_hom_comparison_num_problems.pdf")),
width = 7, height = 7)
par(mar = c(5, 4+1, 4, 2) + 0.1) # increase space (for y axis label)
plot(d, VaR[1,], type = "l", ylim = ylim, log = "y", yaxt = "n", xlab = "d",
ylab = as.expression(substitute(bar(VaR)[a](L^{"+"})~
"directly (solid) or with transformed h (dashed) for Par("*theta*") margins",
list(a = alpha))), col = cols[1])
eaxis(2)
for(k in 2:length(theta)) lines(d, VaR [k,], col = cols[k])
for(k in 1:length(theta)) lines(d, VaR.[k,], col = cols[k], lty = 2, lwd = 2)
legend("right", bty = "n",
col = cols, lty = rep(1, length(theta)),
legend = as.expression( lapply(1:length(theta), function(k)
substitute(theta==th, list(th = theta[k]))) ))
if(doPDF) dev.off()
We see that the results differ slightly if both \(d\) and \(\theta\) are large. It remains numerically challenging to appropriately compute worst value-at-risk in this case. However, both of the above approaches seem to work for dimensions even as large as \(10^5\).
If we order the realizations in each column to be in increasing order
(so the smallest (largest) entry of the first column lies next to the
smallest (largest) entry of the second column of 
We take away that column rearrangements don’t change the marginal
distributions but the dependence between the two columns and thus the
distributions of their sums. To see the latter (in terms of boxplots and
kernel density estimates), consider:
We see that if one of the three components is large, the other two have
to compensate for that (to keep the variance of the sum small, which is
what the (A)RA tries to optimize) and so the locally singular components
(strings to (1, 0, 0), (0, 1, 0), (0, 0, 1)) appear.