Implementing the Rearrangement Algorithm: An Example from Computational Risk Management

Abstract

After a brief overview of aspects of computational risk management, the implementation of the rearrangement algorithm in R is considered as an example from computational risk management practice. This algorithm is used to compute the largest quantile (worst value-at-risk) of the sum of the components of a random vector with specified marginal distributions. It is demonstrated how a basic implementation of the rearrangement algorithm can gradually be improved to provide a fast and reliable computational solution to the problem of computing worst value-at-risk. Besides a running example, an example based on real-life data is considered. Bootstrap confidence intervals for the worst value-at-risk as well as a basic worst value-at-risk allocation principle are introduced. The paper concludes with selected lessons learned from this experience.

1. Introduction

Computational risk management is a comparably new and exciting field of research at the intersection of statistics, computer science and data science. It is concerned with computational problems of quantitative risk management, such as algorithms, their implementation, computability, numerical robustness, parallel computing, pitfalls in simulations, run-time optimization, software development and maintenance. At the end of the (business) day, solutions in quantitative risk management are to be run on a computer of some sort and thus concern computational risk management.

There are various factors that can play a role when developing and implementing a solution to a problem from computational risk management in practice, for example:

(1)

Theoretical hurdles. Some solutions cannot be treated analytically anymore. Suppose $\mathit{X}=(X_1,\cdots ,X_d)$ is a d-dimensional random vector of risk factor changes, for example negative log-returns, with marginal distribution functions $F_1,\cdots ,F_d$ and we are interested in computing the probability $\mathbb{P}(X_j\in (F_j^{\leftarrow }\left(p_j\right),F_j^{\leftarrow }\left(q_j\right)],\mathrm{for}\mathrm{all}j)$ of the jth risk factor change to take on values between its $p_j$- and its $q_j$-quantile for all j. Even if we knew the joint distribution function H of $\mathit{X}$, computing such a probability analytically involves evaluating H at $2^d$ different values. For $d\approx 32$, this can already be time-consuming (not to speak of numerical inaccuracies appearing, which render the formula as good as useless) and for $d\approx 260$ the number of points at which to evaluate H is roughly equal to the estimated number of atoms in the universe.

(2)

Model risk. The risk of using the wrong model or a model in which assumptions are not fulfilled; every solution we implement is affected by model risk to some degree. In the aforementioned example, not knowing the exact H leads to model risk when computing the given probability by Monte Carlo integration.

(3)

The choice of software. The risk of using the wrong software; in the above example, a software not suitable for simulating from joint distribution functions H. Another example is to use a programming language too low-level for the problem at hand, requiring to implement standard tasks (sampling, optimization, etc.) manually and thus bearing the risk of obtaining unreliable results because of many possible pitfalls one might encounter when developing the tools needed to implement the full-blown solution. These aspects become more crucial nowadays since efficient solutions to different problems are often only available in one software each but companies need to combine these solutions and thus the corresponding software in order to solve these different problems at hand.

(4)

Syntax errors. These are compilation errors or execution errors by an interpreter because of a violation of the syntax of the programming language under consideration. Syntax errors are typically easy to detect as the program simply stops to run. Also, many programming languages provide tools for debugging to find the culprits.

(5)

Run-time errors. These are errors that appear while a program runs. Quite often, run-time errors are numerical errors, errors that appear, for example, because of the floating-point representation of numbers in double precision. Run-time errors can sometimes be challenging to find, for example, when they only happen occasionally in a large-scale simulation study and introduce a non-suspicious bias in the results. Run-time errors can sometimes also be hard to reproduce.

(6)

Semantic errors. These are errors where the code is syntactically correct and the program runs, but it does not compute what is intended (for example, due to a flipped logical statement). If results do not look suspicious or pass all tests conducted, sometimes measuring run time can be the only way to find such problems (for example, if a program finishes much earlier than expected).

(7)

Scaling errors. These are errors of an otherwise perfectly fine running code that fails when run at large scale. This could be due to a lack of numerical robustness, the sheer run time of the code, an exploding number of parameters involved, etc.

(8)

User errors. These are errors caused by users of the software solution; for example, when calling functions with wrong arguments because of a misinterpretation of their meaning. The associated risk of wrongly using software can, by definition, often be viewed as part of operational risk.

(9)

Warnings. Warnings are important signs of what to watch out for; for example, a warning could indicate that a numerical optimization routine has not converged after a predetermined number of steps and only the current best value is returned, which might be far off a local or the global optimum. Unfortunately, especially in large-scale simulation studies, users often suppress warnings instead of collecting and considering them.

(10)

Run time. The risk of using or propagating method A over B because the run time of A beats the one of B by a split second or second, not realizing that run time depends on factors such as the hardware used, the current workload, the  algorithm implemented, the programming language used, the implementation style, compiler flags, whether garbage collection was used, etc. There is not even a unique concept of time (system vs. user vs. elapsed time).

(11)

Development and maintenance. It is significantly more challenging to provide a robust, well developed, maintained and documented class of solutions as part of a bigger, coherent software package rather than a single script with hard-coded values to solve a special case of the same problem. Although stand-alone scripts get specific tasks done by 4 p.m., having a software package available is typically beneficial mid- to long-term. It can significantly reduce the risk of the aforementioned errors by reusing code well tested and applied by the users of the package or it can avoid the temptation of introducing errors long after the code was developed just because it suddenly looks suspicious although it is actually correct.

At the core of all solutions to problems from computational risk management lies an algorithm, a well-defined (unambiguous) finite set of instructions (steps) for solving a problem. An implementation of an algorithm in a programming language allows one to see how a problem is actually solved, unlike a formulation in pseudo-code, which is often vague and thus opens the door for misinterpretation. A classical example is if pseudo-code says “Minimize the function…” without mentioning how initial values or intervals can be obtained. Another example is “Choose a tolerance $\epsilon >0$…” but one is left in the dark concerning what suitable tolerances $\epsilon$ are for the problem at hand; they often depend on the unknown output one is interested to compute in the first place, see Step (1) of Algorithm 1 later. Oftentimes, these are the hardest parts to solve of the whole problem and it is important for research and scientific journals to accept corresponding contributions as important instead of brushing them off as “not innovative” or representing “no new contribution”.

