2.1 Notational preliminaries

Let the random variable X represent the change in assets less liabilities over a one-year horizon^{Footnote 4}. Negative values of X will represent losses and positive values will represent gains. The 1 in 200 loss is given by the 0.5th percentile of the distribution. The capital requirement is the amount required to be held by the firm to cover this potential loss.

The internal model associates a given risk factor $r\in\mathbb{R}^m$ with a loss $x(r)\in\mathbb{R}$ . We call x(r) the exact loss calculation since it contains no proxy error. The risk factor is considered as a random variable, denoted R, distributed according to a known multivariate distribution, so that the loss X satisfies the composition relationship:

(1) \begin{equation}X \sim x(R).\end{equation}

Here the notation $A\sim B$ means A and B have equal distribution. A proxy function for x is a function p of the risk factor space that is inexpensive to evaluate. Where lower and upper error bounds on the exact loss are known, they can be considered as functions of the risk factor space satisfying:

(2) \begin{equation}l(r) \leq x(r) \leq u(r).\end{equation}

The distribution function of X is denoted by F and its density function by f, if it exists, so that

(3) \begin{equation}F(s) = \mathbb{P}(X\leq s) = \int_{-\infty}^{s} f(t)\,\mathrm{d}t, \quad -\infty < s < \infty.\end{equation}

The function F is called the loss distribution. For a given percentile $\alpha\in(0,1)$ , the quantile $\xi_{\alpha}$ of the loss distribution is defined by the generalised inverse of F:

(4) \begin{equation}\xi_{\alpha} = F^{-1}(\alpha) \,:\!=\, \inf\{s : F(s)\geq \alpha\},\quad 0 < \alpha <1.\end{equation}

To avoid notational confusion, note that the symbol p is used only to refer to the proxy, whereas probability percentiles are given the symbol $\alpha$ .

The capital requirement is defined as the one in two hundred loss given by $-\xi_{0.005}$ . Note that a minus sign has been introduced so that the capital requirement is stated as a positive number. The value $\xi_{\alpha}$ is not easily calculable since f and F are not known precisely. Estimation techniques are therefore required. In the following, statistical estimation techniques involving Monte Carlo sampling are introduced.

2.2 Basic estimation

Distributional parameters of a random variable can be estimated from empirical samples. Estimators can target parameters such as the mean, variance, and quantiles. They are designed to converge to the parameter’s value as the sample size is taken ever larger. Moreover, it is a common phenomenon that the distribution of the error between the estimator and the parameter, when suitably scaled, converges towards normal for large sample sizes^{Footnote 5}. An elementary example of this is the central limit theorem (Chapter 8 of Stirzaker, Reference Stirzaker2003) where under mild conditions on the independent identically distributed random variables $X_i$ with mean $\mu$ and variance $\sigma^2$ :

(5) \begin{equation}\mathbb{P}\left( \frac{1}{\sigma \sqrt{n}}\sum_{i=1}^n (X_i-\mu) \leq s\right) \rightarrow \Phi(s) \,:\!=\, \int_{-\infty}^s \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}t^2}\,\mathrm{d}t,\quad \text{as }n\rightarrow \infty.\end{equation}

Equivalently, it may be said that the error between the estimator $n^{-1}\sum_{i=1}^n X_i$ and parameter $\mu$ is approximately normal and, when suitably scaled, converges in distribution^{Footnote 6}, written

(6) \begin{equation}\sqrt{n}\left( \frac{1}{n}\sum_{i=1}^n X_i\,\,-\mu \right)\xrightarrow{d} N\left(0,{\sigma^2}\right),\quad\text{as }n\rightarrow\infty.\end{equation}

When the distribution of the estimator (suitably scaled) converges to normal, it is said that the estimator is asymptotically normal. The difference between the estimator and the parameter is called statistical error. The standard deviation of the statistical error is called the standard error^{Footnote 7}. We discuss the estimation of statistical error in section 2.4. These quantities are of interest since they allow for the construction of confidence intervals for the parameter.

For a given sample size n, the kth largest sample is called the kth order statistic. Given n independent identically distributed random variables $X_i$ , the kth largest is denoted by $X_{k:n}$ . Order statistics are important for our framework since they are very closely related to the quantiles of the distribution. A simple case of asymptotic normality of order statistics is given by Corollary B, section 2.3.3 of Serfling (Reference Serfling2009), where it is shown that if F possesses a positive density f in a neighbourhood of $\xi_{\alpha}$ , and f is continuous at $\xi_{\alpha}$ , then

(7) \begin{equation}\sqrt{n}\left( X_{\lceil \alpha n\rceil:n}-\xi_{\alpha}\right) \xrightarrow{d} N\left(0,\frac{\alpha(1-\alpha)}{f^2(\xi_{\alpha})} \right), \quad\text{as }n\rightarrow\infty.\end{equation}

Estimators taking the form of a weighted sum of ordinals are referred to as L-estimators^{Footnote 8} and can be written

(8) \begin{equation}T_n = \sum_{i=1}^n c_iX_{i:n}\end{equation}

where $c_i$ are constants dependent on the percentile $\alpha$ and the number of scenarios n and vary according to the chosen estimator. Depending on the choice of weights, L-estimators can target different parameters of the distribution. An elementary example of estimating the mean $\mu$ of X. By placing $c_i=1/n$ in (8) and applying (6), noting that the sum of the ordinals is equal to the sum of the sample, we have $\sqrt{n}\left(T_n-\mu\right) \xrightarrow{d} N(0,\sigma^2)$ .

When the underlying distribution of X is known, the distribution and density functions of $X_{k:n}$ can be derived through elementary probability^{Footnote 9}:

(9) \begin{equation}F_{X_{k:n}}(s) = \sum_{i=k}^n \left(\begin{array}{l}n\\i\end{array}\right) F(s)^i(1-F(s))^{n-i}, \quad -\infty < s< \infty, \end{equation}

(10) \begin{equation} f_{X_{k:n}}(s) = \frac{n!}{(n-k)!(k-1)!} f(s)F(s)^{k-1}(1-F(s))^{n-k}.\end{equation}

