VaR_ES_bounds: Worst and Best Value-at-Risk and Best Expected Shortfall for Given Marginals

Description

Compute the worst and best Value-at-Risk (VaR) and the best expected shortfall (ES) for given marginal distributions.

Usage

## Homogeneous case
crude_VaR_bounds(level, qF, d = NULL, ...)
VaR_bounds_hom(level, d, method = c("Wang", "Wang.Par", "dual"),
               interval = NULL, tol = NULL, ...)
dual_bound(s, d, pF, tol = .Machine$double.eps^0.25, ...)

## Inomogeneous case
## Workhorses
rearrange(X, tol = 0, tol.type = c("relative", "absolute"),
          n.lookback = ncol(X), max.ra = Inf,
          method = c("worst.VaR", "best.VaR", "best.ES"),
	  sample = TRUE, is.sorted = FALSE, trace = FALSE, ...)
block_rearrange(X, tol = 0, tol.type = c("absolute", "relative"),
                n.lookback = ncol(X), max.ra = Inf,
                method = c("worst.VaR", "best.VaR", "best.ES"),
                sample = TRUE, trace = FALSE, ...)
## User interfaces
## Rearrangement Algorithm
RA(level, qF, N, abstol = 0, n.lookback = length(qF), max.ra = Inf,
   method = c("worst.VaR", "best.VaR", "best.ES"), sample = TRUE)
## Adaptive Rearrangement Algorithm
ARA(level, qF, N.exp = seq(8, 19, by = 1), reltol = c(0, 0.01),
    n.lookback = length(qF), max.ra = 10*length(qF),
    method = c("worst.VaR", "best.VaR", "best.ES"),
    sample = TRUE)
## Adaptive Block Rearrangement Algorithm
ABRA(level, qF, N.exp = seq(8, 19, by = 1), absreltol = c(0, 0.01),
     n.lookback = NULL, max.ra = Inf,
     method = c("worst.VaR", "best.VaR", "best.ES"),
     sample = TRUE)

Arguments

level

confidence level $\alpha$ for VaR and ES (e.g., 0.99).

dimension (number of risk factors; $\ge 2$). For crude_VaR_bounds(), d only needs to be given in the homogeneous case in which qF is a function.

d-list containing the marginal quantile functions. For crude_VaR_bounds(), qF can also be function (homogeneous case, dimension d).

method

character string. For

VaR_bounds_hom():: method = "Wang" and method = "Wang.Par" apply the approach of McNeil et al. (2015, Proposition 8.32) for computing best (i.e., smallest) and worst (i.e., largest) VaR. The latter method assumes Pareto margins and thus does not require numerical integration. method = "dual" applies the dual bound approach as in Embrechts et al. (2013, Proposition 4) for computing worst VaR (no value for the best VaR can be obtained with this approach and thus NA is returned for the best VaR).
rearrange(), block_rearrange(), RA(), ARA():: method indicates whether bounds for the worst/best VaR or the best ES should be computed. These bounds are termed $\underline{s}_N$ and $\overline{s}_N$ in the literature (and below) and are theoretically not guaranteed bounds of worst/best VaR or best ES; however, they are treated as such in practice and are typically in line with results from VaR_bounds_hom() in the homogeneous case, for example.

interval

initial interval (a numeric(2)) for computing worst VaR. If not provided, these are the defaults chosen:

method = "Wang":: initial interval is $[0,(1-\alpha)/d]$.
method = "Wang.Par":: initial interval is $[c_l,c_u]$, where $c_l$ and $c_u$ are chosen as in Hofert et al. (2015).
method = "dual":: in this case, no good defaults are known. Note that the lower endpoint of the initial interval has to be sufficiently large in order for the the inner root-finding algorithm to find a root; see Details.

tol

VaR_bounds_hom():: tolerance for uniroot() for computing worst VaR. This defaults (for tol = NULL) to $2.2204*10^{-16}$ for method = "Wang" or method = "Wang.Par" (where a smaller tolerance is crucial) and to uniroot()'s default .Machine$double.eps^0.25 otherwise. Note that for method = "dual", tol is used for both the outer and the inner root-finding procedure.
rearrange(), block_rearrange():: (absolute or relative) tolerance to determine (the individual) convergence. This should normally be a number greater than or equal to 0, but rearrange() also allows for tol = NULL which means that columns are rearranged until each column is oppositely ordered to the sum of all other columns.