A fundamental principle when developing an implementation is that of a minimal working example. A minimal working example is source code that is a working example in the sense that it allows someone else to reproduce a problem (sufficiency) and it is minimal in the sense that it is as small and as simple as possible, without non-relevant code, data or dependencies (necessity). Minimal working examples are often constructed from existing larger chunks of code by divide and conquer, that is by breaking down the problem into sub-problems and commenting out non-relevant parts until the code becomes simple enough to show the problem (which is then easier to grasp and solve or can be sent to an expert in the field to ask for help).

In this paper, we consider the problem of implementing the rearrangement algorithm (RA) of Embrechts et al. (2013) as an example from computational risk management. The RA allows one to compute, for example, the worst value-at-risk of the sum of d risks of which the marginal distributions are known but the dependence is unknown; it can also be used to compute the best value-at-risk or expected shortfall, but we focus on the worst value-at-risk. Section 2 contains the necessary details about the RA. Section 3 addresses how the major workhorse underlying this algorithm can be implemented in a straightforward way in R. This turns out to be inefficient and Section 4 presents ways to improve the implementation. Section 5 utilizes the implementations in our R package $\mathtt{qrmtools}$ for tracing and to motivate the chosen default tolerances, and Section 6 considers a real data example, presents a bootstrap confidence interval for worst value-at-risk and introduces a basic capital allocation principle based on worst value-at-risk. The lessons learned throughout the paper are summarized in Section 7.

2. The Rearrangement Algorithm in a Nutshell

Let $L_1\sim F_1,\cdots ,L_d\sim F_d$ be continuously distributed loss random variables over a predetermined period and let $L^{+}={\sum }_{j=1}^d{L}_j$ be the total loss over that time period. From a risk management perspective we are interested in computing value-at-risk (${VaR}_{\alpha }\left(L^{+}\right)$, ${VaR}_{\alpha }$ or VaR) at confidence level $\alpha \in (0,1)$, that is the $\alpha$-quantile $F_{L^{+}}^{\leftarrow }\left(\alpha \right)=inf\{x\in \mathbb{R}:F_{L^{+}}\left(x\right)\ge \alpha \}$ of the distribution function $F_{L^{+}}$ of $L^{+}$; in typical applications, $\alpha \in [0.99,1)$. If we had enough joint observations, so realizations of the random vector $\mathit{L}=(L_1,\cdots ,L_d)$, we could estimate ${VaR}_{\alpha }\left(L^{+}\right)$ empirically. A major problem is that one often only has realizations of each of $L_j$ individually, non-synchronized. This typically allows to pin down $F_1,\cdots ,F_d$ but not the joint distribution function H of $\mathit{L}$. By Sklar’s Theorem, such H can be written as

$$\begin{array}{c}\hfill H\left(\mathit{x}\right)=C(F_1\left(x_1\right),\cdots ,F_d\left(x_d\right)),\mathit{x}=(x_1,\cdots ,x_d)\in {\mathbb{R}}^d,\end{array}$$

for a unique copula C. In other words, although we often know or can estimate $F_1,\cdots ,F_d$, we typically neither know the copula C nor have enough joint realizations to be able to estimate it. However, the copula C determines the dependence between $L_1,\cdots ,L_d$ and thus the distribution of $L^{+}$ as the following example shows.

Example 1

(Motivation for rearrangements).

Consider $d=2$ and $L_j\sim F_j\left(x\right)=1-x^{-1/2}$, $x\ge 1$, $j=1,2$, that is both losses are Pareto Type I distributed with parameter $1/2$. The left-hand side of Figure 1 shows $n=1000$ realizations of $\mathit{L}=(L_1,L_2)=(F_1^{\leftarrow }\left(U_1\right),F_2^{\leftarrow }\left(U_2\right))$, once under independence, so for $(U_1,U_2)\sim {\mathrm{U}(0,1)}^2$, and once under comonotonicity, so for $U_1=U_2\sim \mathrm{U}(0,1)$. We see that the different dependence structures directly influence the shape of the realizations of $\mathit{L}$.

As the middle of Figure 1 shows, the dependence also affects ${VaR}_{\alpha }\left(L^{+}\right)$. This example is chosen to be rather extreme, (perhaps) in (stark) contrast to intuition, ${VaR}_{\alpha }\left(L^{+}\right)$ under independence is even larger than under comonotonicity, for all α; this can also be shown analytically in this case in a rather tedious calculation; see (Hofert et al. 2020, Exercise 2.28).

The right-hand side of Figure 1 shows the same plot with the marginal distributions now being gamma. In this case, comonotonicity leads to larger values of ${VaR}_{\alpha }\left(L^{+}\right)$ than independence, but only for sufficiently large α; the corresponding turning point seems to get larger the smaller the shape parameter of the second margin, for example.

Note that we did not utilize the stochastic representation $(L_1,L_2)=(F_1^{\leftarrow }\left(U\right),F_2^{\leftarrow }\left(U\right))$, $U\sim \mathrm{U}(0,1)$, for sampling $\mathit{L}=(L_1,L_2)$. Instead, we sorted the realizations of $L_2$ so that their ranks are aligned with those of $L_1$, in other words, so that the ith-smallest realization of $L_2$ lies next to the ith-smallest realizations of $L_1$. The rows of this sample thus consist of $(L_{\left(i\right)1},L_{\left(i\right)2})$, where $L_{\left(i\right)j}$ denotes the ith order statistic of the n realizations $L_{1j},\cdots ,L_{nj}$ of $L_j$. Such a rearrangement mimics comonotonicity between the realizations of $L_1$ and $L_2$. Note that this did not change the realizations of $L_1$ or $L_2$ individually, so it did not change the marginal realizations, only the joint ones.