The expected value of the kth order statistic is given by

(11) \begin{equation}\mathbb{E}\left[X_{k:n}\right] = \frac{1}{B(k,n-k-1)}\int_0^1 F^{-1}(t)t^{k-1}(1-t)^{n-k}\,\mathrm{d}t\end{equation}

where B(x,y) is defined in (16c).

Next we consider a special case of L-estimators designed to target quantiles. This is of particular interest since capital requirements are determined by quantile estimates.

2.3 Quantile L-estimators

In our application, we are concerned with the special case of L-estimators that target quantiles $\xi_{\alpha}$ , and importantly the special case of $\alpha=0.005$ representing the 1 in 200 year loss. We call these quantile L-estimators and write their general form as

(12) \begin{equation}\hat{\xi}_{\alpha} = \sum_{i=1}^n c_iX_{i:n}\end{equation}

where we use the hat notation $\hat{\xi}_{\alpha}$ to highlight that the estimator is targeting $\xi_{\alpha}$ . For a given empirical sample $\{x_i\}_{i=1}^N$ from x(R), write the ordered values by $\{x_{(i)}\}_{i=1}^N$ so that $x_{(i)}\leq x_{(i+1)}$ for $i=1,\ldots,N-1$ and the value of the statistic for the sample is written $\sum_{i=1}^n c_ix_{(i)}$ .

We give here two basic examples of quantile L-estimators. In practice, the choice of the estimator involves considerations of bias and variance, as well as the number of non-zero weights involved in its calculation.

2.3.1 Basic L-estimator

To understand the form of quantile L-estimators, consider a basic example where the quantile is approximated by a single order statistic:

(13) \begin{equation}\hat{\xi}_{\alpha} = X_{\lceil \alpha n \rceil:n}.\end{equation}

This can be seen as having the form of the L-estimator (12) by setting

(14) \begin{equation}c_i = \delta_{i,\lceil \alpha n\rceil}, \end{equation}

where $\delta_{i,j}$ is the Kronecker delta notation satisfying $\delta_{i,j}=0$ whenever $i\neq j$ and $\delta_{i,i}=1$ . From (7), we see that under mild assumptions the estimator converges to the quantile with asymptotically normal statistical error:

(15) \begin{equation}\sqrt{n}\left(\hat{\xi}_{\alpha} - \xi_{\alpha}\right) \xrightarrow{d} N\left(0, \frac{\alpha(1-\alpha)}{f^2(\xi_{\alpha})}\right)\quad \text{as }n\rightarrow\infty.\end{equation}

One reason to use the basic estimator is its simplicity. However, estimators with lower standard error may be preferred. An example is the following.

2.3.2 Harrell-Davis L-estimator

Harrell & Davis (Reference Harrell and Davis1982) introduced the quantile estimator that weights ordinals according to the target percentile $\alpha$ in a neighbourhood around the $n\alpha$ -th ordinal. It is given by

(16a) \begin{equation}\hat{\xi}_{\alpha} = \sum_{i=1}^n w_{i}X_{i:n}\end{equation}

with coefficients $w_{i}$ dependent on n and $\alpha$ given by

(16b) \begin{equation}w_{i} = I_{i/n}(\alpha(n+1),(n+1)(1-\alpha)) - I_{(i-1)/n}(\alpha(n+1),(n+1)(1-\alpha))\end{equation}

where $I_x(a,b)$ denotes the incomplete-Beta function

(16c) \begin{equation}I_x(a,b) = \frac{B(x;\,a,b)}{B(1;\, a,b)}, \quad B(x;\, a,b) = \int_0^x t^{a-1}(1-t)^{b-1} \mathrm{d}t, \quad x>0\end{equation}

formed from the Beta function $B(x;\, a,b)$ . When $x=1$ , write $B(x,y)\equiv B(1;\,x,y)$ .

Asymptotic normality of the Harrell-Davis estimator was established by Falk (Reference Falk1985) under general conditions on the distribution involving smoothness of F and positivity of f. A rate of convergence was also derived under somewhat stronger conditions. Numerical properties of the estimator were investigated by Harrell & Davis (Reference Harrell and Davis1982).

2.4 Bootstrap estimation of statistical error

Estimators with low statistical errors are desirable due to them yielding small statistical error bounds on capital estimates. In section 2.3.1, it is shown that in the case of the estimator using a single order statistic, the standard error is approximately $\sqrt{\frac{\alpha(1-\alpha)}{nf^2(\xi_{\alpha})}}$ . This analytic expression for the standard error requires knowledge of $f(\xi_{\alpha})$ , the density of the loss distribution at the quantile $\xi_{\alpha}$ . In practical examples, f is not known analytically at any quantile. We therefore need a method of approximating the statistical error of estimators.

Efron (Reference Efron1979) introduced a technique called the bootstrap^{Footnote 10} that can be used in practical situations to investigate the distribution of statistical errors when sampling from unknown distributions. This allows, for example, the standard error to be approximated, for example, the parameter $\sigma$ in (6), and the expression $\sqrt{\frac{\alpha(1-\alpha)}{nf^2(\xi_{\alpha})}}$ in (15).

Suppose there is a sample $\{x_i\}_{i=1}^n$ from an unknown distribution F. A statistic, say $\hat{\xi}_{\alpha}$ , has been calculated on the basis of the sample, and the question is how representative is it of $\xi_{\alpha}$ ?

Calculating many such samples from F to analyse statistical errors may be prohibitively expensive or even not possible in real-life empirical settings. The bootstrap technique instead uses resampling with replacement to mimic new independent samples. The technique involves drawing an independent sample of size n from $\{x_i\}_{i=1}^n$ with replacement, denoted with star notation $\left\{x_i^*\right\}_{i=1}^n$ . Resampling from $\{x_i\}_{i=1}^n$ is computationally cheap and so this can be repeated many times. Each of our new samples gives a realisation of the estimator $\hat{\xi}_{\alpha}^*$ , and we combine these to form an empirical distribution, from which the variance (and other properties) can be estimated. In applications, the bootstrap variance $\mathrm{var}\left(\hat{\xi}_{\alpha}^*\right)$ is often considered an estimator for the original variance $\mathrm{var}(\hat{\xi}_{\alpha})$ .