tol.type

character string indicating the type of convergence tolerance function to be used ("relative" for relative tolerance and "absolute" for absolute tolerance).

n.lookback

number of rearrangements to look back for deciding about numerical convergence. Use this option with care.

dual bound evaluation point.

marginal loss distribution function.

(N, d)-matrix of quantiles (to be rearranged). If is.sorted it is assumed that the columns of X are sorted in increasing order.

max.ra

maximal number of (considered) column rearrangements of the underlying matrix of quantiles (can be set to Inf).

number of discretization points.

N.exp

exponents of the number of discretization points (a vector) over which the algorithm iterates to find the smallest number of discretization points for which the desired accuracy (specified by abstol and reltol) is attained; for each number of discretization points, at most max.ra-many column rearrangements are of the underlying matrix of quantiles are considered.

abstol

absolute convergence tolerance $\epsilon$ to determine the individual convergence, i.e., the change in the computed minimal row sums (for method = "worst.VaR") or maximal row sums (for method = "best.VaR") or expected shortfalls (for method = "best.ES") for the lower bound $\underline{s}_N$ and the upper bound $\overline{s}_N$. abstol is typically $\ge0$; it can also be NULL, see tol above.

reltol

vector of length two containing the individual (first component; used to determine convergence of the minimal row sums (for method = "worst.VaR") or maximal row sums (for method = "best.VaR") or expected shortfalls (for method = "best.ES") for $\underline{s}_N$ and $\overline{s}_N$) and the joint (second component; relative tolerance between the computed $\underline{s}_N$ and $\overline{s}_N$ with respect to $\overline{s}_N$) relative convergence tolerances. reltol can also be of length one in which case it denotes the joint relative tolerance; the individual relative tolerance is taken as NULL (see tol above) in this case.

absreltol

vector of length two containing the individual (first component; used to determine convergence of the minimal row sums (for method = "worst.VaR") or maximal row sums (for method = "best.VaR") or expected shortfalls (for method = "best.ES") for $\underline{s}_N$ and $\overline{s}_N$) absolute and the joint (second component; relative tolerance between the computed $\underline{s}_N$ and $\overline{s}_N$ with respect to $\overline{s}_N$) relative convergence tolerances. absreltol can also be of length one in which case it denotes the joint relative tolerance; the individual absolute tolerance is taken as 0 in this case.

sample

logical indicating whether each column of the two underlying matrices of quantiles (see Step 3 of the Rearrangement Algorithm in Embrechts et al. (2013)) are randomly permuted before the rearrangements begin. This typically has quite a positive effect on run time (as most of the time is spent (oppositely) ordering columns (for rearrange()) or blocks (for block_rearrange())).

is.sorted

logical indicating whether the columns of X are sorted in increasing order.

trace

logical indicating whether the underlying matrix is printed after each rearrangement step. See vignette("VaR_bounds", package = "qrmtools") for how to interpret the output.

…

crude_VaR_bounds():: ellipsis argument passed to (all provided) quantile functions.
VaR_bounds_hom():: case method = "Wang" requires the quantile function qF() to be provided and additional arguments passed via the ellipsis argument are passed on to the underlying integrate(). For method = "Wang.Par" the ellipsis argument must contain the parameter shape (the shape parameter $\theta>0$ of the Pareto distribution). For method = "dual", the ellipsis argument must contain the distribution function pF() and the initial interval interval for the outer root finding procedure (not for d = 2); additional arguments are passed on to the underlying integrate() for computing the dual bound $D(s)$.
dual_bound():: ellipsis argument is passed to the underlying integrate().
rearrange():: additional arguments passed to the underlying optimization function. Currently, this is only used if method = "best.ES" in which case the required confidence level $\alpha$ must be provided as argument level.

Value

crude_VaR_bounds() returns crude lower and upper bounds for VaR at confidence level $\alpha$ for any $d$-dimensional model with marginal quantile functions specified by qF.

VaR_bounds_hom() returns the best and worst VaR at confidence level $\alpha$ for $d$ risks with equal distribution function specified by the ellipsis ....

dual_bound() returns the value of the dual bound $D(s)$ as given in Embrechts, Puccetti, R<U+00FC>schendorf (2013, Eq. (12)).

rearrange() and block_rearrange() return a list containing