As we learned from Example 1, by rearranging the marginal loss realizations, we can mimic different dependence structures between the losses and thus influence ${VaR}_{\alpha }\left(L^{+}\right)$. In practice, the worst ${VaR}_{\alpha }\left(L^{+}\right)$, denoted by ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$, is of interest, that is the largest ${VaR}_{\alpha }\left(L^{+}\right)$ for given margins $L_1\sim F_1,\cdots ,L_d\sim F_d$. The following remark motivates an objective when rearranging marginal loss realizations in order to increase ${VaR}_{\alpha }\left(L^{+}\right)$ and thus approximate ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$.

Remark 1

(Objective of rearrangements).

(1) 

By Example 1, the question which copula C maximizes ${VaR}_{\alpha }\left(L^{+}\right)$ is thus, at least approximately, the question of which reordering of realizations of each of $L_1,\cdots ,L_d$ leads to ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$. For any C, the probability of exceeding ${VaR}_{\alpha }\left(L^{+}\right)$ is (at most, but for simplicity let us assume exactly) $1-\alpha$. How that probability mass is distributed beyond ${VaR}_{\alpha }\left(L^{+}\right)$ depends on C. If we find a C such that $L^{+}$ has a rather small variance $var\left(L^{+}\right)$, more probability mass will be concentrated around a single point which can help to pull ${VaR}_{\alpha }\left(L^{+}\right)$ further into the right tail and thus increase ${VaR}_{\alpha }\left(L^{+}\right)$; Figure 2 provides a sketch in terms of a skewed $t_3$ density. If that single point exists and if $var\left(L^{+}\right)=0$, then it must be expected shortfall${ES}_{\alpha }\left(L^{+}\right)=\frac{1}{1-\alpha }{\int }_{\alpha }^1{VaR}_u\left(L^{+}\right)\mathrm{d}u$, which provides an upper bound to ${VaR}_{\alpha }\left(L^{+}\right)$.

(2) 

It becomes apparent from Figure 2 that the distribution of $L^{+}$ below its α-quantile is irrelevant for the location of ${VaR}_{\alpha }\left(L^{+}\right)$. It thus suffices to consider losses $L_1,\cdots ,L_d$ beyond their marginal α-quantiles $F_1^{\leftarrow }\left(\alpha \right),\cdots ,F_d^{\leftarrow }\left(\alpha \right)$, so the conditional losses $L_j|L_j>F_j^{\leftarrow }\left(\alpha \right)$, $j=1,\cdots ,d$, and their copula $C^{\alpha }$, called worst VaR copula. Note that since α is typically rather close to 1, the distribution of $L_j|L_j>F_j^{\leftarrow }\left(\alpha \right)$ (the tail of $L_j$) is typically modeled by a continuous distribution; this can be justified theoretically by the Pickands–Balkema–de-Haan Theorem, see (McNeil et al. 2015, Section 5.2.1).

(3) 

In general, if  $var\left(L^{+}\right)=0$ cannot be attained, the smallest point of the support of $L^{+}$ given that $L_j>F_j^{\leftarrow }\left(\alpha \right)$, $j=1,\cdots ,d$, is taken as an approximation to ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$.

The RA aims to compute ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$ by rearranging realizations of $L_j|L_j>F_j^{\leftarrow }\left(\alpha \right)$, $j=1,\cdots ,d$, such that their sum becomes more and more constant, in order to approximately obtain the smallest possible variance and thus the largest ${VaR}_{\alpha }\left(L^{+}\right)$. There are still fundamental open mathematical questions concerning the convergence of this algorithm, but this intuition is sufficient for us to understand the RA and to study its implementation. To this end, we now address the remaining underlying concepts we need.

For each margin $j=1,\cdots ,d$, the RA uses two discretizations of $F_j^{\leftarrow }$ beyond the probability $\alpha$ as realizations of $L_j|L_j>F_j^{\leftarrow }\left(\alpha \right)$ to be rearranged. In practice, one would typically utilize the available data from $L_j$ to estimate its distribution function $F_j$ and then use the implied quantile function $F_j^{\leftarrow }$ of the fitted $F_j$ to obtain such discretizations; see Section 6. Alternatively, $F_j^{\leftarrow }$ could be specified by expert opinion. The RA uses two discretizations to obtain upper and lower bounds to ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$, which, when sufficiently close, can be used to obtain an estimate of ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$, for example, by taking the midpoint. The two discretizations are stored in the matrices

$$\begin{array}{c}\hfill {\underline{X}}^{\alpha }={\left(F_j^{-}\left(\alpha +(1-\alpha )\frac{i-1}{N}\right)\right)}_{i=1,\cdots ,N}^{j=1,\cdots ,d}\mathrm{and}{\overline{X}}^{\alpha }={\left(F_j^{-}\left(\alpha +(1-\alpha )\frac{i}{N}\right)\right)}_{i=1,\cdots ,N}^{j=1,\cdots ,d}.\end{array}$$

(1)

These matrices are the central objects the RA works with. For simplicity of the argument, we write $X^{\alpha }=\left(X_{ij}^{\alpha }\right)$ if the argument applies to ${\underline{X}}^{\alpha }$ or ${\overline{X}}^{\alpha }$. The RA iterates over the columns of $X^{\alpha }$ and successively rearranges each column in a way such that the row sums ${\left({\sum }_{j=1}^d{X}_j^{\alpha }\right)}_{i=1,\cdots ,N}$ of the rearranged matrix have a smaller variance; compare with Remark 1.

In each rearrangement step, the RA rearranges the jth column ${\mathit{X}}_j^{\alpha }$ of $X^{\alpha }$ oppositely to the vector of row sums ${\mathit{X}}_{-j}^{\alpha }={\sum }_{k=1,k\ne j}^d{\mathit{X}}_k^{\alpha }$ of all other columns, where ${\mathit{X}}_k^{\alpha }=(X_{1k}^{\alpha },\cdots ,X_{Nk}^{\alpha })$. Two vectors $\mathit{a}=(a_1,\cdots ,a_N)$ and $\mathit{b}=(b_1,\cdots ,b_N)$ are called oppositely ordered if $(a_{i_1}-a_{i_2})(b_{i_1}-b_{i_2})\le 0$ for all $i_1,i_2=1,\cdots ,N$. After oppositely ordering, the (second-, third-, etc.) smallest component of $\mathit{a}$ lies next to the (second-, third-, etc.) largest of $\mathit{b}$, which helps decrease the variance of their componentwise sum. To illustrate this, consider the following example.