An algorithm for calculating the bootstrap estimate of standard errors in described in Algorithm 6.1 of Efron & Tibshirani (Reference Efron and Tibshirani1994). Analytical expressions for the bootstrap sample means, variances and covariances were established by Hutson & Ernst (Reference Hutson and Ernst2000). They show that the bootstrap sample is equivalent to generating a random sample of size n drawn from the uniform distribution U(0,1) with order statistics $(U_{1:n},\ldots,U_{n:n})$ and applying the sample quantile function $F_n^{-1}$ defined by

(17) \begin{equation}F_n^{-1}(u) = x_{(k)}\quad\text{where } (k-1)/n < u \leq k/n, k\in\mathbb{N}.\end{equation}

Denote the mean and variance of the ordinal $X_{k:n}$ by $\mu_k$ and $\sigma_k$ respectively, so that

(18) \begin{equation}\mu_k = \mathbb{E}(X_{k:n}), \end{equation}

(19) \begin{equation}\sigma_k^{2} = \mathrm{var}\left(X_{k:n}\right).\end{equation}

Theorem 1 of Hutson & Ernst (Reference Hutson and Ernst2000) establishes that the bootstrap estimators of $\mu_k$ and $\sigma_k$ denoted $\mu_k^*$ and $\sigma_k^*$ are given by

(20a) \begin{equation}\mu_k^* = \mathbb{E}\left(F_n^{-1}\left(U_{k:n}\right)\right) = \sum_{i=1}^n w_i(k)x_{(i)}, \text{where} \end{equation}

(20b) \begin{equation}w_i(k) \,:\!=\, I_{i/n}(k, n-k+1) - I_{(i-1)/n}(k, n-k+1), \end{equation}

and $I_x(a,b)$ denotes the incomplete-Beta function defined in (16c). From (10) and the form of $F_n$ , we see that when $\alpha(n+1)=k$ , the Harrell-Davis weights $w_i$ , defined in (16b), exactly coincide with the bootstrap weights $w_i(k)$ defined in (20b).

The bootstrap variance $\sigma_k^{*2}$ is established in Theorem 2 of Hutson & Ernst (Reference Hutson and Ernst2000):

(21) \begin{equation}\sigma_k^{*2} \,:\!=\, \mathrm{var}\left(F_n^{-1}\left(U_{k:n}\right)\right) = \sum_{i=1}^n w_i(k)\left(x_{(i)} - \mu_k^*\right)^2\end{equation}

where $w_i(k)$ is given by (20b).

An interpretation of (2.4) is that for a given loss scenario, the relative chance (in the sense of the bootstrap) amongst all ordinals that the i-th ordinal represents the true quantile ( $\alpha=k/n)$ is given by the weight $w_i$ .

2.5 Assumption on proxy errors

The term proxy error refers to approximations introduced by the use of proxies. For a given risk factor scenario r, we denote the output of a proxy model by p(r) and the exact loss calculation x(r). Fundamental to our analysis is the following assumption regarding knowledge of proxy error bounds.

Assumption 1. The proxy model is calculated with known proxy error bounds: given a risk factor scenario r, the (possibly unknown) exact loss scenario x(r) satisfies

(22) \begin{equation}l(r) \leq x(r) \leq u(r)\end{equation}

for some known values l(r),u(r).

In practical applications, Assumption 1 can be satisfied by choosing prudent error bounds, as may be found as part of the design and validation of the proxy model. Proxy design within life insurance capital applications is discussed in Caroll & Hursey (Reference Caroll and Hursey2011), Hursey (Reference Hursey2012), and Hursey et al. (Reference Hursey, Cocke, Hannibal, Jakhria, MacIntyre and Modisett2014). Androschuck et al. (Reference Androschuck, Gibbs, Katrakis, Lau, Oram, Raddall, Semchyshyn, Stevenson and Waters2017) suggest during design and validation to consider, amongst other measures, the maximum absolute error amongst in and out-of-sample testing across the whole distribution. Under this approach, prudent choices of error bounds are larger than the maximum error observed.

For a given sample of risk factors $\{r_i\}_{i=1}^n$ , the associated error sample, also called the residuals, is given by $\{p(r_i)-x(r_i)\}_{i=1}^n$ . The residuals may also be considered as a random sample from $p(R)-x(R)$ . The case when the error is bounded is covered by our framework. Suppose the residual, considered as a random variable, is bounded by some $\Delta>0$ :

(23) \begin{equation}p(R) - x(R) \sim E, \quad |E|\leq \Delta.\end{equation}

Then by placing $l(r)=p(r)-\Delta$ and $u(r)=p(r)+\Delta$ , the conditions of Assumption 1 are satisfied.

Assumption 1 asserts no distributional properties on the error sample other than that the residuals are bounded for a given risk factor by known values. Assumption 1 is also trivially satisfied in the situation where a proxy function p has a constant error bound $\Delta>0$ so that

(24) \begin{equation}p(r)-\Delta \leq x(r) \leq p(r) + \Delta.\end{equation}

Similar to case (23), this can be seen by setting $l(r)=p(r)-\Delta$ and $u(r)=p(r)+\Delta$ . Assumption 1 is however much more general since it allows the error interval to be arbitrarily large and to vary with risk factor.

Hursey et al. (Reference Hursey, Cocke, Hannibal, Jakhria, MacIntyre and Modisett2014) describe an error function, formed from fitting a function from the risk factor space to the residuals, to inform the design of the proxy function. This error curve cannot be used directly to form upper and lower error bounds since it is formed from least-squares fitting, and therefore represents an averaged error. However, Cardiel (Reference Cardiel2009) shows how least-squared fitting can be adapted to create regressions of the upper and lower boundaries of data sets. In this application, a regression of the upper boundary of the residuals could be used to propose an upper bound function u(r). Similarly, a regression of the lower boundary of the residuals can be used to propose a lower bound function l(r).

