require(qrmtools)
## Pareto setup
alpha <- 0.99 # VaR confidence level
th <- 2 # Pareto parameter theta
qF <- function(p, theta=th) qPar(p, theta=theta) # Pareto quantile function
pF <- function(q, theta=th) pPar(q, theta=theta) # Pareto distribution function
## d=2: Compute best/worst VaR explicitly (hom. case) and compare with (A)RA ###
d <- 2 # dimension
## Explicit
VaRbounds <- VaR_bounds_hom(alpha, d=d, qF=qF) # (best VaR, worst VaR)
## Adaptive Rearrangement Algorithm (ARA)
set.seed(271) # set seed (for reproducibility)
ARAbest <- ARA(alpha, qF=rep(list(qF), d), method="best")
ARAworst <- ARA(alpha, qF=rep(list(qF), d))
## Rearrangement Algorithm (RA) with N as in ARA()
RAbest <- RA(alpha, qF=rep(list(qF), d), N=ARAbest$N.used, method="best")
RAworst <- RA(alpha, qF=rep(list(qF), d), N=ARAworst$N.used)
## Compare
stopifnot(all.equal(c(ARAbest$bounds[1], ARAbest$bounds[2],
RAbest$bounds[1], RAbest$bounds[2]),
rep(VaRbounds[1], 4), tolerance=0.004, check.names=FALSE))
stopifnot(all.equal(c(ARAworst$bounds[1], ARAworst$bounds[2],
RAworst$bounds[1], RAworst$bounds[2]),
rep(VaRbounds[2], 4), tolerance=0.003, check.names=FALSE))
## d=8: Compute best/worst VaR (hom. case) and compare with (A)RA ##############
d <- 8 # dimension
## Compute VaR bounds with various methods
I <- crude_VaR_bounds(alpha, qF=rep(list(qF), d)) # crude bound
VaR.W <- VaR_bounds_hom(alpha, d=d, method="Wang", qF=qF)
VaR.W.Par <- VaR_bounds_hom(alpha, d=d, method="Wang.Par", theta=th)
VaR.dual <- VaR_bounds_hom(alpha, d=d, method="dual", interval=I, pF=pF)
## Adaptive Rearrangement Algorithm (ARA) (with different relative tolerances)
set.seed(271) # set seed (for reproducibility)
ARAbest <- ARA(alpha, qF=rep(list(qF), d), reltol=c(0.001, 0.01), method="best")
ARAworst <- ARA(alpha, qF=rep(list(qF), d), reltol=c(0.001, 0.01))
## Rearrangement Algorithm (RA) with N as in ARA and abstol (roughly) chosen as in ARA
RAbest <- RA(alpha, qF=rep(list(qF), d), N=ARAbest$N.used,
abstol=mean(tail(abs(diff(ARAbest$m.row.sums$low)), n=1),
tail(abs(diff(ARAbest$m.row.sums$up)), n=1)),
method="best")
RAworst <- RA(alpha, qF=rep(list(qF), d), N=ARAworst$N.used,
abstol=mean(tail(abs(diff(ARAworst$m.row.sums$low)), n=1),
tail(abs(diff(ARAworst$m.row.sums$up)), n=1)))
## Compare
stopifnot(all.equal(c(VaR.W[1], ARAbest$bounds, RAbest$bounds),
rep(VaR.W.Par[1],5), tolerance=0.004, check.names=FALSE))
stopifnot(all.equal(c(VaR.W[2], VaR.dual[2], ARAworst$bounds, RAworst$bounds),
rep(VaR.W.Par[2],6), tolerance=0.003, check.names=FALSE))
## Using (some of) the additional results computed by (A)RA() ##################
xlim <- c(1, max(sapply(RAworst$m.row.sums, length)))
ylim <- range(RAworst$m.row.sums)
plot(RAworst$m.row.sums[[2]], type="l", xlim=xlim, ylim=ylim,
xlab="Number or rearranged columns",
ylab=paste0("Minimal row sum per rearranged column"),
main=substitute("Worst VaR minimal row sums ("*alpha==a.*","~d==d.*"and Par("*
th.*"))", list(a.=alpha, d.=d, th.=th)))
lines(1:length(RAworst$m.row.sums[[1]]), RAworst$m.row.sums[[1]], col="royalblue3")
legend("bottomright", bty="n", lty=rep(1,2),
col=c("black", "royalblue3"), legend=c("upper bound", "lower bound"))
## => One should use ARA() instead of RA()
## "Reproducing" examples from Embrechts et al. (2013) #########################
## "Reproducing" Table 1 (but seed and eps are unknown)
## Left-hand side of Table 1
N <- 50
qPar <- rep(list(qF), 3)
p <- alpha + (1-alpha)*(0:(N-1))/N # for 'worst' (= largest) VaR
X <- sapply(qPar, function(qF) qF(p))
cbind(X, rowSums(X))
## Right-hand side of Table 1
set.seed(271)
res <- RA(alpha, qF=qPar, N=N)
row.sum <- rowSums(res$X.rearranged$low)
cbind(res$X.rearranged$low, row.sum)[order(row.sum),]
## "Reproducing" Table 3 for alpha=0.99 (but seed is unknown)
N <- 2e4 # we use a smaller N here to save run time
eps <- 0.1 # absolute tolerance
xi <- c(1.19, 1.17, 1.01, 1.39, 1.23, 1.22, 0.85, 0.98)
beta <- c(774, 254, 233, 412, 107, 243, 314, 124)
qF.lst <- lapply(1:8, function(j){ function(p) qGPD(p, xi=xi[j], beta=beta[j])})
set.seed(271)
res.best <- RA(0.99, qF=qF.lst, N=N, abstol=eps, method="best")
print(format(res.best$bounds, scientific=TRUE), quote=FALSE) # close to first value of 1st row
res.worst <- RA(0.99, qF=qF.lst, N=N, abstol=eps)
print(format(res.worst$bounds, scientific=TRUE), quote=FALSE) # close to last value of 1st rowRun the code above in your browser using DataLab