bound:: computed $\underline{s}_N$ or $\overline{s}_N$.
tol:: reached tolerance (i.e., the (absolute or relative) change of the minimal row sum (for method = "worst.VaR") or maximal row sum (for method = "best.VaR") or expected shortfall (for method = "best.ES") after the last rearrangement).
converged:: logical indicating whether the desired (absolute or relative) tolerance tol has been reached.
opt.row.sums:: vector containing the computed optima (minima for method = "worst.VaR"; maxima for method = "best.VaR"; expected shortfalls for method = "best.ES") for the row sums after each (considered) rearrangement.
X.rearranged:: (N, d)-matrix containing the rearranged X.

RA() returns a list containing

bounds:: bivariate vector containing the computed $\underline{s}_N$ and $\overline{s}_N$ (the so-called rearrangement range) which are typically treated as bounds for worst/best VaR or best ES; see also above.
rel.ra.gap:: reached relative tolerance (also known as relative rearrangement gap) between $\underline{s}_N$ and $\overline{s}_N$ computed with respect to $\overline{s}_N$.
ind.abs.tol:: bivariate vector containing the reached individual absolute tolerances (i.e., the absolute change of the minimal row sums (for method = "worst.VaR") or maximal row sums (for method = "best.VaR") or expected shortfalls (for mehtod = "best.ES") for computing $\underline{s}_N$ and $\overline{s}_N$; see also tol returned by rearrange() above).
converged:: bivariate logical vector indicating convergence of the computed $\underline{s}_N$ and $\overline{s}_N$ (i.e., whether the desired tolerances were reached).
num.ra:: bivariate vector containing the number of column rearrangments of the underlying matrices of quantiles for $\underline{s}_N$ and $\overline{s}_N$.
opt.row.sums:: list of length two containing the computed optima (minima for method = "worst.VaR"; maxima for method = "best.VaR"; expected shortfalls for method = "best.ES") for the row sums after each (considered) column rearrangement for the computed $\underline{s}_N$ and $\overline{s}_N$; see also rearrange().
X:: initially constructed (N, d)-matrices of quantiles for computing $\underline{s}_N$ and $\overline{s}_N$.
X.rearranged:: rearranged matrices X (for $\underline{s}_N$ and $\overline{s}_N$).

ARA() and ABRA() return a list containing

bounds:: see RA().
rel.ra.gap:: see RA().
tol:: trivariate vector containing the reached individual (relative for ARA(); absolute for ABRA()) tolerances and the reached joint relative tolerance (computed with respect to $\overline{s}_N$).
converged:: trivariate logical vector indicating individual convergence of the computed $\underline{s}_N$ (first entry) and $\overline{s}_N$ (second entry) and indicating joint convergence of the two bounds according to the attained joint relative tolerance (third entry).
N.used:: actual N used for computing the (final) $\underline{s}_N$ and $\overline{s}_N$.
num.ra:: see RA(); computed for N.used.
opt.row.sums:: see RA(); computed for N.used.
X:: see RA(); computed for N.used.
X.rearranged:: see RA(); computed for N.used.

Details