Whatever technique is used to pose error bound functions, it is still necessary to appropriately validate the functions before use in capital requirement modelling.

In the analysis of Murphy & Radun (Reference Murphy and Radun2021), an assumption that the error term is normally distributed is made so that

(25) \begin{equation}p(R) - x(R) \sim N(\mu, \sigma), \quad \mu\in\mathbb{R},\sigma>0.\end{equation}

This setting is excluded from our framework since it is not possible to assert any finite lower and upper bounds on the proxy error, and therefore requirement (22) of Assumption 1 is not satisfied. We note that deconvolution techniques may be applicable when an unbounded error model is assumed; see for example Ghatak & Roy (Reference Ghatak and Roy2018).

2.6 Assumptions on computational feasibility

Practical computational considerations motivate the use of approximation methods in capital requirement modelling. To reflect this, computational restrictions are introduced into the framework. By doing so, practical limitations faced by firms on the volume of computations they can perform are captured. Such limitations may arise for example through time constraints, computational budgets or practical considerations arising from the use of third party systems.

The notation $n\ll N$ is used to mean that N may be several orders of magnitude bigger than n. Androschuck et al. (Reference Androschuck, Gibbs, Katrakis, Lau, Oram, Raddall, Semchyshyn, Stevenson and Waters2017) cite a 2014 survey by Deloitte showing life insurance firms using between 20 and 1,000 exact loss calculations for calibration purposes. From this, we infer the likely number of feasible exact computations n for UK based life insurance firms to be at least in the range 20–1,000 and likely much higher given advances in computing capacity since 2014. They also state that to create a credible distribution of profit and loss, firms would expect to run hundreds of thousands, or even millions of scenarios, leading us to infer that N is of the order of one million. For our exposition, we consider $n=1,000$ and $N=1,000,000$ to be realistic values for some firms. Firms interested in understanding how proxy errors enter their capital estimates can adjust error bounds, n and N to their circumstances.

A consequence of Assumption 2 is that the residual error sets of the form $\{p(r_i)-x(r_i)\}_{i=1}^m$ , considered in section 2.5, are feasible to compute for $m=n$ but infeasible for $m=N$ .

Next, we consider the construction of computationally feasible bounds on quantile L-estimators.

2.7 Estimating quantiles of the loss distribution

Our interest is to estimate a quantile $\xi_{\alpha}$ , and especially $\xi_{0.005}$ , of the loss distribution without the presence of proxy errors. With a lemma about ordered sequences, we will apply quantile L-estimators to computable quantities to create lower and upper bounds on the L-estimator.

Lemma 1. Let $\{a_i\}_{i=1}^N\subset\mathbb{R}$ and $\{b_i\}_{i=1}^N\subset\mathbb{R}$ be arbitrary and unordered sequences of real numbers of length N satisfying the component-wise inequality:

(26) \begin{equation} a_i\leq b_i\quad \text{for}\quad i=1,\ldots,N. \end{equation}

Denote the ordered sequences by $\{a_{(i)}\}_{i=1}^N$ and $\{b_{(i)}\}_{i=1}^N$ so that $a_{(i)}\leq a_{(i+1)}$ and $b_{(i)}\leq b_{(i+1)}$ for $i=1,\ldots,N-1$ . Then, the following inequality holds between the ordinals:

(27) \begin{equation}a_{(i)} \leq b_{(i)}\quad \text{for}\quad i=1,\ldots,N.\end{equation}

Note that the permutation that sorts the two sequences is not assumed to be equal.

Lemma 1, along with assumptions on the existence of proxy error bounds (Assumption 1), establishes upper and lower bounds on L-estimators.

Lemma 2. Suppose there exists functions of the risk factor space $l,x,u\colon\mathbb{R}^m\mapsto\mathbb{R}$ satisfying

(28) \begin{equation}l(r) \leq x(r) \leq u(r)\end{equation}

for all $r\in\mathbb{R}^m$ (that is, we suppose Assumption 1 holds). Let $\{r_i\}_{i=1}^N$ be an independent sample from the risk factor distribution R and write $x_i \,:\!=\,x(r_i), l_i\,:\!=\,l(r_i)$ and $u_i\,:\!=\,u(r_i)$ for $i=1,\ldots,N$ . Then, whenever $c_i\geq 0$ for $i=1,\ldots,N$ , the value of the L-statistic $\hat{\xi}_{\alpha}$ given by

(29) \begin{equation}\hat{\xi}_{\alpha}\,:\!=\,\sum_{i=1}^Nc_i x_{(i)}\end{equation}

satisfies

(30) \begin{equation} \sum_{i=1}^N c_il_{(i)} \leq \hat{\xi}_{\alpha} \leq \sum_{i=1}^N c_i u_{(i)}\end{equation}

where $\{x_{(i)}\}_{i=1}^N$ is the ordered sequence of $\{x_i\}_{i=1}^N$ so that $x_{(i)}\leq x_{(i+1)}$ for all $i=1,\ldots N-1$ , and similarly, $\{l_{(i)}\}_{i=1}^N$ is the ordered sequence of $\{l_i\}_{i=1}^N$ , and $\{u_{(i)}\}_{i=1}^N$ is the ordered sequence of $\{u_i\}_{i=1}^N$ .

Under Assumption 2, the computation of the lower and upper bounds on $\hat{\xi}_{\alpha}$ in (30) is feasible.

Proof. By definition, $l_i=l(r_i)$ and $x_i=x(r_i)$ , and by (28) it follows that

(31) \begin{equation}l_i = l(r_i) \leq x(r_i) = x_i.\end{equation}

Therefore by writing $a_i\,:\!=\,l_i$ and $b_i\,:\!=\,x_i,$ we observe $a_i \leq b_i$ and so the conditions of Lemma 1 are satisfied to give

(32) \begin{equation}l_{(i)}\leq x_{(i)}.\end{equation}