Example 2

(Motivation for oppositely ordering columns).

We start by writing an auxiliary function to compute a matrix $X^{\alpha }$ with the modification that for ${\overline{X}}^{\alpha }$, 1-quantiles $F_j^{-}\left(1\right)$ appearing in the last row are replaced by $F_j^{-}\left(\alpha +(1-\alpha )\frac{N-1/2}{N}\right)=F_j^{-}\left(1-\frac{1-\alpha }{2N}\right)$ to avoid possible infinities; we come back to this point later. It is useful to start writing functions at this point, as we can reuse them later and rely on them for computing partial results, inputs needed, etc. In R, the function $\mathtt{stopifnot}\left(\right)$ can be used to check inputs of functions. If any of the provided conditions fails, this would produce an error. For more elaborate error messages, one can work with $\mathtt{stop}\left(\right)$. Input checks are extremely important when functions are exported in a package (so visible to a user of the package) or even made available to users by simply sharing R scripts. Code is often run by users who are not experts in the underlying theory or mechanisms and good input checks can prevent them from wrongly using the shared software; see Section 1 Point (8). Also, documenting the function (here: Roxygen-style documentation) is important, providing easy to understand information about what the function computes, meaning of input parameters, return value and who wrote the function.

Consider $d=2$ and $L_j\sim F_j\left(x\right)=1-{(1+x)}^{-{\theta }_j}$, $x\ge 0$, $j=1,2$, that is each of the two losses is Pareto distributed; we choose ${\theta }_1=3/2$ and ${\theta }_2=2$ here. We start by building the list of marginal quantile functions and then compute the corresponding (here: upper) quantile matrix $X^{\alpha }(={\overline{X}}^{\alpha })$ for $\alpha =0.9$ and $N=10000$.

The columns in $X^{\alpha }$ are sorted in increasing order, so mimic comonotonicity in the joint tail, which leads to the following variance estimate of the row sums.

If we randomly shuffle the first column, so mimicking independence in the joint tail, we obtain the following variance estimate.

Now if we oppositely order the first column with respect to the sum of all others, here the second column, we can see that the variance of the resulting row sums will decrease; note that we mimic countermonotonicity in the joint tail in this case.

The two marginal distributions largely differ in their heavy-tailedness, which is why the variance of their sum, even when oppositely reordered is still rather large. For ${\theta }_1={\theta }_2=2$, one obtains $138.6592$, for two standard normal margins $0.0808$ and for two standard uniform distributions 0.

As a last underlying concept, the RA utilizes the minimal row sum operator

$$\begin{array}{c}\hfill s\left(X\right)=\underset{1\le i\le N}{min}\sum _{1\le j\le d}{X}_{ij}\mathrm{for}X\in {\mathbb{R}}^{N\times d},\end{array}$$

which is motivated in Remark 1 Part (3). We are now ready to provide the RA.