For d = 2, VaR_bounds_hom() uses the method of Embrechts et al. (2013, Proposition 2). For method = "Wang" and method = "Wang.Par" the method presented in McNeil et al. (2015, Prop. 8.32) is implemented; this goes back to Embrechts et al. (2014, Prop. 3.1; note that the published version of this paper contains typos for both bounds). This requires one uniroot() and, for the generic method = "Wang", one integrate(). The critical part for the generic method = "Wang" is the lower endpoint of the initial interval for uniroot(). If the (marginal) distribution function has finite first moment, this can be taken as 0. However, if it has infinite first moment, the lower endpoint has to be positive (but must lie below the unknown root). Note that the upper endpoint $(1-\alpha)/d$ also happens to be a root and thus one needs a proper initional interval containing the root and being stricticly contained in $(0,(1-\alpha)/d$. In the case of Pareto margins, Hofert et al. (2015) have derived such an initial (which is used by method = "Wang.Par"). Also note that the chosen smaller default tolerances for uniroot() in case of method = "Wang" and method = "Wang.Par" are crucial for obtaining reliable VaR values; see Hofert et al. (2015).

For method = "dual" for computing worst VaR, the method presented of Embrechts et al. (2013, Proposition 4) is implemented. This requires two (nested) uniroot(), and an integrate(). For the inner root-finding procedure to find a root, the lower endpoint of the provided initial interval has to be “sufficiently large”.

Note that these approaches for computing the VaR bounds in the homogeneous case are numerically non-trivial; see the source code and vignette("VaR_bounds", package = "qrmtools") for more details. As a rule of thumb, use method = "Wang" if you have to (i.e., if the margins are not Pareto) and method = "Wang.Par" if you can (i.e., if the margins are Pareto). It is not recommended to use (the numerically even more challenging) method = "dual".

Concerning the inhomogeneous case, rearrange() is an auxiliary function (workhorse). It is called by RA() and ARA(). After a column rearrangement of X, the tolerance between the minimal row sum (for the worst VaR) or maximal row sum (for the best VaR) or expected shortfall (obtained from the row sums; for the best ES) after this rearrangement and the one of $n.lookback$ rearrangement steps before is computed and convergence determined. For performance reasons, no input checking is done for rearrange() and it can change in future versions to (futher) improve run time. Overall it should only be used by experts.

block_rearrange(), the workhorse underlying ABRA(), is similar to rearrange() in that it checks whether convergence has occurred after every rearrangement by comparing the change to the row sum variance from n.lookback rearrangement steps back. block_rearrange() differs from rearrange in the following ways. First, instead of single columns, whole (randomly chosen) blocks (two at a time) are chosen and oppositely ordered. Since some of the ideas for improving the speed of rearrange() do not carry over to block_rearrange(), the latter should in general not be as fast as the former. Second, instead of using minimal or maximal row sums or expected shortfall to determine numerical convergence, block_rearrange() uses the variance of the vector of row sums to determine numerical convergence. By default, it targets a variance of 0 (which is also why the default tol.type is "absolute").

For the Rearrangement Algorithm RA(), convergence of $\underline{s}_N$ and $\overline{s}_N$ is determined if the minimal row sum (for the worst VaR) or maximal row sum (for the best VaR) or expected shortfall (obtained from the row sums; for the best ES) satisfies the specified abstol (so $\le\epsilon$) after at most max.ra-many column rearrangements. This is different from Embrechts et al. (2013) who use $<\epsilon$ and only check for convergence after an iteration through all columns of the underlying matrix of quantiles has been completed.

For the Adaptive Rearrangement Algorithm ARA() and the Adaptive Block Rearrangement Algorithm ABRA(), convergence of $\underline{s}_N$ and $\overline{s}_N$ is determined if, after at most max.ra-many column rearrangements, the (the individual relative tolerance) reltol[1] is satisfied and the relative (joint) tolerance between both bounds is at most reltol[2].

Note that RA(), ARA() and ABRA() need to evalute the 0-quantile (for the lower bound for the best VaR) and the 1-quantile (for the upper bound for the worst VaR). As the algorithms, due to performance reasons, can only handle finite values, the 0-quantile and the 1-quantile need to be adjusted if infinite. Instead of the 0-quantile, the $\alpha/(2N)$-quantile is computed and instead of the 1-quantile the $\alpha+(1-\alpha)(1-1/(2N))$-quantile is computed for such margins (if the 0-quantile or the 1-quantile is finite, no adjustment is made).

rearrange(), block_rearrange(), RA(), ARA() and ABRA() compute $\underline{s}_N$ and $\overline{s}_N$ which are, from a practical point of view, treated as bounds for the worst (i.e., largest) or the best (i.e., smallest) VaR or the best (i.e., smallest ES), but which are not known to be such bounds from a theoretical point of view; see also above. Calling them “bounds” for worst/best VaR or best ES is thus theoretically not correct (unless proven) but “practical”. The literature thus speaks of $(\underline{s}_N, \overline{s}_N)$ as the rearrangement gap.

More details not provided here can be found in the references listed below.

References

Embrechts, P., Puccetti, G., R<U+00FC>schendorf, L., Wang, R. and Beleraj, A. (2014). An Academic Response to Basel 3.5. Risks 2(1), 25--48.

Embrechts, P., Puccetti, G. and R<U+00FC>schendorf, L. (2013). Model uncertainty and VaR aggregation. Journal of Banking \& Finance 37, 2750--2764.

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Hofert, M., Memartoluie, A., Saunders, D. and Wirjanto, T. (2017). Improved Algorithms for Computing Worst Value-at-Risk. Statistics \& Risk Modeling or, for an earlier version, http://arxiv.org/abs/1505.02281.

Bernard, C., R<U+00FC>schendorf, L. and Vanduffel, S. (2013). Value-at-Risk bounds with variance constraints. See http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2342068.

Bernard, C. and McLeish, D. (2014). Algorithms for Finding Copulas Minimizing Convex Functions of Sums. See http://arxiv.org/abs/1502.02130v3.

Examples

Run this code

# NOT RUN {
### 1 Reproducing selected examples of McNeil et al. (2015; Table 8.1) #########

## Setup
alpha <- 0.95
d <- 8
theta <- 3
qF <- rep(list(function(p) qPar(p, shape = theta)), d)

## Worst VaR
N <- 5e4
set.seed(271)
system.time(RA.worst.VaR <- RA(alpha, qF = qF, N = N, method = "worst.VaR"))
RA.worst.VaR$bounds
stopifnot(RA.worst.VaR$converged,
          all.equal(RA.worst.VaR$bounds[["low"]],
                    RA.worst.VaR$bounds[["up"]], tol = 1e-4))

## Best VaR
N <- 5e4
set.seed(271)
system.time(RA.best.VaR <- RA(alpha, qF = qF, N = N, method = "best.VaR"))
RA.best.VaR$bounds
stopifnot(RA.best.VaR$converged,
          all.equal(RA.best.VaR$bounds[["low"]],
                    RA.best.VaR$bounds[["up"]], tol = 1e-4))

## Best ES
N <- 5e4 # actually, we need a (much larger) N here (but that's time consuming)
set.seed(271)
system.time(RA.best.ES <- RA(alpha, qF = qF, N = N, method = "best.ES"))
RA.best.ES$bounds
stopifnot(RA.best.ES$converged,
          all.equal(RA.best.ES$bounds[["low"]],
                    RA.best.ES$bounds[["up"]], tol = 5e-1))


### 2 More Pareto examples (d = 2, d = 8; hom./inhom. case; explicit/RA/ARA) ###

alpha <- 0.99 # VaR confidence level
th <- 2 # Pareto parameter theta
qF <- function(p, theta = th) qPar(p, shape = theta) # Pareto quantile function
pF <- function(q, theta = th) pPar(q, shape = theta) # Pareto distribution function


### 2.1 The case d = 2 #########################################################

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.VaR")
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.VaR")
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))