Similarly, with $a_i\,:\!=\,x_i$ and $b_i\,:\!=\,u_i$ , it follows that $a_i\leq b_i$ by (28). Then, a direct application of Lemma 1 gives

(33) \begin{equation}x_{(i)} \leq u_{(i)}.\end{equation}

Since inequalities are preserved under multiplication of non-negative numbers, from (32) and (33) it follows that

(34) \begin{equation} \sum_{i=1}^N c_il_{(i)} \leq \sum_{i=1}^N c_i x_{(i)} \leq \sum_{i=1}^N c_i u_{(i)}.\end{equation}

The result now follows trivially from the definition of the L-statistic $\hat{\xi}_{\alpha}$ in (29).

If the upper and lower bound proxies have a constant error bound $\Delta>0$ , so that $p(r)-\Delta \leq x(r) \leq p(r)+\Delta$ for all risk factor scenarios r, and we assume the common case of the quantile L-estimator weights summing to one, $\sum_i c_i = 1$ , then a trivial consequence of Lemma 2 is that for a given risk factor scenario set $\{r_i\}_{i=1}^N$ with $p_i\,:\!=\,p(r_i)$ , we have

(35) \begin{equation}\left| \sum_{i=1}^N c_ip_{(i)} - \hat{\xi}_{\alpha} \right| \leq \Delta.\end{equation}

For both the Basic and Harrell-Davis estimators, the weights $c_i$ sum to 1, therefore in either of these cases the value of the exact statistic lies within $\pm \Delta$ of the statistic calculated with proxies: $\sum_i c_i p_{(i)}$ .

The assumption that the coefficients of the L-estimator are non-negative in Lemma 2, $c_i\geq 0$ , can be removed at the cost of accounting for how the sign impacts the upper and lower bounds:

(36) \begin{equation}\left.\begin{array}{l}c_i l_i\\ c_i u_i \end{array}\right\} \leq c_ix(r_i) \leq \begin{cases} c_i u_i, & c_i\geq 0, \\ c_i l_i, &c_i\leq 0. \end{cases}\end{equation}

The consideration of signs in (36) allows for estimators with negative weights, such as those targeting the inter-quartile range, to be considered.

Next, the question of how quantile estimates can be improved by targeted exact computation is considered.

2.8 Targeted exact losses

As the use of the heavy model to create exact loss scenarios over a large set of risk factors is assumed infeasible, we introduce the term targeted exact computation to emphasise when computations are performed with the heavy model on only a selected subset.

When exact computations are performed, the knowledge of the lower and upper proxy bounds can be updated. Consider a set of risk factor scenarios $\{r_i\}_{i=1}^N$ with lower and upper proxy bounds satisfying

(37) \begin{equation}l(r_i) \leq x(r_i) \leq u(r_i).\end{equation}

Suppose an exact calculation is performed for a subsect of scenarios $\mathcal{A}\subset\{1,2,\ldots,N\}$ . Then, writing $x_i \,:\!=\,x(r_i), l_i\,:\!=\,l(r_i)$ and $u_i\,:\!=\,u(r_i)$ for $i=1,\ldots,N$ , updated knowledge of lower and upper bounds can be indicated with a prime notation so that

(38a) \begin{equation}l_i \leq l^{\prime}_i= x_i = u^{\prime}_i \leq u_i\quad \text{ if }\quad i\in\mathcal{A},\end{equation}

(38b) \begin{equation}l_i = l^{\prime}_i \leq x_i \leq u^{\prime}_i = u_i\quad \text{ if }\quad i\not\in\mathcal{A}.\end{equation}

targeted exact computations improve our knowledge by removing proxy errors from the selected scenarios. By Assumption 2, it is possible to perform targeted exact computations on a set of risk factor scenarios $\{r_i\}_{i=1}^n$ whenever n is feasibly small. After a targeted exact computation has been performed, the updated ordered lower and upper bounds are denoted by $\{l^{\prime}_{(i)}\}_{i=1}^N,\{u^{\prime}_{(i)}\}_{i=1}^N$ .

In some circumstances, performing exact computations on arbitrary scenarios may not help improve our knowledge of specific loss scenarios. Consider, for example, the objective to find the value of the (unknown) kth largest loss scenario $x_{(k)}$ . Suppose there exists an upper bound $u_i$ with $u_i< l_{(k)}$ . Then, the scenario $x_i$ cannot correspond to the kth largest scenario since $x_i\leq u_i<l_{(k)}\leq x_{(k)}$ . Therefore, if $u_i<l_{(k)}$ then $x_i\neq x_{(k)}$ . Similarly, if $u_{(k)}< l_i$ then $x_i\neq x_{(k)}$ . Therefore, the index $i_k$ of the scenario $x_{i_k}$ corresponding to $x_{(k)}$ must satisfy

(39) \begin{equation}i_k \in \mathcal{A}_k\,:\!=\, \{1,2,\ldots,N\} \setminus \{j\,:\, u_j<l_{(k)} \text{ or }u_{(k)} < l_j , 1\leq j \leq N\}.\end{equation}

Note that

(40a) \begin{equation}u_j<l_{(k)} \text{ or }u_{(k)} < l_j \iff l_j\leq u_j<l_{(k)} \text{ or } u_{(k)}<l_j\leq u_j \end{equation}

(40b) \begin{equation}\qquad\qquad\quad\ \ \iff [l_j,u_j ]\cap [l_{(k)},u_{(k)} ]= \emptyset\end{equation}

and so equivalently

(41) \begin{equation}i_k\in \mathcal{A}_k=\left\{j:[l_j,u_j ]\cap [l_{(k)},u_{(k)} ]\neq \emptyset,1\leq j\leq N\right\}.\end{equation}

The indices $j\in\mathcal{A}_k$ therefore have corresponding intervals $[l_j,u_j]$ that must overlap, or touch, the interval $[l_{(k)},u_{(k)}]$ . At least one of these intervals must contain $x_{(k)}$ . If $|\mathcal{A}_k|=1$ , then $\mathcal{A}_k=\{i_k\}$ and the kth loss scenario has been identified. Otherwise, it is not known which intervals contain the kth loss scenario since the intervals are not necessarily pairwise disjoint, and so there may be many possible orderings of the exact losses.