The RA as introduced in Embrechts et al. (2013) states that ${\underline{s}}_N\le {\overline{s}}_N$ and in practice ${\underline{s}}_N\approx {\overline{s}}_N\approx {\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$. Furthermore, the initial randomizations in Steps (2.2) and (3.2) are introduced to avoid convergence problems of ${\overline{s}}_N-{\underline{s}}_N\to 0$.

3. A First Implementation of the Basic Rearrangement Step

When implementing an algorithm such as the RA, it is typically a good idea to divide and conquer, which means breaking down the problem into sub-problems and addressing those first. Having implemented solutions to the sub-problems (write functions!), one can then use them as black boxes to implement the whole algorithm. In our case, we can already rely on the function $\mathtt{quantile}\_\mathtt{matrix}\left(\right)$ for Steps (2.1) and (3.1) of the RA and we have already seen how to oppositely order two columns in Example 2. We also see that apart from which matrix is rearranged, Steps (2) and (3) are equal. We thus focus on this step in what follows.

Example 3

(Basic rearrangements).

We start with a basic implementation for computing one of the bounds ${\underline{s}}_N,{\overline{s}}_N$ of the RA.

In comparison to Algorithm 1, our basic implementation already uses a relative rather than an absolute convergence tolerance ε (more intuitively named $\mathtt{tol}$ in our implementation), the former is more suitable here as we do not know ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$ and therefore cannot judge whether an absolute tolerance of, say, $0.001$ is reasonable. Also, we terminate if the attained relative tolerance is less than or equal to $\mathtt{tol}$ since including equality allows us to choose $\mathtt{tol}=\mathtt{0}$; see Steps (2.4) and (3.4) of Algorithm 1. Furthermore, we allow the tolerance to be $\mathtt{NULL}$ in which case the columns are rearranged until all of them are oppositely ordered with respect to the sum of all other. The choice $\mathtt{tol}=\mathtt{NULL}$ is good for testing, but typically results in much larger run times and rarely has an advantage over $\mathtt{tol}=\mathtt{0}$.

Algorithm 1: Rearrangement algorithm for computing bounds on ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$

(1)  Fix a confidence level $\alpha \in (0,1)$, marginal quantile functions $F_1^{-},\cdots ,F_d^{-}$, a number of discretization points $N\in \mathbb{N}$ and an absolute convergence tolerance $\epsilon \ge 0$. (2)  Compute the lower bound: (2.1.)  Define the matrix ${\underline{X}}^{\alpha }$ as in (1). (2.2.)     Permute randomly the elements in each column of ${\underline{X}}^{\alpha }$. (2.3.)     Set ${\underline{Y}}^{\alpha }={\underline{X}}^{\alpha }$. For $1\le j\le d$, rearrange the jth column of the matrix ${\underline{Y}}^{\alpha }$ so that it becomes oppositely ordered to the sum of all other columns. (2.4.)     While $s\left({\underline{Y}}^{\alpha }\right)-s\left({\underline{X}}^{\alpha }\right)\ge \epsilon$, set ${\underline{X}}^{\alpha }$ to ${\underline{Y}}^{\alpha }$ and repeat Step (2.3). (2.5.)     Set ${\underline{s}}_N=s\left({\underline{Y}}^{\alpha }\right)$. (3)  Compute the upper bound: (3.1.)  Define the matrix ${\overline{X}}^{\alpha }$ as in (1). (3.2.)     Permute randomly the elements in each column of ${\overline{X}}^{\alpha }$. (3.3.)     Set ${\overline{Y}}^{\alpha }={\overline{X}}^{\alpha }$. For $1\le j\le d$, rearrange the jth column of the matrix ${\overline{Y}}^{\alpha }$ so that it becomes oppositely ordered to the sum of all other columns. (3.4.)     While $s\left({\overline{Y}}^{\alpha }\right)-s\left({\overline{X}}^{\alpha }\right)\ge \epsilon$, set ${\overline{X}}^{\alpha }$ to ${\overline{Y}}^{\alpha }$ and repeat Step (3.3). (3.5.)     Set ${\overline{s}}_N=s\left({\overline{Y}}^{\alpha }\right)$. (4)  Return $({\underline{s}}_N,{\overline{s}}_N)$.

    The opposite ordering step in our basic implementation contains the argument $\mathtt{ties}.\mathtt{method}=\mathtt{"}\mathtt{first}\mathtt{"}$ of $\mathtt{rank}\left(\right)$, which specifies how ties are handled, namely those ties with smaller index get assigned the smaller rank. Although this was not a problem in Example 2, when N is large and $d>2$, ties can quickly arise numerically. It can be important to use a sorting algorithm that has deterministic behavior in this case, as otherwise, the RA might not terminate or not as quickly as it could; see the vignette $\mathtt{VaR}\_\mathtt{bounds}$ of the Rpackage $\mathtt{qrmtools}$ for more details.

Finally, we can call $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ on both matrices in Equation (1) to obtain the lower and upper bounds ${\underline{s}}_N$, ${\overline{s}}_N$ and thus implement the RA. We would thus reuse the same code for Steps (2) and (3) of the RA, which means less code to check, improve or maintain.

Example 4

(Using the basic implementation).

To call $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ in a running example, we use the following parameters and build the input matrix ${\overline{X}}^{\alpha }$.

We can now call $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ and also measure the elapsed time in seconds of this call.

4. Improvements

It is always good to have a first working version of an algorithm, for example, to compare against in case one tries to improve the code and it becomes harder to read or check. However, as we saw in Example 4, our basic implementation of Example 3 is rather slow even if we chose a rather small N in this running example. Our goal is thus to present ideas how to improve $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ as a building block of the RA.

4.1. Profiling

As a first step, we profile the last call to see where most of the run time is spent.

Profiling writes out the call stack every so-many split seconds, so checks where the execution currently is. This allows one to measure the time spent in each function call. Note that some functions do not create stack frame records and thus do not show up. Nevertheless, profiling the code is often helpful. It is typically easiest to consider the measurement in percentage (see second and fourth column in the above output) and often in total, so the run time spent in that particular function or the functions it calls (see the fourth column). As this column reveals, two large contributors to run time are the computation of row sums and the sorting.

4.2. Avoiding Summations

First consider the computation of row sums. In each rearrangement step, $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ computes the row sums of all but the current column. As we explained, this rearrangement step is mathematically intuitive but it is computationally expensive. An important observation is that the row sums of all but the current column is simply the row sums of all columns minus the current column; in short, ${\mathit{X}}_{-j}^{\alpha }=\left({\sum }_{k=1}^N{\mathit{X}}_k^{\alpha }\right)-{\mathit{X}}_j^{\alpha }$. Therefore, if we, after randomly permuting each column, compute the total row sum once, we can subtract the jth column to obtain the row sums of all but the jth column. After having rearranged the jth column, we then add the rearranged jth column to the row sums of all other columns in order to obtain the updated total row sums.

Example 5

(Improved rearrangements).

The following improved version of $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ incorporates said idea; to save space, we omit the function header.

Next, we compare $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$ with $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$. We check if we obtain the same result and also measure run time as a comparison.

Run time is substantially improved by about 85%. This should be especially helpful for large d, so let us investigate the percentage improvement; here only done for rather small d to show the effect. Figure 3 shows the relative speed-up in percentage of $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$ over $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ as a function of d for the running example introduced in Example 4.

We can infer that the percentage improvement of $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$ in comparison to $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ converges to 100% for large d. This is not a surprise since determining the row sums of all but a fixed column requires $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ to compute $N(d-1)$-many sums whereas $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$ requires only N-many. This is an improvement by $O\left(d\right)$, which especially starts to matter for large d.

A quick note on run-time measurement is in order here. In a more serious setup, one would take repeated measurements of run time (at least for smaller d), determine the average run time for each d and perhaps empirical confidence intervals for the true run time based on the repeated measurements. Also, we use the same seed for each method, which guarantees the same shuffling of columns and thus allows for a fair comparison. Not doing so can result in quite different run times and would thus indeed require repeated measurements to get more reliable results.

Finally, the trick of subtracting a single column from the total row sums in order to obtain the row sums of all other columns comes at a price. Although the RA as in Algorithm 1 works well if some of the 1-quantiles appearing in the last row of ${\overline{X}}^{\alpha }$ are infinity, this does not apply to $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$ anymore since if the last entry in the jth column is infinity, the last entry in the vector of total row sums is infinity and thus the last entry in the vector of row sums of all but the jth column would not be defined ($\mathtt{NaN}$ in R). This is why $\mathtt{quantile}\_\mathtt{matrix}\left(\right)$ avoids computing 1-quantiles; see Example 2.

4.3. Avoiding Accessing Matrix Columns and Improved Sorting

There are two further improvements we can consider.

The first concerns the opposite ordering of columns. It so far involved a call to $\mathtt{sort}\left(\right)$ and one to $\mathtt{rank}\left(\right)$. It turns out that these two calls can be replaced by a nested $\mathtt{order}\left(\right)$ call, which is slightly faster than $\mathtt{rank}\left(\right)$ alone. We start by demonstrating that we can replicate $\mathtt{rank}(,\mathtt{ties}.\mathtt{method}=\mathtt{"}\mathtt{first}\mathtt{"})$ by a nested $\mathtt{order}\left(\right)$ call based on standard uniform data.

Providing $\mathtt{rank}\left(\right)$ with a method to deal with ties is important here as the default assigns average ranks on ties and thus can produce non-integer numbers. One would not guess to see ties in such a small set of standard uniform random numbers, but that is not true; see Hofert (2020) for details.

To see how much faster the nested $\mathtt{order}\left(\right)$ call is than $\mathtt{rank}\left(\right)$, we consider a small simulation study. For each of the two methods, we measure elapsed time in seconds when applied 20 times to samples of random numbers of different sizes.

Figure 4 provides a summary in terms of a box plot. We see that the nested $\mathtt{order}\left(\right)$ call is faster than $\mathtt{rank}\left(\right)$; percentages by how much the latter is worse are displayed in the labels.

As a second improvement, recall that the RA naturally works on the columns of matrices. Matrices or arrays are internally stored as long vectors with attributes indicating when to “wrap around”. Accessing a matrix column thus requires to determine the beginning and the end of that column in the flattened vector. We can thus speed-up $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$ by working with lists of columns instead of matrices.

The following example incorporates these two further improvements and some more.

Example 6

(Advanced rearrangements).

The following implementation saves column-access time by working with lists rather than matrices. We also use the slightly faster $.\mathtt{rowSums}\left(\right)$ (as we know the dimensions of the input matrix) instead of $\mathtt{rowSums}\left(\right)$ for computing the initial row sums once in the beginning, and we incorporate the idea of a nested $\mathtt{order}\left(\right)$ call instead of $\mathtt{rank}\left(\right)$. Furthermore, before shuffling the columns, we also store their original sorted versions, which allows us to omit the call of $\mathtt{sort}\left(\right)$ in $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$.

We now check for correctness of $\mathtt{advanced}\_\mathtt{rearrange}\left(\right)$ and measure its run time in our running example. Concerning correctness, for the more involved $\mathtt{advanced}\_\mathtt{rearrange}\left(\right)$ it is especially useful to have simpler versions available to check against; we use the result of $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$ here as we already checked it, but could have equally well used the one of $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$.

Although the gain in run time is not as dramatic as before, we still see an improvement of about 62% in comparison to $\mathtt{improved}\_\mathtt{rearrange}\left(\right)$.

4.4. Comparison with the Actual Implementation

The computational improvements so far have already lead to a percentage improvement of 94% in comparison to $\mathtt{basic}\_\mathtt{rearrange}\left(\right)$ in our running example. The following example compares our fastest function $\mathtt{advanced}\_\mathtt{rearrange}\left(\right)$ so far with our implementation $\mathtt{rearrange}\left(\right)$ in the R package $\mathtt{qrmtools}$ whose development was motivated by a risk management team at a Swiss bank.

Example 7

(Comparison withrearrange().)

To save even more run time, the function $\mathtt{rearrange}\left(\right)$ from the Rpackage $\mathtt{qrmtools}$ splits the input matrix into its columns at C level. It also uses an early stopping criterion in the sense that after all columns have been rearranged once, stopping is checked after every column rearrangement. Because of the latter, the result is not quite comparable to our previous versions.

We did not measure run time here since, similar to Heisenberg’s uncertainty principle, the difference in run time becomes more and more random; this is in line with the point we made about run-time measurement in the introduction. In other words, the two implementations become closer in run time. This can be demonstrated by repeated measurements.

As we see from the box plot in Figure 5, on average, $\mathtt{rearrange}\left(\right)$ is slightly faster than $\mathtt{advanced}\_\mathtt{rearrange}\left(\right)$. Although $\mathtt{rearrange}\left(\right)$ could even be made faster, it also computes and returns more information about the rearrangements than $\mathtt{advanced}\_\mathtt{rearrange}\left(\right)$; for example, the computed minimal row sums after each column rearrangement and the row of the final rearranged matrix corresponding to the minimal row sum in our case here. Given that this is a function applied by users who are not necessarily familiar with the rearrangement algorithm, it is more important to provide such information. Nevertheless, $\mathtt{rearrange}\left(\right)$ is still faster than $\mathtt{advanced}\_\mathtt{rearrange}\left(\right)$ on average.

The function $\mathtt{rearrange}\left(\right)$ is the workhorse underlying the function $\mathtt{RA}\left(\right)$ that implements the rearrangement algorithm in the R package $\mathtt{qrmtools}$. The adaptive rearrangement algorithm (ARA) introduced in Hofert et al. (2017) improves the RA in that it introduces the joint tolerance, that is the tolerance between ${\underline{s}}_N$ and ${\overline{s}}_N$; the tolerance used so far is called individual tolerance. The ARA applies column rearrangements until the individual tolerance is met for each of the two bounds and until the joint tolerance is met or a maximal number of column rearrangements has been reached. If the latter is the case, N is increased and the procedure repeated. The ARA implementation $\mathtt{ARA}\left(\right)$ in $\mathtt{qrmtools}$ also returns useful information about how the algorithm proceeded. Overall, $\mathtt{ARA}\left(\right)$ neither requires to choose a tolerance nor the number of discretization points, which makes this algorithm straightforward to apply in practice where the choice of such tuning parameters is often unclear. Finding good tuning parameters (adaptively or not) is often a significant problem in newly suggested procedures and not rarely a question open for future research by itself.

5. The Functions $\mathtt{rearrange}\left(\right)$ and $\mathtt{ARA}\left(\right)$

In this section we utilize the implementations $\mathtt{rearrange}\left(\right)$ and $\mathtt{ARA}\left(\right)$ of the rearrangement step and the adaptive rearrangement algorithm.

Example 8

(Tracingrearrange().)

The implementation $\mathtt{rearrange}\left(\right)$ has a tracing feature. It can produce a lot of output but we here consider the rather minimal example of rearranging the matrix

$$\begin{array}{c}\hfill X=\left(\begin{array}{ccc}1& 1& 1\\ 2& 2& 2\\ 3& 3& 3\end{array}\right).\end{array}$$

The following code enables tracing ($\mathtt{trace}=\mathtt{TRUE}$) to see how the column rearrangements proceed; to this end we rearrange until all columns are oppositely ordered to the sum of all other columns ($\mathtt{tol}=\mathtt{NULL}$), disable random permutation of the columns ($\mathtt{sample}=\mathtt{FALSE}$) and provide the information that the columns are already sorted here ($\mathtt{is}.\mathtt{sorted}=\mathtt{TRUE}$).

A “|” or “=” symbol over the just worked on column indicates whether this column was changed (“|”) or not (“=”) in this rearrangement step. The last two printed columns provide the row sums of all other columns as well as the new updated total row sums after the rearrangement. As we can see from this example, the rearranged matrix has minimal row sum 5 and is given by

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}3& 1& 1\\ 2& 2& 2\\ 1& 3& 3\end{array}\right),\mathit{whereas}\mathit{the}\mathit{matrix}\left(\begin{array}{ccc}1& 2& 3\\ 2& 3& 1\\ 3& 1& 2\end{array}\right)\end{array}$$