### 2.2 The case d = 8 #########################################################

d <- 8 # dimension

## Compute VaR bounds with various methods
I <- crude_VaR_bounds(alpha, qF = qF, d = 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", shape = 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.VaR")
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$opt.row.sums$low)), n = 1),
                            tail(abs(diff(ARAbest$opt.row.sums$up)), n = 1)),
              method = "best.VaR")
RAworst <- RA(alpha, qF = rep(list(qF), d), N = ARAworst$N.used,
              abstol = mean(tail(abs(diff(ARAworst$opt.row.sums$low)), n = 1),
                            tail(abs(diff(ARAworst$opt.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$opt.row.sums, length)))
ylim <- range(RAworst$opt.row.sums)
plot(RAworst$opt.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$opt.row.sums[[1]]), RAworst$opt.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()


### 3 "Reproducing" examples from Embrechts et al. (2013) ######################

### 3.1 "Reproducing" Table 1 (but seed and eps are unknown) ###################

## Left-hand side of Table 1
N <- 50
d <- 3
qPar <- rep(list(qF), d)
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),]


### 3.2 "Reproducing" Table 3 for alpha = 0.99 #################################

## Note: The seed for obtaining the exact results as in Table 3 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, shape = xi[j], scale = beta[j])})
set.seed(271)
res.best <- RA(0.99, qF = qF.lst, N = N, abstol = eps, method = "best.VaR")
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 row


### 4 Further checks ###########################################################

## Calling the workhorses directly
set.seed(271)
ra <- rearrange(X)
bra <- block_rearrange(X)
stopifnot(ra$converged, bra$converged,
          all.equal(ra$bound, bra$bound, tolerance = 6e-3))

## Checking ABRA against ARA
set.seed(271)
ara  <- ARA (alpha, qF = qPar)
abra <- ABRA(alpha, qF = qPar)
stopifnot(ara$converged, abra$converged,
          all.equal(ara$bound[["low"]], abra$bound[["low"]], tolerance = 2e-3),
          all.equal(ara$bound[["up"]],  abra$bound[["up"]],  tolerance = 6e-3))
# }

Run the code above in your browser using DataLab