An exact calculation of a scenario outside of this set does not provide any information about $x_{(k)}$ . This motivates the study of the set $\mathcal{A}_k$ in the context of targeted computation. The following result shows exactly how targeted exact loss computation over the set $\mathcal{A}_k$ can be used to remove proxy errors from basic quantile estimators.

Theorem 1. Let $\mathcal{A}_k$ be the set of indices j for which the interval $[l_j,u_j]$ of lower and upper proxy error bounds intersects the interval $[l_{(k)},u_{(k)} ]$ :

(42) \begin{equation}\mathcal{A}_k\,:\!=\,\left\{j:[l_j,u_j ]\cap [l_{(k)},u_{(k)} ]\neq \emptyset,1\leq j\leq N\right\}. \end{equation}

Suppose the exact loss $x_j$ is calculated for each $j \in \mathcal{A}_k$ and denote the updated knowledge of the lower and upper bounds of $\{x_i\}_{i=1}^N$ by $\{l^{\prime}_i\}_{i=1}^N$ and $\{u^{\prime}_i\}_{i=1}^N$ . Denote the sorted updated bounds by $\{l^{\prime}_{(i)}\}_{i=1}^N$ and $\{u^{\prime}_{(i)}\}_{i=1}^N$ . Then, the kth value of the sorted bounds are equal, and the kth exact loss (without proxy errors) $x_{(k)}$ is given by

(43) \begin{equation}x_{(k)} = l^{\prime}_{(k)}=u^{\prime}_{(k)}.\end{equation}

Assumption 2 implies that the computation of the exact losses is feasible if $|\mathcal{A}_k|\leq n$ .

Proof. Suppose $l_i\leq x_i\leq u_i$ for $i=1,2,\ldots,N$ , and let $i_k$ be an index satisfying $x_{i_k}=x_{(k)}$ . Following the equivalent formulation of $\mathcal{A}_k$ from (39), define the sets of indices

(44) \begin{equation}\mathcal{I} \,:\!=\, \{1,2,\ldots,N\},\end{equation}

(45) \begin{equation}\mathcal{L}_k \,:\!=\, \{i\in \mathcal{I}\,:\, u_i < l_{(k)}\},\end{equation}

(46) \begin{equation}\mathcal{R}_k \,:\!=\, \{i\in \mathcal{I}\,:\, u_{(k)} < l_i \}, \end{equation}

(47) \begin{equation}\mathcal{A}_k \,:\!=\, \mathcal{I}\setminus(\mathcal{L}_k\cup \mathcal{R}_k).\end{equation}

Note that $i_k\in\mathcal{A}_k$ . Suppose that the exact value of $x_i$ is calculated for each $i\in\mathcal{A}_k$ . Then, denoting the updated knowledge of upper and lower bounds with the prime notation, we have

(48a) \begin{equation} l_i\leq l^{\prime}_i = x_i = u^{\prime}_i\leq u_i\text{ whenever } i\in\mathcal{A}_k\end{equation}

(48b) \begin{equation} l_i=l^{\prime}_i \leq x_i \leq u^{\prime}_i = u_i \text{ whenever } i\not\in\mathcal{A}_k.\end{equation}

By Lemma 1 and (48) it follows that

(49) \begin{equation}l^{\prime}_{(i)}\leq x_{(i)} \leq u^{\prime}_{(i)} \text{ for all }i\in\mathcal{I}.\end{equation}

We wish to show that (49) holds with equality when $i=k$ . Consider first the inequality $l^{\prime}_{(k)} \leq x_{(k)}$ and suppose for a contradiction that $l^{\prime}_{(k)} < x_{(k)}$ .

If $l^{\prime}_{(k)}<x_{(k)}$ then there exists k distinct indices $\{i_1,\ldots,i_k\}\subset\mathcal{I}$ with $l^{\prime}_{i_k} < x_{(k)}$ . However, by the definition of $x_{(k)}$ , there exists at most $k-1$ distinct indices $\{j_1,\ldots,j_{k-1}\}\subset\mathcal{I}$ with $x_{j_k}<x_{(k)}$ . Since $l^{\prime}_{i}\leq x_{i}$ for any i, we have $\{j_1,\ldots,j_{k-1}\}\subset\{i_1,\ldots,i_k\}$ , and so there exists an index $i\in \{i_1,\ldots,i_k\}\setminus\{j_1,\ldots,j_{k-1}\}$ satisfying $l^{\prime}_{i}<x_{(k)}\leq x_i\leq u^{\prime}_i$ .

In particular, i is such that $l_{(k)}\leq l^{\prime}_{(k)} < x_{(k)} \leq u^{\prime}_i$ , and so $i\not\in \mathcal{L}_k$ . Also, $l^{\prime}_i < x_{(k)}\leq u^{\prime}_{(k)} \leq u_{(k)}$ , hence $i\notin \mathcal{R}_k$ . Therefore $i\in\mathcal{A}_k$ , and by (48a) we have $l^{\prime}_{i}=x_i=u^{\prime}_i$ . This contradiction, along with (49), proves $l^{\prime}_{(k)} = x_{(k)}$ .

Finally, by similar logic it holds that $x_{(k)}= u^{\prime}_{(k)}$ .

Theorem 1 establishes that we can identify loss ordinals exactly and therefore remove proxy errors from the basic quantile estimator whenever the computation is feasible.

It is a straightforward consequence that proxy errors can be removed from general L-estimators, when all non-zero weighted loss ordinals have been calculated exactly. To see this, consider a general L-estimator of the form (29):

(50) \begin{equation}\hat{\xi}_{\alpha}\,:\!=\,\sum_{i=1}^N c_i x_{(i)}\end{equation}

and define the set

(51) \begin{equation}\mathcal{B} \,:\!=\, \bigcup_{j: c_j>0} \mathcal{A}_j.\end{equation}