would have led to the larger minimal row sum 6; this rather constructed example for when the greedy column rearrangements of the RA can fail to provide the maximal minimal row sum was provided by Haus (2015). As more interesting example to trace is to rearrange the matrix

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}1& 1& 1\\ 2& 3& 2\\ 3& 5& 4\\ 4& 7& 8\end{array}\right);\end{array}$$

see the vignette $\mathtt{VaR}\_\mathtt{bounds}$ of the Rpackage $\mathtt{qrmtools}$.

Example 9

(CallingARA())

We now call $\mathtt{ARA}\left(\right)$ in the case of our running example in this paper.

As we can see, the implementation returns a list with the computed bounds ${\underline{s}}_N,{\overline{s}}_N$, the relative rearrangement gap $({\overline{s}}_N-{\underline{s}}_N)/{\overline{s}}_N$, the two individual and one joint tolerance reached on stopping, a corresponding vector of logicals indicating whether the required tolerances have been met, the number N of discretization points used in the final step, the number of considered column rearrangements for each bound, vectors of computed optimal (for worst VaR this means “minimal”) row sums after each column rearrangement considered, a list with the two input matrices ${\underline{X}}^{\alpha },{\overline{X}}^{\alpha }$, a list with the corresponding final rearranged matrices ${\underline{Y}}^{\alpha },{\overline{Y}}^{\alpha }$ and the rows corresponding to their optimal row sums. An estimate of ${\overline{\mathrm{VaR}}}_{0.99}\left(L^{+}\right)$ in our running example can then be computed as follows.

