{'fraction_non_alphanumeric': 0.04523876012287071, 'fraction_numerical': 0.016419994414967887, 'mean_word_length': 3.8682254848649626, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computation of expected shortfall by fast detection of worst scenarios

We consider a multi-step algorithm for the computation of the historical expected shortfall such as defined by the Basel Minimum Capital Requirements for Market Risk. At each step of the algorithm, we use Monte Carlo simulations to reduce the number of historical scenarios that potentially belong to the set of worst scenarios. The number of simulations increases as the number of candidate scenarios is reduced and the distance between them diminishes. For the most naive scheme, we show that the L p-error of the estimator of the Expected Shortfall is bounded by a linear combination of the probabilities of inversion of favorable and unfavorable scenarios at each step, and of the last step Monte Carlo error associated to each scenario. By using concentration inequalities, we then show that, for sub-gamma pricing errors, the probabilities of inversion converge at an exponential rate in the number of simulated paths. We then propose an adaptative version in which the algorithm improves step by step its knowledge on the unknown parameters of interest: mean and variance of the Monte Carlo estimators of the different scenarios. Both schemes can be optimized by using dynamic programming algorithms that can be solved off-line. To our knowledge, these are the first non-asymptotic bounds for such estimators. Our hypotheses are weak enough to allow for the use of estimators for the different scenarios and steps based on the same random variables, which, in practice, reduces considerably the computational effort. First numerical tests are performed.


Introduction
The Basel Minimum Capital Requirements for Market Risk [4] has brought two main changes in the way that investment banks need to compute their capital requirements. Expected Shortfall (ES) replaces Value at Risk (VaR) as the main risk indicator for the computation of capital requirements. The advantages of ES over VaR have been brought forward in Artzner et al. [2], and Expected Shortfall is now considered by most researchers and practitioners as superior to VaR as a risk measure, because it respects the subadditivity axiom, see [1,2,23]. The second main change is that the number of required daily computations of ES has been multiplied by a factor of up to 90. Where banks used to need to compute one VaR per day, they now need to compute up to three ES per liquidity horizon and risk class, as well as three ES per liquidity horizon for all risk classes combined. The above has triggered several works on the fast computation of ES.
Mathematically, if V is a random variable modeling the level of loss 1 of a portfolio that will be known at a future time, and 0 < α < 1, the expected shortfall of level α ∈ (0, 1) is defined by where VaR γ is the Value at Risk at level γ, i.e.
VaR γ (V ) := max {x ∈ R : P[V ≥ x] > γ}. (2) Nearly all of the literature concentrates on studying the ES by using parametric, non-parametric or semiparametric approaches to approximate the distribution of V based on historical data. See in particular [9,11,12,15,16,19,20,21,24,25,26,27]. Another approach consists in using the fact that V is the risk neutral value of a book, and therefore of the form E[P |S] in which S is a random variable associated to market parameters and P represents the future (discounted) payoffs of the book. This suggests using a nested Monte Carlo approach : simulate a set of values in the distribution of S (outer scenarios), and, for each simulation of S, compute a Monte Carlo estimator of E[P |S] by using simulations in the conditional distribution (inner scenarios). This is for instance the approach of [8,13]. But, as defined in the regulatory document of Basel [4], the expected shortfall is based on n s = 253 scenarios of market parameters s = (s i ) i≤ns that are generated in a deterministic way. Therefore, S is just uniformly distributed in the sequence s and there is no need for simulating outer scenarios. Since V is defined by a pricing formula E[P |S] that is fully described by the value of S, there is also no room for approximating the law of V based on historical data, if we are only interested by the requirements of [4]. The only issue is to compute in an efficient way the loss impacts (µ i ) i≤ns of the book, in which s 0 is the current value of the market parameters, and then compute the average over the n w = 6 worst impacts, say if, for ease of notations, we assume that Methods that are in line with the above have also been studied, in particular in [17,22] in which the authors define a distance on the space of scenarios induced by the distance between their risk factors. Starting with the original outer-level scenarios (called "prediction points"), they determine "design points" that are included in their convex hull. Inner-level paths are simulated in order to evaluate the portfolio value at the design points. These values are then used to establish a metamodel of the portfolio price with respect to the risk factors, and this metamodel is then used to select among the prediction points those that are most likely to be part of the worst scenarios set. They are then added to the design points, and evaluated by using inner-level simulations, after which the metamodel is updated. These methods are very smart but neglect one important point for practitioners: the cost of launching a pricer is high, as it typically entails instanciating thousands of objects at initialization, as well as volatility surface calibrations and sometimes even graphical interfaces. Furthermore, these pricers usually do not have the flexibility to add dynamically, at each inner-level pricing, new paths to a given scenario. Therefore, we do not allow ourselves to adapt our strategies at such a level of granularity. Instead, we will consider strategies that only entail L-levels of sets of simulations, where L is typically quite low, so as not to pay too many times the overhead of launching the pricer and/or calibrating the required volatility surfaces. We also do not use any concept of distance between scenarios induced by their risk factors. Although this enables [17] and [22] to obtain better empirical convergence rates, we see at least one problem with this approach: at the scale of a bank, the space of risk factors is both of a very high dimension (a few thousands) and with a very complex geometry (the payoffs of the portfolio's derivative products are usually non-convex, and path-dependent), so that it is very difficult to establish a model describing the proximity of scenarios in a robust way.
We thus study a relatively simple procedure that also has the advantage of allowing us to establish non-asymptotic bounds on the L p -error of our estimator, in the spirit of the simplest ranking by mean procedures, see e.g. [3,5,6,14]. It consists in using a first set of simulated paths to provide a crude estimation of the impact factors µ i . These first estimators are ranked to select the q 1 < n s outer-level scenarios with the highest estimated impact values. Then, only the impact values of these q 1 pre-selected scenarios are estimated again by using the previous estimators together with a new set of simulated paths. Among these new estimators we select the scenarios with the q 2 < q 1 highest estimated impact factors. And so on. After L ≥ 2 steps, L being small in practice, we just keep the mean of the six highest estimated impacts. The rationale behind this is that a first crude estimation should be sufficient to rule out a large part of the scenarios from the candidates of being in the 6 worst ones, because the corresponding values should be far enough. While the number of candidates reduces, one can expect that the differences between the corresponding impacts diminish as well and that more Monte Carlo simulations are needed to differentiate them. Under an indifference zone hypothesis, similar to the one used in the above mentioned paper, and a sub-gamma distribution assumption, the convergence is exponential in the number of simulations used at the different steps and of order 1/2 in the total number of simulations. See Proposition 2.2 and Corollary 2.3 below. The optimal number of additional paths that should be used at each step to minimize the strong estimation error, given a maximal computational cost, can be determined by a simple dynamic programming algorithm, that can be solved off-line, see Section 2.4. In theory, this requires the a priori knowledge of the means and covariances of our estimators, which are obviously not known in practice. However, one can easily define a version based on a robust specification of the error. One can also take advantage of the different simulation sets to improve our prior on the true hidden parameters. This leads to a natural adaptative algorithm, see Section 3, for which convergence is also proved, see Proposition 3.3. Estimating the optimal policy associated to this adaptative algorithm is costly but can be done off-line by using a neural network approximation combined with a backward dynamic programming algorithm. We explain how this can be done in Section 3.3 (further details are in Appendix B). The rest of the paper is organized as follows. Section 2 is dedicated to the most naive deterministic algorithm. In particular, Section 2.5 gives a very easy to use two levels algorithm for the case where the impacts decrease linearly in the scenarios' rank order. The adaptative version of the algorithm is presented in Section 3. Finally, we perform first numerical tests in Section 4.

Algorithm with a deterministic effort allocation
In this section, we describe the simplest version of the algorithm. It uses a pre-defined deterministic number of simulations. We establish a strong L p -error bound and discuss several ways of choosing the optimal strategy for minimizing this error.

The algorithm
From now on, we assume that E[P |S = s 0 ] is known perfectly and set it to 0 for ease of notations. As explained above, the algorithm relies on the idea of selecting progressively the scenarios that will be used for the computation of the Expected Shortfall. Namely, let P |s := (P |s i ) i≤ns be a n s -dimensional random variable such that each P |s i has the law of P given S = s i . We first simulate independent copies (P 1 j , . . . , P ns j ) j≥1 of P |s and compute the Monte Carlo estimators of E[P |S = s i ], i ≤ n s : for some N 1 ≥ 1. Among these random variables, we then select the ones that are the most likely to coincide with the worst scenarios s 1 , . . . , s nw , for some 1 ≤ n w < n s . To do this, one considers the (random) permutation m 1 on [[1, n s ]] such that the components of μ m1(i) 1 i≤ns are in decreasing order: and only keep the indexes (m 1 ( )) ≤q1 of the corresponding q 1 ≥ n w highest values, i.e. the indexes belonging to We then iterate the above procedure on the scenarios in I 1 and so on. Namely, we fix L ≥ 1 different thresholds (q ) =0,...,L−1 such that n w =: q L−1 ≤ · · · ≤q 0 := n s .
Assuming that I −1 is given, for some 1 ≤ − 1 ≤ L − 1, we compute the estimators 2 and only keep the elements in for the next step. If = L, we just compute the final estimator of the ES given by Note that only the L − 1-first steps are used to select the worth scenarios, the step L is a pure Monte Carlo step. Again, the general idea is to reduce little by little the number of candidate scenarios to be part of the worst ones. As the number of candidates diminishes, one increases the number of simulated paths so as to reduce the variance of our Monte Carlo estimators and be able to differentiate between potentially closer true values of the associated conditional expectations.
Remark 2.1. Note that, given j, we do not assume that the P i j , i ≤ n s , are independent. The simulations associated to different scenarios are in general not independent. Moreover, theμ i , ≤ L, use the same simulated paths, only the number of used simulations changes. Both permit to reduce the computational cost, by allowing the use of the same simulations of the underlying processes across scenarios and steps.

General a-priori bound on the L p error
In this section, we first provide a general L p estimate of the error. A more tractable formulation will be provided in Corollary 2.3 under an additional sub-gamma distribution assumption.
Before providing the proof of this general estimate, let us make some comments. The last two terms in (10) are natural as they are due to the Monte Carlo error made on the estimation of the various conditional expectations that can enter, after the (L−1)-levels selection procedure, in the estimation of ES. Note that it corresponds to the estimation errors using the cumulated number of Monte Carlo simulations N L−1 of step L−1 and the number N L −N L−1 of simulations used only for the last step. In practice, these numbers should be sufficiently large. The first term involves the quantities max Each term corresponds to the situation in which an element i ∈ [[1, n w ]] gets out the set of selected indexes I exactly at the -th step. In the worst situation, it is replaced by an element of index k larger than q and this can happen only ifμ k >μ i . The probability of this event is controlled by the number of Monte Carlo simulations N used at the step but also by the distance between the two scenarios. More specifically, for small, one expects that P[μ k >μ i ] is small because the law of P |s k is concentrated far away from where the law of P |si is. This quantity potentially increases with , as we reduce the number of selected indexes. This should be compensated by an increase in the number of used Monte Carlo simulations. Otherwise stated, we expect to balance the various terms of (10) by considering a suitable increasing sequence (N ) ≤L .
Obviously, (10) implies that the algorithm converges as N → ∞ for all ≤ L, see Proposition 3.3 below for a proof in a more general framework.
Proof of Proposition 2.2. We split the error into a permutation and a Monte Carlo error: Let us first look at the second term which corresponds to a Monte Carlo error. We have To discuss the first term in the right-hand side of (11), the permutation error, let us first define S q such that |A| = q . Then, it follows from (4), Hölder's inequality and (12) that

Error bound for Sub-Gamma distributions
To illustrate how the general error bound of Proposition 2.2 can be used in practice to decide of the sequence (q , N ) , we now consider the case where the components of P |s have sub-gamma distributions, and apply Bernstein's inequality in (10), see e.g. [7,Chapter 2]. This requires the following assumption.
Assumption 1. There exists c ∈ R + such that the random variables Z[i, k] := (P |s i − µ i ) − (P |s k − µ k ), i, k ≤ n s , satisfy Bernstein's condition : From now on, we shall assume that the constant c is known. It can usually be estimated in practice.
Corollary 2.3. Assume that Assumption 1 holds. Then, for all p ≥ 1, in which and C p,c := 4 p Γ (p) (15) where Γ is the Gamma function defined by The upper-bound of Corollary 2.3 has two advantages on Proposition 2.2. First, the dependence on (q , N ) ≥0 is more explicit. It depends on unknown quantities, but we can estimate (at least rough) confidence intervals for them, see e.g. Section 2.4 below. Second, as we will see in the next section, it allows one to define a tractable deterministic optimal control problem satisfying a dynamic programming principle, or even simple heuristics (see Section 2.5), to select an appropriate sequence (q , N ) ≥0 .
Proof of Corollary 2.3. The first term in (14) is an upper-bound for the first term in the right-hand side of (10), see [7, Theorem 2.1]. As for the two other terms in (10), we use the usual argument, for i ≤ n s , and then appeal to [7, Theorem 2.1] again to deduce that Remark 2.4. If the (μ i ) i≤ns and (δμ i ) i≤ns are Gaussian, which is the case asymptotically, then the bound of Corollary 2.3 remains valid with c = 0. This fact will be used later on for simplifying our numerical algorithms.

Optimal a-priori allocation by deterministic dynamic programming based on fixed a-priori bounds
Given N := (N ) 0≤ ≤L and q = (q ) 0≤ ≤L−1 , the total computation cost is with the convention N 0 := 0. Let N denote the collection of non-decreasing sequences N := (N ) 0≤ ≤L with values in N such that N 0 = 0, and let Q denote the collections of non-increasing sequences 4 q = (q ) 0≤ ≤L with values in [[n w , n s ]] satisfying (5). In this section, we fix a total effort K > 0 and recall how F p (q, N ), as defined in (14), can be minimized over the collection A of sequences (N, q) ∈ N × Q satisfying C(N, q) ≤ K by using a standard dynamic programming approach. where and defineF Then, the dynamic programming principle implies that This reduces the search for an optimal selection of (N ,q) to L − 1 one-step optimization problems, which is much simpler to solve than the optimization problem associated to the left-hand side of (13).
In practice, the exact values of (µ i , σ 2 i ) i≤ns and (σ 2 ik ) i,k≤ns are not known. However, one can consider robust versions of the above. For instance, if we know that there exists some (δ q , δ q ) q≤ns and σ 2 such that 4 We write (q ) 0≤ ≤L for convenience also q L will never play any role. 5 In the following, we only write A(0) for = 0 as it does not depend on (q,N ).
then one can similarly minimize the upper-bound of F p defined as This corresponds to a worst case scenario, when only the a priori bounds (δ q , δ q ) q≤ns and σ 2 are known. In the above, one can also impose that q takes values in a given subset of Q of Q. In this case, we will only need to know (δ q , δ q ) q∈Q .
We refer to Section 4 below for numerical tests that show that such an algorithm seems to perform pretty well. Note that the optimization can be done off-line.

Simplified 2-levels algorithm for a linear indifference zone's size
Inspired by [3,5,6,14], we assume here that we know the value of a constant δ 0 > 0 such that the impacts of the n w worst scenarios have values that are separated by at least (k − n w )δ 0 from the k-th worst scenario, for k > n w : To illustrate this, we plot on Figures 1-4 the curves k → |µ nw − µ k | for different formerly used test books of Natixis. We see that they are more flat on the interval [100, 120], so that a rather conservative value would be the minimum (over the different books) of (µ 100 − µ 120 )/20. Another choice in practice could be to take the ratio (µ nw − µ 100 )/(100 − n w ) which amounts to considering only the first part of the curve, and neglecting points that are anyway far from the worst scenarios. Figure 1: k → |µ nw − µ k | for Book #1 Figure 2: k → |µ nw − µ k | for Book #2 We now consider a simplified version of the algorithm of Section 2.1 where we only do one intermediate "fast pricing" (meaning N 1 rather small) and one final "full pricing" (meaning N 2 large). In theory, this corresponds to L = 3 with q 2 = n w , δN 3 = 0 and δN 2 → ∞. As δN 2 → ∞, the second and third terms in (14) vanish, as well as the component of the first term corresponding to = 2. We therefore neglect them. In practice, we only take N 2 large enough (and given) from the point of view of the bank, and minimize over (q 1 , N 1 ) the remaining term in (14): in whichσ is estimated to be as in (17), under the computation cost constraint For N 1 (or K) large enough, the condition (18) leads to minimizing over The optimal q * 1 can then be found easily by performing a one-dimensional numerical minimization. Upon replacing n s −q 1 by n s in the denominator of the exponential term, which provides a further upper-bound, the optimum can even be computed explicitly, see Appendix A. This provides a very easy to use algorithm.
Considering the case p = 1, let us now perform first numerical tests to see if the proxy based on h 1 0 is far from F ∞ 1 . We use the parameters of Tables 1 and 2 below and µ i = −iδ 0 , i ≤ n s , where δ 0 := (µ nw − µ 100 )/(100 − n w ) for the µ i s of Figure 6. In particular, we take c = 0, see Remark 2.4.
In Figure 5, the two increasing curves show the optimum q * 1 (right axis) as found when applying the deterministic dynamic programming algorithm (dashed line) of Section 2.4 associated to the real sample book curve of Figure 6, and the heuristic (solid line) based on (19). The two decreasing curves show the corresponding F ∞ 1 (q * 1 ) (left axis) found when applying the deterministic dynamic programming algorithm (dotted line) and the heuristic (dashdot line). We see that the heuristic and the real minimizer are extremely close. The noise in the lines associated to the dynamic programming algorithm are due to grid effects.     Although the true value θ • = (µ • , Σ • ) of the vector of means and of the covariance matrix of P |s are unknown, we can set on it a prior distribution, e.g. based on previous Monte Carlo experiments, rather than just working on robust bounds as in the end of Section 2.4. Since the estimation of ES uses Monte Carlo simulations of P |s , the knowledge of these quantities can be improved along the different steps of our estimation procedure. This suggests an adaptative algorithm for the optimization of the numerical effort allocation, in which we learn progressively the true value of these parameters, or part of them. From now on, we therefore view the true value of the parameters as a random variableθ := (μ,Σ) on which a prior law ν 0 is set. At each step , new Monte Carlo simulations will allow us to update this prior, and our strategy for the next steps accordingly.

Error bounds and convergence for predictable strategies
Let us first adapt the proof of Proposition 2.2 and Corollary 2.3 to the case where controls are not deterministic but stochastic processes. Given a stochastic process α with values in Q × N , we set (q α , N α ) := α where q α and N α are respectively Q and N -valued. We then defineμ α = (μ α ) ≤L , (I α , m α ) ≤L as in Section 2.1 except that we seeμ α as a q α -dimensional random variables with entries given by (μ α,i ) i∈I α . We use the same convention for δμ α , recall (9). We say that α is admissible if it is . We call A ad the collection of such processes. Then, one defines The true value of the expected shortfall is now also written as a random variable in whichm is the random permutation such that We let M be a collection of laws on R ns × S ns , where S ns denotes the collection of covariance matrices of size n s . Given ν ∈ M, we denote by E ν the expectation operator given thatθ admits the law ν. When ν is a Dirac mass, we retrieve the situation of Section 2 (up to re-ordering in a deterministic way the components of µ). We first provide a natural extension of Proposition 2.2.
with the convention max ∅ = 0 and in whichm α −1 is defined asm but on the subset I α −1 instead of . Proof. We proceed as in the proof of Proposition 2.2 to obtain that We define k α as k in the proof of Proposition 2.2 for the strategy α, with R replaced by R α := Remark 3.2. Note that, when α is deterministic and ν is concentrated on a Dirac, the right-hand side of (20) is bounded from above by which coincides with the bound of Proposition 2.2 The above guarantees the convergence of the algorithm.
Proposition 3.3. Let (K n ) n≥1 ⊂ N * be a sequence converging to infinity and let (α n ) n≥1 be a sequence in A ad such that C(q α n , N α n ) ≤ K n for each n ≥ 1. Assume further that min 1≤ ≤L N α n → ∞ a.s. Let ν be concentrated on the Dirac mass on θ • . Then, Proof. It suffices to use the fact that, for some C p > 0, for all > 1 and i ≤ n s . By induction, this implies that for all ≤ L and i ≤ n s . Moreover, for some C > 0, for all i < k and ≤ L − 1.
Using the fact that a control α ∈ A ad is predictable, one can then proceed as in the proof of Corollary 2.3 to derive a more tractable upper-bound. It appeals to the following version of Assumption 1.

A generic progressive learning algorithm
Let us now describe how the result of Corollary 3.4 can be turned into a (stochastic) dynamic programming algorithm, in the spirit of Section 2.4, that can be implemented in practice.
By Jensen's inequality, the upper-bound of Corollary 3.4 can be rewritten as where F ad to which we can associate the optimal control problem 6 It admits a dynamic programming principle that involves a Bayesian update of the prior law onθ at each step of the algorithm, see e.g. [10]. Let us first observe that, from step on, our bound only involves the components ofθ associated to the indexes in I α . We therefore set This means that the update of the prior can be restricted to a reduced number of components ofθ. This explains why we will concentrate on minimizing this upper-bound rather than directly the left-hand side of (22), which would lead to a very high-dimensional optimal control problem, at each step . This way, we expect to reduce very significantly the computation cost of the corresponding "optimal" strategy.
In order to make the updating rule explicit, we use the following assumption.
Under this assumption, we can compute the law ν α, . It is given by 6 Only the conditional law given F α of the components ofθ corresponding to indexes in I α play a role in the definition ofF ad p ( , α, ν) and F ad p ( , α, ν) below. To avoid introducing new complex notations, we shall indifferently take ν or only the conditional law of the corresponding components as an argument, depending on the context. 7 As for measurability, we identify I α −1 to the element of R ns with i-th entry given by for a Borel set A of D α −1 := T α −1 (R ns × S ns ). From this, one can deduce the law ν α ofθ α = T α (θ) given F α , in the form ν α = U( , α, ν α −1 ), by simply integrating on the components corresponding to indexes that are not in I α (meaning that U is explicit in terms of U 1 2 ).
We are now in position to state our dynamic programming principle, see e.g. [10]. Again, note that the law of f ad p ( +1, α ,θ) given F α depends onθ only throughθ α . For ease of notations, we identify all measures to an element of M (even if it supported by a space smaller than R ns × S ns ).
Proposition 3.5. Let Assumption 3 hold. Then, for all α ∈ A ad , 0 ≤ ≤ L − 2 and ν ∈ M, In principle, this dynamic programming algorithm allows one to estimate numerically the optimal policy α in a feed-back form, off-line. Importantly, solving this problem given an initial prior ν 0 is very different from first estimating the parameterθ and then solving the control problem as ifθ was given. In the first case, we take into account the risk due to the uncertainly on the true value ofθ, not in the second one.
Remark 3.6. In practice, the algorithm requires estimating and manipulating the law of a high-dimensional parameter, at least at the first steps. But the above can be modified by changing the filtration ) with τ α := inf{l ≤ L : q α l ≤ ρ} for some ρ > 0. In this case, no additional information is considered up to step τ α , the update of the prior only takes place from step τ α on and it only concernsθ α τ α whose dimension is controlled by ρ. As for the first steps of the algorithm, namely before τ ρ , one can replace f ad p by a robust version in the spirit of Section 2.4.
Remark 3.7. The algorithm also requires knowing the conditional density g α . Although, P |s can be simulated, its conditional density is not known in general. However, one can use a proxy and/or again modify the flow of information to reduce to a more explicit situation. Let us consider the situation in which (F α ) ≤L is replaced by (F α ) ≤L withF α =F α −1 ∨ σ(δμ α,i , i ∈ I α ) andF α 0 is trivial. Then, conditionally toF α −1 ∨σ(θ α −1 ), δN α (Σ α ) −1 (δμ α −μ α ) is asymptotically Gaussian as δN α increases to infinity. In practice, we can do as if δN α (Σ α ) −1 (δμ α −μ α ) was actually following a standard Gaussian distribution, conditionally toθ α and F α −1 , which provides an explicit formula for the conditional densitȳ g α of δμ α givenθ α −1 and F α −1 , to be plugged into (23). Namely, the updating procedure takes the form whereǓ is explicit. Then, if the initial prior ν 0 is such that (μ,Σ) is a Normal-inverse-Wishart distribution, all the posterior distribution ν α , ≤ L, are such that (μ,Σ) remains in the class of Normal-inverse-Wishart distributions with parameters that can computed explicitly from our simulations. Namely, if, givenF α ,Σ has the distribution 8 W −1 i α (Σ α ) andμ has the distribution N (m α ,Σ/k α ) givenΣ, then the coefficients corresponding to the law givenF α +1 are see e.g. [18,Section 9]. Later on, we shall write the corresponding law as N W −1 (p α +1 ) with

Example of numerical implementation using neural networks
In this section, we aim at solving the version of the dynamic programming equation of Proposition 3.5, using an initial Normal-inverse-Wishart prior and the approximate updating procedure suggested in Remark 3.7:F ad p ( , α, ν) = ess inf withǓ as in Remark 3.7 anď It would be tempting to use a standard grid-based approximation. However, to turn this problem in a Markovian one, one needs to let the value function at step depend on q α , N α , C α ,μ α and p α , where C α is the running cost of strategy α up to level , defined for = 0 by C α = −1 l=0 q α l δN α l+1 and C α 0 = 0. The dimension is then 1 + 1 + 1 + q α + (1 + q α + 1 + (q α ) 2 ). Even for q α = 20, the corresponding space is already much too big to construct a reasonable grid on it. We therefore suggest using a neural network approximation. Let us consider a family of bounded continuous functions {φ x , x ∈ X}, X being a compact subset of R dX for some d X ≥ 1, such that, for all q, δq ≤ n s and N, δN ≥ 1, We then fix a family {α k } k≤k of deterministic paths of A(0) and simulate independent copies {θ j } j≤j ofθ according to ν 0 , a Normal-inverse-Wishart distribution N W −1 (p 0 ). For each j, we consider an i.i.d. sequence (P j,1 j , . . . , P j,ns j ) j ≥1 in the law N (μ j ,Σ j ) withθ j =: (μ j ,Σ j ). We take these sequences independent and independent ofθ. For each k and j, we denote by (μ k,j ) ≤L , (p k,j ) ≤L and (I k,j ) ≤L the paths (μ α k ) ≤L , (p α k ) ≤L and (I α k ) ≤L associated to the j-th sequence (P j,1 j , . . . , P j,ns j ) j ≥1 and the control α k . Similarly, we write f ad,k,j p ( , ·) to denote the function f ad p ( , ·) defined as in (21) but in terms of I k,j −1 in place of I α −1 . Given an integer r ≥ 1, we first computex L−1 as the argmin over x ∈ X of k k=1j j=1 L−1 means that the expectation is taken only overμ according to the law ν k,j L−1 , i.e. N W −1 (p k,j L−1 ), and (·) i means that we take the i-th component of the vector in the brackets. Then, for any α ∈ A ad , we setφ where E ν k,j means that the expectation is computed over (μ α k +1 , p α k +1 ,θ,m α k ) given (μ k , p k ) = (μ k,j , p k,j ) and using the prior ν k,j onθ associated to p k,j . Then, we seť φ (q α , N α , C α , ·) := min δN ≤ K}, and so on until obtaining φ 0 (n s , 0, 0, 0, p 0 ). By continuity of φ · (·) and compactness of X and A( , α) for α given, the minimum is achieved in the above, possibly not unique, and one can choose a measurable map a such that a (q α , N α , C α , ·) ∈ arg min for all α ∈ A ad . Then, given the parameter p 0 of our initial prior ν 0 , our estimator of the optimal policy is given by α = (q , N ) defined by induction by (δq 1 , δN 1 ) = a 0 (n s , 0, 0, 0, p 0 ) and (δq +1 , δN +1 ) = a (q , N , C ,μ α , p α ) for 0 < < L.
Note that the above algorithm for the estimation of the optimal control only requires off-line simulations according to the initial prior ν 0 . It is certainly costly but does not require to evaluate the real financial book, it can be trained on a proxy, and can be done off-line. It can be combined with the approach of Remark 3.6 to reduce the computation time. In order to prepare for the use of a different initial prior, one can also slightly adapt the above algorithm by considering different initial values of p 0 (e.g. drawn from another distribution around p 0 ), so as to estimateφ 0 not only at the point p 0 . When applied to the real book, the update of the prior according to (24) leads to an additional cost that is negligible with respect to the simulation of the book. It leads to the computation of new priors associated to the financial book at hand, that can be used for a new estimation of the optimal policy or simply as a new initial prior for the next computation of the ES. An example of a simple practical implementation is detailed in Appendix B, while numerical tests are performed in Section 4.