Theorem 1 shows that if the exact loss $x_i$ is calculated for every $i\in\mathcal{A}_j$ , then $x_{(j)}$ is known (without proxy errors). Clearly then, if the exact losses $x_i$ are calculated for all $i\in\mathcal{B}$ then $x_{(j)}$ is known exactly since $\mathcal{A}_j\subseteq\mathcal{B}$ .

However, $|\mathcal{B}|$ may be large and the calculation to remove proxy errors from general L-estimators infeasible. For example, in the case of the Harrell-Davis estimator, $c_i=w_i(k)>0$ for $i=1,\ldots,N$ , and therefore $\left|\mathcal{B}\right|=N$ , hence the exact calculation of all scenarios in this set is infeasible by Assumption 2. In the following, we motivate the introduction of sets $\mathcal{B}_{k,\varepsilon}\subseteq\mathcal{B}$ , which are sufficiently small to enable exact losses $x_i$ to be calculated for all $i\in\mathcal{B}_{k,\varepsilon}$ (i.e. $|\mathcal{B}_{k,\varepsilon}|\leq n$ ), but are chosen such that potential proxy errors on losses $x_i$ with $i\notin\mathcal{B}_{k, \varepsilon}$ have a small, measurable impact (controlled by the parameter $\varepsilon$ ) on approximations to certain calculations, for example, $\hat{\xi}_{\alpha}$ , $\sigma^*_k$ .

Consider the case where values required in calculations are of the form $w_i(k) x_{(i)}$ . By Lemma 1, the value is bounded: $|w_i(k)x_{(i)}|\leq w_i(k)\max(|l_{(i)}|,|u_{(i)}|)$ . If $|\mathcal{A}_i|\leq n$ , the computation to remove proxy errors from $x_{(i)}$ is feasible (Assumption 2 and Theorem 1) and so $w_i(k)x_{(i)}$ can be calculated exactly. Alternatively, only the bound is known. Since the weights $\{w_i(k)\}_{i=1}^N$ decay rapidly to zero away from $i=k$ , we identify the set

(52) \begin{equation}\mathcal{B}_{k, \varepsilon} \,:\!=\, \bigcup_{j: w_j(k)>\varepsilon} \mathcal{A}_j\end{equation}

as having potential practical uses. The decay of the weights may enable a choice of $\varepsilon$ such that $|\mathcal{B}_{k,\varepsilon}|\leq n$ making the removal of proxy errors from losses $x_i$ with $i\in\mathcal{B}_{k,\varepsilon}$ feasible and may be such that the bounds $|w_i(k)x_{(i)}|\leq \varepsilon\max(|l_{(i)}|,|u_{(i)}|)$ for $i\notin\mathcal{B}_{k, \varepsilon}$ are sufficiently small to support the applicability of approximations.

In section 4, we show through a prototypical numerical example that in some circumstances it may be appropriate for practitioners to approximate the bootstrap estimate of standard errors, (21), using

(53) \begin{equation} \hat{\sigma}_{k,\varepsilon}^{*2}\,:\!=\,\sum_{j:w_j(k)>\varepsilon} w_j(k)\left(x_{(j)} - \sum_{i:w_i(k)>\varepsilon} w_i(k)x_{(i)}\right)^2,\end{equation}

for some $\varepsilon>0$ , where in particular $|\mathcal{B}_{k,\varepsilon}|$ is sufficiently small for exact loss calculations to be performed over the entire set.

Next, an explanatory example is used to present a calculation recipe for removing proxy errors from loss ordinals.

4.1 Computational assumptions

The capital requirement is assumed to be of the form $-x_{(k)}$ making the evaluation of $x_{(k)}$ the primary objective of the study. By Theorem 1, it is sufficient to perform $|\mathcal{A}_k|$ targeted exact computations to find the kth loss ordinal $x_{(k)}$ exactly. The secondary objective is to find an approximation to the bootstrap estimate of standard error of $x_{(k)}$ – the expression $\hat{\sigma}_{k,\varepsilon}^{*}$ is posed as a candidate for this. As discussed in section 2.8, it is sufficient to perform $|\mathcal{B}_{k,\varepsilon}|$ targeted exact computations to compute the value $\hat{\sigma}_{k,\varepsilon}^{*}$ .

In what follows, we study the computational feasibility of eliminating proxy errors over a range of N and over a range of proxy error bounds $\Delta$ . In section 2.6, it is established that a plausible value for the number of proxy loss scenarios N is given by $N=1,000,000$ , and a plausible number of feasible exact loss calculations n is given by $n=1,000$ . Therefore, for demonstration purposes, it is assumed that it is feasible to perform targeted exact computations of scenarios in the sets $\mathcal{A}_k$ or $\mathcal{B}_{k,\varepsilon}$ whenever $\left|\mathcal{A}_k\right|\leq1,000$ or $\left|\mathcal{B}_{k,\varepsilon}\right|\leq1,000$ , respectively.

4.2 Conclusions of numerical experiment

Figures 2 and 3 show example calculations based on the model described in section 4.3. We conclude that, in the prototypical setting:

• It may be possible to eliminate proxy errors from capital estimates for large N if proxy error bounds are sufficiently small.

• In this example, proxy errors can be eliminated from the capital estimate when $N=1,000,000$ if the additive error bound is approximately £ 60m.

• The number of calculations required to approximate the bootstrap estimate of the standard error of the basic L-estimator to within $0.1\%$ is small relative to N.

• For a fixed number of targeted exact computations, there is a choice of N which minimises error bounds (proxy plus statistical) around the true loss.

• When the number of targeted exact computations n is computationally limited, increasing the Monte Carlo sample size N may not improve total error bounds due to the reintroduction of proxy errors.