Figure 6 shows the minimal row sum as a function of the number of rearranged columns for each of the two bounds ${\underline{s}}_N$ and ${\overline{s}}_N$ on ${\overline{\mathrm{VaR}}}_{0.99}$ in our running example. We can see that already after all columns were rearranged once, the actual value of the two bounds does not change much anymore in this case, yet the algorithm still proceeds until the required tolerances are met.

6. Applications

We now consider a real data example and introduce bootstrap confidence intervals for ${\overline{\mathrm{VaR}}}_{0.99}$ as well as a basic worst VaR allocation principle.

Example 10

(A real-life example).

We consider negative log-returns of Google, Apple and Microsoft stocks from 2006-01-03 to 2010-12-31.

We treat the negative log-returns as (roughly) stationary here and consider mean excess plots in Figure 7 to determine suitable thresholds for applying the peaks-over-threshold method.

We then fit generalized Pareto distributions (GPDs) to the excess losses over these thresholds.

We consider the corresponding quantile functions as approximate quantile functions to our empirical losses; see (McNeil et al. 2015, Section 5.2.3).

Based on these quantile functions, we can then compute ${\overline{\mathrm{VaR}}}_{0.99}\left(L^{+}\right)$ with $\mathtt{ARA}\left(\right)$.

Given the additional information obtained from $\mathtt{ARA}\left(\right)$, we can also visualize a (pseudo-)sample from the worst VaR copula $C^{\alpha }$, see Figure 8. Such a sample is obtained by computing the pseudo-observations, that is componentwise ranks scaled to $(0,1)$, of the rearranged matrix ${\overline{Y}}^{\alpha }$ corresponding to the upper bound ${\overline{s}}_N$; we could have also considered ${\underline{Y}}^{\alpha }$ here.

Also, we can investigate how much the ${\overline{\mathrm{VaR}}}_{0.99}\left(L^{+}\right)$ bounds ${\underline{s}}_N$ and ${\overline{s}}_N$ change with changing individual relative tolerance. Figure 9 motivates the use of 0 as individual relative tolerance; this is often quite a bit faster than $\mathtt{NULL}$ and provides similar ${\overline{\mathrm{VaR}}}_{0.99}\left(L^{+}\right)$ bounds.

Example 11

(Bootstrap confidence interval for worst VaR).

As in Example 10 and mentioned in Section 2, the marginal quantile functions used as input for the rearrangement algorithm may need to be estimated from data. The corresponding estimation error translates to a confidence interval for ${\overline{\mathrm{VaR}}}_{0.99}\left(L^{+}\right)$, which we can obtain from a bootstrap. To this end, we start by building B bootstrap samples from the negative log-returns; we choose a rather small $B=20$ here to demonstrate the idea.

Next, we build excess losses over 90% thresholds (a typical broad-brush choice in practice), fit GPDs to the d component samples for each of the B bootstrap samples and construct the corresponding quantile functions.

Now we can call $\mathtt{ARA}\left(\right)$ with each of the B sets of marginal quantile functions.

And then extract the computed estimates of ${\overline{\mathrm{VaR}}}_{0.99}\left(L^{+}\right)$.

A bootstrap 95%-confidence interval for ${\overline{\mathrm{VaR}}}_{0.99}\left(L^{+}\right)$ is thus given as follows.