Numerical Experiments
This section is dedicated to first numerical tests of the different algorithms presented in the previous sections. The settings of the experiments are as follows. We first choose a Normal-inverse-Wishart prior distribution ν 0 with parameters p 0 := (m 0 , k 0 , i 0 , Σ 0 ). The vector m 0 is represented on Figure 6 with m i 0 = µ i , i ≤ n s , and Σ 0 = (i 0 − n s − 1)Σ where Σ has entries with ρ = 0.6 or ρ = 0 depending on the experiments below. As for k 0 and i 0 , they are chosen equal to 300, meaning that we have a low confidence in our prior. The computing power is K = 10 7 .
We apply the four different algorithms on 5 000 runs (i.e. 5 000 independent implementations of each algorithm). For each run, we • first simulate a value for the real scenarios and covariance matrices μ,Σ ∼ N W −1 (p 0 ), • apply each of the four algorithms, with simulated prices following P |s ∼ N μ,Σ , • for each algorithm, we measure the relative error The four algorithms that we compare are: • A Uniform Pricing Algorithm: All the scenarios are priced with K/n s Monte Carlo simulations, and the estimator ES is the average of the n w = 6 worst scenarios. This is the most naive method, with only one step and where all scenarios are priced with an equal number of Monte Carlo simulations. It serves as a benchmark.
• The Heuristic Algorithm: We use the 2-levels strategy described in Section 2.5 with the book sample parameters of Table 1 and the computation parameters of Table 2. We do not evaluate the constant c of Assumption 1 but simply set it to 0, see Remark 2.4. The optimal strategy is given by (q 0 , q 1 , N 1 , N 2 ) = (253, 68, 17 297, 100 000).
• The Deterministic Algorithm: We run the deterministic algorithm of Section 2.4 optimized with µ = m 0 as the values of the scenarios, Σ with ρ = 0.6 as the covariance matrix and L = 4. Note that using the real mean parameter as an entry for optimization is quite favorable for this algorithm, although the "true" parameter of each run will actually deviate from this mean value. This gives us the strategy (q 0 , q 1 , q 2 , q 3 , N 0 , N 1 , N 2 , N 3 , N 4 ) = (253, 35, 10, 6, 0, 6 000, 44 000, 44 000, 1 235 666), which we apply to each run.
• The Adaptative Algorithm: We do the training part of the adaptative algorithm using our prior p 0 := (m 0 , k 0 , i 0 , Σ 0 ), with ρ = 0.6, as parameters and L = 4. We use a very simple one hidden-layer neural network. It could certainly be improved by using a more sophisticated multi-layers neural network, but this version will be enough for our discussion. Details on the implementation are given in the Appendix B. Once this is done, we apply the optimal adaptative strategy on each run.