Figure 2 The method of targeted exact computation is applied to a prototypical loss distribution. Panel A: Number of exact calculations sufficient to remove proxy errors from 0.5th percentile basic L-estimator statistics, given by $\left|\mathcal{A}_{0.005\times N}\right|$ , with N proxy scenarios and proxy error bound $\Delta$ . The contours indicate feasible combinations of $\Delta$ and N for a given number of targeted exact computations n. The region of infeasible combinations for $n>1,000$ is shaded. Panel B: Confidence interval due to statistical and proxy errors for the true 0.5th percentile loss, plotted for varying N for fixed $\Delta=100$ and $n=1,000$ . Dashed line shows maximum N at which proxy errors are removed after 1,000 targeted exact computations. In both panels, loss is a normal inverse Gaussian (NIG) random variable with parameters $a=0.6/750$ , $b=-0.2/750$ , $\delta=750$ , and $\mu=200$ .

Figure 3 Panel A: Sufficient number of exact calculations, given by $\left|\mathcal{B}_{k,\varepsilon}\right|$ , to eliminate proxy errors from the bootstrap approximation of standard error, when the loss is a normal inverse Gaussian (NIG) random variable with parameters $a=0.6/750$ , $b=-0.2/750$ , $\delta=750$ , and $\mu=200$ , and the 0.5th percentile is calculated using the basic L-estimator with N scenarios and proxy error bound $\Delta$ . For a given N we choose $\varepsilon=\varepsilon(N)$ such that $\sum_{j:w_j(k)>\varepsilon(N)}w_j(k)\geq0.999$ . The contours indicate feasible combinations of $\Delta$ and N for a given number of targeted exact computations n. The region of infeasible combinations for $n>1,000$ is shaded. Panel B: Percentage error of the approximation of the bootstrap estimate of standard error, (53), relative to the actual bootstrap estimate (21) given by $(\sigma^*_k-\hat{\sigma}_{k,\varepsilon}^*)(\sigma^*_k)^{-1}$ . Parameter $\varepsilon=\varepsilon(N)$ in (53) is chosen as in Panel A.

In Figure 2 Panel A and in Figure 3 Panel A, the boundary of feasible and infeasible combinations of N and $\Delta$ are shown, when the feasible number of exact computations n satisfies $n\leq1,000$ . In both cases, the number of targeted scenarios increases with N and $\Delta$ . Smaller error bounds achieve both the elimination of proxy errors from capital estimates, and the possibility of reducing statistical error by increasing the sample size of proxy loss scenarios.

Figure 2 Panel B shows confidence intervals around the true loss from proxy and statistical errors for a fixed proxy error $\Delta=100$ and fixed number of targeted exact computations $n=1,000$ . The dashed vertical line shows that proxy errors are introduced when $N>600,000$ . As N increases, the proxy error initially grows slower than the statistical error reduces, and so the confidence interval becomes tighter. Past an optimal point of $N\approx1,250,000$ , the proxy error grows more quickly and the confidence interval increases. This shows that when the number of possible targeted exact computations is limited (computationally), it may be sub-optimal to increase the number of Monte Carlo runs N past a certain value. On the other hand, for a fixed N, the full number of targeted exact computations $\left|\mathcal{A}_{0.005\times N}\right|$ may not be required to achieve near-optimal confidence bounds.

Figure 3 Panel B plots the percentage error of using $\hat{\sigma}_{k,\varepsilon(N)}^{*}$ , from (53), to approximate the bootstrap estimate of the standard error, $\sigma_{k,}^{*}$ , given in (21), for varying N. The percentage error is below $0.1\%$ with $N=300,000$ simulations. By comparison with Figure 2 Panel A and Figure 3 Panel A, in this example it is computationally feasible to eliminate proxy errors from capital estimates and approximate the bootstrap estimate of the standard error to within $0.065\%$ with a fixed proxy error $\Delta=100$ and fixed maximum number of targeted exact computations $n=1,000$ .

4.3 Distributional assumptions

Define the prototypical loss random variable X as being a normal inverse Gaussian (NIG) random variable, introduced in Barndorff-Nielsen (Reference Barndorff-Nielsen1997), having a density function given by

(60) \begin{equation}f(t)=\frac{a}{\pi} \mathrm{exp}\left(\delta\sqrt{a^2-b^2}+b(x-\mu)\right)\frac{K_1\left(a\delta\sqrt{1+\frac{(x-\mu)^2}{\delta^2}}\right)}{ \sqrt{1+\frac{(x-\mu)^2}{\delta^2}}}\end{equation}

for $a,b, \delta, \mu,t\in\mathbb{R}$ , with $0\leq|b|<a$ , and $\delta>0$ , where $K_1(\cdot)$ is a modified Bessel function defined by

(61) \begin{equation}K_1(z) = z \int_1^\infty e^{-zt}\sqrt{t^2-1}\,\mathrm{d}t,\quad z>0.\end{equation}

We choose $a=0.6/750$ , $b=-0.2/750$ , $\delta=750$ , and $\mu=200$ , so that the 0.5th percentile loss is $\xi_{0.005}=-4465.22$ , see Figure 4 for a plot of the distribution.

Figure 4 Prototypical loss distribution chosen to exhibit large tails and skew taking the form of a normal inverse Gaussian (NIG) distribution. The density and distribution functions are shown in Panels A and B, corresponding to NIG parameters: $a=0.6/750$ , $b=-0.2/750$ , $\delta=750$ , and $\mu=200$ .

In order to have lower and upper bounds that satisfy Assumption 1, suppose for any given risk factor scenario r,

(62) \begin{equation}l(r) = x(r) - \Delta, \quad u(r) = x(r) + \Delta,\end{equation}

for some $\Delta>0$ which we choose to vary as part of the numerical investigation. Note, by construction, $l(r) \leq x(r) \leq u(r)$ for all r.

For a given N, we generate a sample $\{x_i\}_{i=1}^N$ using the pseudo-random NIG generator of Scipy, from which we also calculate $\{l_i\}_{i=1}^N$ , $\{u_i\}_{i=1}^N$ , $\{l_{(i)}\}_{i=1}^N$ , and $\{u_{(i)}\}_{i=1}^N$ . The sufficient number of calculations required in Theorem 1, $\left|\mathcal{A}_{0.005\times N}\right|$ , is then found for different values of $\Delta$ .