Note that our run-time improvements of Section 4 are especially useful here as a bootstrap would otherwise be relatively time-consuming.

Example 12

(Basic worst VaR allocation).

Capital allocation concerns the allocation of a total capital K to d business lines so that $K={AC}_1+\cdots +{AC}_d$ (the full allocation property), where ${AC}_j$ denotes the capital allocated to business line j. As the RA uses the minimal row sum as ${\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$ estimate, the rows in the final rearranged matrices ${\underline{Y}}^{\alpha },{\overline{Y}}^{\alpha }$ provide an allocation of $K={\overline{\mathrm{VaR}}}_{\alpha }\left(L^{+}\right)$ one can call worst VaR allocation. There could be multiple rows in ${\underline{Y}}^{\alpha },{\overline{Y}}^{\alpha }$ leading to the minimal row sum in which case the componentwise average of these rows is considered. The resulting two rows are returned by ($\mathtt{rearrange}\left(\right)$ and) $\mathtt{ARA}\left(\right)$ in $\mathtt{qrmtools}$ version $\mathtt{0.0.13}$. By componentwise averaging, the latter two rows we obtain an approximate worst VaR allocation. Given that we already have a bootstrap sample, we can average the B computed allocations componentwise to obtain the bootstrap mean of the worst VaR allocation.

We can also compare the corresponding estimate of the total allocation ${AC}_1+\cdots +{AC}_d$ with the average of the B optimal row sums computed by the ARA to check the full allocation property; to this end, we first compute the B averages of the optimal row sums over the lower and upper bounds for each replication.

Bootstrap 95%-confidence intervals for each of ${AC}_1,\cdots ,{AC}_d$ can be obtained as follows.

For other applications of the RA, see, for example, Embrechts and Jakobsons (2016), Bernard et al. (2017), Bernard et al. (2018) and, most notably, Ramsey and Goodwin (2019) where the RA is applied in the context of crop insurance.

7. Selected Lessons Learned

After using R to motivate column rearrangements and opposite ordering for finding the worst VaR for given marginal distributions in Section 2, we considered a basic implementation of such rearrangement in Section 3. We profiled the code in Section 4 and improved the implementation step-by-step by avoiding unnecessary summations, avoiding accessing matrix columns and improving the sorting step. In Section 5 we then used the implementations $\mathtt{rearrange}\left(\right)$ and $\mathtt{ARA}\left(\right)$ in our R package $\mathtt{qrmtools}$ to trace the rearrangement steps and to motivate the use of the default tolerances used by $\mathtt{ARA}\left(\right)$. A real data example was considered throughout Section 6 where we computed the worst VaR with $\mathtt{ARA}\left(\right)$, visualized a sample from the corresponding worst VaR copula and assessed the sensitivity of the computed worst VaR with respect to the individual convergence tolerance. To incorporate the uncertainty due to the estimation error of the marginal distributions, we used a bootstrap idea to obtain a confidence interval for the true worst VaR. We also introduced worst VaR allocation as a capital allocation principle and computed bootstrap confidence intervals for the allocated capitals.

Our experiments reveal lessons to learn that frequently appear in problems from the realm of computational risk management:

■

We can and should use software to learn about new problems and communicate their solutions. Implementing an algorithm often helps to fully comprehend a problem, experiment with it and get new ideas for practical solutions. It also allows one to communicate a solution since an implementation is unambiguous even if a theoretical presentation or pseudo-code of an algorithm leaves questions unanswered (such as how to choose tuning parameters); the latter is often the case when algorithms are provided in academic research papers.

■

Both for understanding a problem and for communicating its solution, it is paramount to use visualizations. Creating meaningful graphs is unfortunately still an exception rather than the rule in academic publications where tables are frequently used. Tables mainly allow one to see single numbers at a time, whereas graphs allow one to see the “bigger picture”, so how results are connected and behave as functions of the investigated inputs.

■

Given an algorithm, start with a basic, straightforward implementation and make sure it is valid, ideally by testing against known solutions in special cases; we omitted that part here, but note that there are semi-analytical solutions available for computing worst VaR in the case of homogeneous margins.

■

Learn from the basic implementation, experiment with it and improve it, for example, by improving numerical stability, run time, maintenance, etc. If you run into problems along the way, use minimal working examples to solve them; they are typically constructed by divide and conquer. As implementations get more involved, compare them against the basic implementation or previously checked improved versions.

■

When improving code, be aware of the issues mentioned in the introduction and exploit the fact that mathematically unique solutions often allow for different computational solutions. Also, computational problems (ties, parts of domains close to limits, etc.) can arise when there are none mathematically. These are often the hardest problems to solve, which is why it is important to be aware of computational risk management problems and their solutions. Implementing a computationally tractable solution is often one “dimension” more complicated than finding a mathematical solution—also literally, since a random variable (random vector) is often replaced by a vector (matrix) of realizations in a computational solution.

■

Results in academic research papers are often based on an implementation of a newly presented algorithm that works only under specific conditions with parameters found by experimentation (often hard-coded) under such conditions in the examples considered. A publicly available implementation in a software package needs to go well beyond this point; for example, as users typically do not have expert knowledge of problems the implementation intends to solve and often come from different backgrounds altogether in the hope to find solutions to challenging problems in their field of application. On the one hand, the amount of work spent on development and maintenance of publicly available software solutions is largely underestimated. On the other hand, one can sometimes benefit from getting to know different areas of application or valuable input for further theoretical or computational investigation through user feedback. Also, ones own implementation often helps one to explore potential further improvements or new solutions altogether.

■

Optimizing run time to be able to apply solutions in practical situations can be important but is not all there is. Providing a readable, comprehensible and numerically stable solution is equally important, for example. Code optimized solely for run time typically becomes less transparent and harder to maintain, is thus harder to adapt to future conceptual improvements and more prone to semantic errors. Also, a good solution typically not only provides a final answer but useful by-products and intermediate results computed along the way, just like a good solution in a mathematical exam requires intermediate steps to be presented in order to obtain full marks (which also simplifies to find the culprit if the final answer is erroneous).