Positively correlated scenarios ρ = 0.6
In this first experiment, the simulated runs use the values ρ = 0.6 and i 0 = k 0 = 300.
To get an idea of how much noise is added to the average scenario values in our simulations, we plot in Figure 7 the prior value m i 0 for each scenario of index i ≤ n s (this is the line) and the first 20μ i j out of the 5 000 runs for each scenario (these are the points). Compared to the deterministic algorithm, we see that the adaptative one uses much less Monte Carlo simulations at the final steps and focuses more on the intermediate steps to select the worst scenarios. The deterministic algorithm is also more aggressive in the choice of q 1 and q 2 . This can be easily explained by the fact that the latter believes that the real distribution is not far from the solid curve on Figure 7 (up to standard deviation) while the adaptative one only knows a much more diffuse distribution corresponding to the cloud of points of Figure 7 since his level of uncertainty is quite high for our choice i 0 = k 0 = 300.
On Figures 8-11, we plot the histograms of the relative errors. We see that the distribution is tightest for the deterministic algorithm, followed quite closely by the adaptative algorithm. Both of them perform very well. As expected, the uniform algorithm is very poor. Note that the heuristic one already very significantly improves the uniform algorithm, although it does not reach the precision of the two most sophisticated algorithms (without surprise). Because of the huge uncertainty mentioned above, the adaptative algorithm is rather conservative while the deterministic algorithm makes full profit of essentially knowing the correct distribution, and performs better. We will see in our second experiment that things will change when we will deviate from the parameters used for optimizing the deterministic algorithm (by simply passing from ρ = 0.6 to ρ = 0 in the simulated runs). In Table 3, we provide the L 1 and relative errors (with standard deviations), the L 2 error and the number of correct selections, that is the number of runs for which a given algorithm has chosen the correct worst 6 scenarios. In terms of L 1 or L 2 error, the relative performance of the algorithms is as above. However, if we look at the number of correct selections, we see that the adaptive algorithm performs better than the other 3 algorithms. Again, by comparing the strategies of the deterministic and the adaptive algorithms, we see that those of the adaptative algorithm are more conservative on the ranking and filtering part versus the final pricing as it puts relatively more Monte Carlo simulations to detect the correct scenarios and relatively less for their estimation.   In Figure 13, we provide, for the first 4 runs, the values and real ranks of the 6 worst scenarios selected by each algorithm. The numbers displayed are the true ranks of the selected scenarios given byμ and their y-coordinate is the value obtained when running the algorithm. "Real" is the real values as sampled.

Uncorrelated scenarios ρ = 0
We now do the numerical test with ρ = 0 as the true correlation. The deterministic and adaptative algorithm are still trained with ρ = 0.6, but P |s is simulated using ρ = 0. On Figures 14-17, we show the histograms of the relative errors. We see that the distribution of the relative errors is now tightest for the adaptative method, followed by the deterministic method, then by the heuristic and the uniform methods. Furthermore, we see that the distribution corresponding to the deterministic method is significantly biased to the left. This is actually true for all algorithms, but at a less significant level. This suggests that we now have a large part of the error that does not come from the final pricing error, but from errors in the selection of scenarios. In Table 4, we provide the L 1 and relative errors (with standard deviations), the L 2 error and the number of correct selections for the 4 algorithms. For all algorithms, compared to the case ρ = 0.6, we see that we have simultaneously a lower number of correct selections of scenarios (which we could expect to increase the errors) and a lower L 1 error. This surprising result is explained by the fact that lowering the correlation has two effects. The filtering and ranking part of the algorithm becomes harder, as can be seen from Corollary 2.3. This explains why the number of correct selections becomes lower. However, we compute at the end an average over the n w worst scenarios and the error on this average is lower when the pricings are uncorrelated compared to the case where they exhibit a positive correlation.
The adaptative algorithm has now simultaneously the lowest L 1 and L 2 errors, as well as the highest number of correct selections. We see that it is especially good in L 2 error, so we expect it to present a very low number of large errors. As, by construction, it has been trained to detect misspecifications of the parameters, it now has a clear advantage on the deterministic algorithm which does not see it. This results in an improvement of almost 20% of the L 2 error.
Following the above reasoning, we understand that, compared to the previous experiment, the final pricing error now plays a smaller role and the ranking and selection error a bigger role, which explains why the histogram of the errors for the determinist algorithm is strongly biased to the left, as it now incorrectly selects scenarios more often.  where X is the absolute error of the algorithm on a run. As was suggested by the L 2 errors of Table 4, we see that the tail distribution of errors is lowest for the adaptative algorithm, followed by the deterministic algorithm (for big errors), and then by the heuristic and uniform algorithms.

Conclusion
We propose in this paper different algorithms for the computation of the expected shortfall based on given historical scenarios. All are multi-steps algorithms that use Monte Carlo simulations to reduce the number of historical scenarios that potentially belong to the set of worst scenarios. We provide explicit error bounds and we test them on simulated data deviating from the true values of the historical impacts used for computing the associated optimal strategies. The first algorithm is a very easy to implement 2-steps algorithm that already provides relatively small errors on our numerical tests. A four step deterministic dynamic programming algorithm performs very well when real datas are not far from the parameters used in the optimization procedure. It seems even to be quite robust, as shown by our numerical test in the case where the true correlation parameter is not the one used for computing the optimal policy. Finally, we propose an adaptative algorithm that aims at learning the true value of the parameters at the different steps of the algorithm. Our first numerical tests suggest that it is more conservative than the deterministic one, but probably more robust to parameters misspecifications, as expected. The version we use is built on a very simple one hidden layer neural network and can certainly be considerably improved for industrial purposes.
A Proxy of the optimal strategy for the heuristic (19) In the case p = 1, (19) can even be further simplified by using the upper-bound The right-hand side of (26) is now tractable for minimization. Given, the optimal policy q h 1 is defined by the following  For simplicity, let us consider the case c = 0, see Remark 2.4. On Figure 19, the square is q 1,2, * 1 = 52.41, the circle is q 2, * 1 = 68.33 and the cross is the real optimum q * 1 = 71 of h 1 0 , for the parameters of Tables 1 and 2. We see that we actually almost reach the correct minimum. It corresponds (up to rounding) to N 1, * 1 = 23723, N 2, * 1 = 17148, N * 1 = 15934. 9 We optimize here over real positive numbers. ). Circle: h 1 0 (q 2, * 1 ). Cross: h 1 0 (q * 1 ).
Using the same set of parameters, we plot on Figure 20 the two functions h 1 0 andh 1 . Although these two functions have values of different orders of magnitude, their shapes are quite close, which explains why we manage to obtain a relatively good approximation for the minimizer.

B Precise implementation of the neural network algorithm
In this Appendix, we describe in more details how the neural network approximation of the optimal policy of the adaptative algorithm is constructed. All the parameters values are given in Tables 6, 7 and 8 below.

B.1 Initialization
• In practice, the neural network input's size depends on the window size q. Therefore, we need to train different neural networks for each window size. In order to get enough points to train each of these neural networks, we have chosen the grid q g = [6,10,15,20,25,30,35,40,45,50,60,70,80,90,100,150,200,253] of possible values for q.
• We simulate independent copies {θ j } j≤j bar = {(μ j ,Σ j )} j≤j bar ofθ, where j bar is given in Table  8. For each 1 ≤ j ≤ j bar,Σ j is an inverse-Wishart of parameters i 0 , Σ 0 , andμ j is a Gaussian random vector of mean m 0 and covariance matrixΣ j /k 0 . The parameters i 0 , k 0 and Σ 0 are defined in Table 8 and (25), while m i 0 = µ i , i ≤ n s , with the µ i 's of Figure 6.

B.2 Strategy Generation
To generate the deterministic strategies (α k ) k≤k bar , where k bar is given in Table 8, we proceed as follows.
• For each 1 ≤ k ≤k bar, we simulate L + 1 uniform random variables (U n ) L n=0 between 0 and 1. We sort them in increasing order U s(n) L n=0 and define a cost K := K(U s( ) − U s( −1) ) when 1 ≤ ≤ L − 1, and K L = K(U s(0) + 1 − U s(L−1) ). The idea is that we select L + 1 points randomly on a cercle of total length K: we choose one of these points, and starting from it, the computational power that we will use at each level 1 ≤ ≤ L − 1 is the length of the arc between the previous and the next point. For the last step, we take K times the length between the points L − 1 and 0, so as to put, in average, twice more computational power on this last step.
• Once we have the computational cost for each step, we can choose the q for each strategy, so that we can deduce δN +1 := K /q . For = 0, we choose q index 0 = 18, where 18 is the number of terms in the grid q g , which therefore gives q 0 = q g [q index 0 ] = n s . For = L − 1, we choose q index L−1 = 0, that is, q L−1 = n w . For 1 ≤ ≤ L − 2, we choose q index as a random integer between [L − , q index −1 − 1]. The choice of q is then q = q g [q index ]. We check that the sequence (N ) 1≤ ≤L is non-decreasing. If this is the case, we keep it, if not, we reject it and do another run.

B.3 Forward Pass
The next step is to generate all prices and execute for each k and each j the strategy k.
• Updating Σ k,j from level − 1 to level can use a lot of memory. Indeed, consists in δN α +1 × |q +1 | 2 terms, which can quickly exceed memory limits. Therefore, we do the sum with only N memory new pricings opt terms at a time, see Table  8 below.

B.4 Computation of f precompute, running costs and admissible sets
• In order to speed up the computation time, we now precompute several values that will be used many times afterwards. First, we compute f precompute( , k, j) as f ad 1 ( , ·) at the point corresponding to (k, j) except that, in the definition of f ad 1 ( , ·), we replace the random permutationm by its estimation from the previous step,μ by its average under the posterior distribution at + 1, andσ by its estimation at step + 1.
• We compute running cost( , k) := C α k of each k at step .
• We restrict the set of possible actions at step , given that we have followed the strategy k so far, to admissible sets( , k) defined as the collection of {(δq k +1 , δN k +1 ), k ≤ k bar}, such that q k + δq k +1 ∈ q g , N k +1 > N k , running cost( , k) + q k δN k +1 ≤ max 1≤k ≤k running cost( + 1, k ).
The last condition avoids inducing a strategy with a running cost that is not present in our data set, when doing the one step optimization.

B.5 Computation of the final expectations
We first pre-compute the quantities by Monte Carlo using N e simulations. As the simulation of an inverse-Wishart random variable is significantly slower than the simulation of a Gaussian random variable, we only simulate 1 inverse-Wishart for N p Gaussians. The values of N e and N p are given by N mu tildes simulated and N wishart proportion of Table 8. The estimation is called expectation k,j L .
B.6 Training of the neural network at level L • We use a neural network with one inner layer with 256 neurons and 1 output layer with 1 neuron to fit (expectation k,j L ) j≤j bar,k≤k bar . The neurons of the inner layer consist of the composition of the softplus function with an affine transformation of the inputs. For each batch, the algorithm used for the training is the standard gradient descent of Tensorflow. We do N Iter training steps in total. The learning rate used is given in Table 7. In order to bring the input values of the parameters close to 0 and 1, we renormalize them according to the values in Table 6.
B.7 Computation of the expectations at level L − 1 We now estimate whereφ L is the fit of (expectation k,j L ) j≤j bar,k≤k bar from the previous step. The most cpu demanding part is no more the simulation of the inverse-Wisharts, but the updates of the parameters of the inverse-Wishart. Therefore, we simulate as many Gaussian random variables as inverse-Wishart random variables, with N e given by N mu tildes simulated non final level of Table 8. For our computations, we need to update Σ k,j L−1 to the corresponding posterior parameter according to (24). This can however lead to an enormous amount of multiplications and additions. Therefore, instead of updating the whole matrix, we only update the diagonal terms according to (24) and estimate non diagonal terms by keeping the correlation terms equal to the ones of Σ k,j L−1 . This enables us to approximately gain a factor of q k L in speed in this critical step.
B.8 Training of the neural network at level L − 1 • To fit the expectation of the previous step, we use a neural network with the same structure as in level L, with the same cost function.
• The initialization, choice of batches, and training of the neural network are the same as for the level L. The number of iteration, learning rate, and renormalization constants are given in Tables 6, 8 and 7.
• We take j batch = min (j batch size, j bar) and k batch = min (k batch size, k bar), where j batch size and k batch size are defined in Table 8.
B.9 Computation of the expectations at levels 0 ≤ ≤ L − 2 • The expectations at step are computed by Monte Carlo after replacing the value function at step + 1 by its neural network approximation, and f ad 1 ( , ·) by f precompute( , ·). • We simulate as many Gaussian random variables as inverse-Wishart random variables, with N e given by N mu tildes simulated non final level of Table 8.
• We not not fully update Σ k,j to the corresponding posterior parameter but proceed as in level L−1.
B.10 Training of neural networks at levels 0 ≤ ≤ L − 2 • We now have to optimize over q k ∈ q g . Therefore, we must now train up to |q g | different neural networks (with different inputs' sizes). In practice, we only train neural networks indexed by q ∈ (q k ) 1≤k≤k bar ⊂ q g , that is, for all the choices of q that are obtained by at least one strategy at level .
• We must also choose a δN that should be added as an entry of the neural network before optimizing. Furthermore, to help the neural networks converge, we decided to add f precompute( , j, k) as an input.
• The loss function and the structure of the neural network is as above, and we still use Xavier initialization, and bring the inputs of the neural networks to reasonable values close to 0 and 1 by renormalizing them using the constants of Table 6.
• Compared to levels L and L − 1, the choice of batches is slightly different. Indeed, to train a neural network associated to q ∈ q g , we only use strategies such that q k = q. To do so, we first define S q = {k ∈ [[1, k bar]] : q k = q}. We then define k batch = min (k batch size, |S q |) and j batch = min (j batch size, j bar). We then proceed nearly identically as for levels L and L − 1.
We select a new batch every N batch change proportion, composed of indices 1 ≤ m a (j) ≤ j batch, 1 ≤ m b (k) ≤ k batch, where m a and m b are uniform random permutations of [[1, j bar]] and S q . For each batch, the algorithm used for the training is again the standard gradient descent of Tensorflow.
• Compared to levels L and L − 1, we found that making the neural networks converge was much harder. In particular, the learning rate had to be really fine tuned. In order to automatize the process, for each q, we proceed as follows. We do not intanciate one, but number of neural networks for learning rate test neural networks. For each of these neural networks, we do N Iter learning rate test training steps, but use different learning rates for each. For the first neural network, we use base learning rate as the learning rate, for the second, base learning rate/10, and for the k-th, base learning rate/10 k−1 . For each of these neural networks, we store at each iteration step the log error. Once the N Iter learning rate test training steps have been done for each of these neural networks, we keep the neural network instance that has the lowest average log error. If it is the k-th neural network, we then train it again for N Iter training steps, using as learning rate base learning rate/10 k .

B.11 Parallelization
In practice, we parallelize the forward pass according to the strategy indices k. We run thread batch size processes in parallel, where thread batch size is defined in Table 8.
At a given level , the computation of expectation k,j can be parallelized according to the sample indices j. In practice, we run number of threads for level expectations number of processes in parallel, where number of threads for level expectations is defined in Table 8.
For a given level, the training of each neural network corresponding to a given q ∈ q g can be done independently. Therefore, at a given level, we multiprocessed our code in order to train all the neural networks in parallel.