Sarah Kaakaï

Anis Matoussi

Achraf Tamtalini

Estimation of Systemic Shortfall Risk Measure using Stochastic Algorithms *

Keywords: Multivariate risk measures, shortfall risk, stochastic algorithms, stochastic root finding, risk allocations

Systemic risk measures were introduced to capture the global risk and the corresponding contagion effects that are generated by an interconnected system of financial institutions. To this purpose, two approaches were suggested. In the first one, systemic risk measures can be interpreted as the minimal amount of cash needed to secure a system after aggregating individual risks. In the second approach, systemic risk measures can be interpreted as the minimal amount of cash that secures a system by allocating capital to each single institution before aggregating individual risks. Although the theory behind these risk measures has been well investigated by several authors, the numerical part has been neglected so far. In this paper, we use stochastic algorithms schemes in estimating Multivariate Shortfall Risk Measure (MSRM) and prove that the resulting estimators are consistent and asymptotically normal. We also test numerically the performance of these algorithms on several examples.

Introduction

The axiomatic theory of risk measures, first initiated by the seminal paper of [START_REF] Artzner | Coherent measures of risk[END_REF], has been widely studied over the past decades. The Value-at-Risk (VaR) is one of the best known and most common risk measures used by practitioners and regulation authorities. However, the VaR lacks one important property: it does not take into account the diversification effect. To circumvent this problem, the VaR was replaced by the Conditional Value-at-Risk (CVaR) and a more general framework of improved risk measures has been introduced: Utility-based Shortfall Risk (SR). Nevertheless, when it comes to a system of financial institutions or portfolios, the question about how to assess the global risk, as well as individual risks, arises. Following the 2008 crisis, the traditional approach of measuring systemic risk that consists in considering each institution as a single entity isolated from other institutions, has shown some limits. Indeed, with this approach, the risk associated with a vector of positions X = (X 1 , ..., X d ) can be written as:

R(X) := d i=1 η i (X i ),
where each η i is a univariate risk measure. Then, C. [START_REF] Chen | An axiomatic approach to systemic risk[END_REF] proposed an approach whose spirit is very close to the axiomatic framework initiated by [START_REF] Artzner | Coherent measures of risk[END_REF]. They showed that any systemic risk measure verifying certain axioms can be written as the composition of a univariate risk measure η and an aggregation Λ, i.e., R(X) = η(Λ(X)).

The previous representation is known as the "Aggregate then Add Cash" approach as it consists first in aggregating the positions X 1 , ..., X d through the aggregation function Λ and then apply a univariate risk measure. One of the most common ways to aggregate the outcomes X i is to simply take the sum, that is to consider Λ(x) = d i=1 x i . It is worth noticing that, while summing up profit and losses might seem reasonable from the point of view of a portfolio manager because portfolios' profits and losses compensate each other, this aggregation rule seems inadequate from the point of view of a regulator where cross-subsidization between institutions is not realistic, since no institution will be willing to pay for the losses of another one.

Motivated by these considerations, [START_REF] Biagini | A unified approach to systemic risk measures via acceptance sets[END_REF] proposed another approach to measure the systemic risk. They first considered the systemic risk as the minimal capital that secures the system by injecting capital into the single institutions, where A is an acceptance set. This approach, known as "Add Cash then Aggregate" consists in adding the amount m i to the financial position X i before the corresponding total loss Λ(X + m) is computed. The systemic risk is then measured as the minimal total amount d i=1 m i injected into the institutions to make it acceptable. With this approach, a joint measure of total risk as well as individual risk contributions to systemic risk is obtained.

If m * = (m * 1 , ..., m * d ) is an optimum, that is R(X) = d i=1 m * i and Λ(X + m * ) ∈ A, one could rank the m * i 's and hence be able to say that institution/portfolio i requires more cash allocation or is riskier than j if m * i ≥ m * j . In this article, we are interested in the numerical computation of the multivariate shortfall risk measure (MSRM) that was introduced in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF]. They are an extension of univariate SR measures, and can be obtained by taking the aggregation function Λ(x) = l(-x) where l is a multivariate loss function (see Section 1) and the acceptance set A = {X ∈ L 0 (R), E[X] ≤ 0}.

In [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF], MSRMs and their optimal allocations are approximated using a two step procedure. The expectation in the aggregation function is first approximated using Monte Carlo, Fourrier or Chebychev schemes. The minimum in equation (0.1) is then found by a built-in minimizing algorithm in python. However, the minimizing algorithm can be sensitive to the starting point of the algorithm, and there are no theoretical results on the convergence of this procedure. In the univariate case, SR can also be characterized as the unique root of a function g : R → R that is expressed as an expectation, which can be estimated by using a deterministic root finding algorithm combined with a Monte Carlo procedure that estimates g(s) for all s ∈ R.

Another idea to compute the VaR and CVaR comes from the fact that they are solutions and the value of the same convex optimization problem as pointed out in [START_REF] Rockafellar | Optimization of conditional value-at-risk[END_REF]. Moreover, as they can be expressed as an expectation, this led [START_REF] Bardou | Computing VaR and CVaR using stochastic approximation and adaptive unconstrained importance sampling[END_REF] to define consistent and asymptotically normal estimators of both quantities using a classical Robbins-Monro (RM) procedure. For SR, a constrained RM algorithm has also been proposed by [START_REF] Dunkel | Stochastic root finding and efficient estimation of convex risk measures[END_REF].

In this paper, we introduce stochastic algorithms for MSRM, and prove their consistency and asymptotic normality. Indeed, the optimal allocations of MSRM can also be characterized as the root of a function that is expressed as an expectation. However, the multivariate setting adds additional constraints: the optimal allocations are characterized as the solution of a first order condition of the Lagrangian associated with the MSRM. Hence, the new problem is to find a saddle point , z)] determined by the first order conditions. RM stochastic algorithms have been the subject of a vast literature (see e.g. [START_REF] Benveniste | Adaptive algorithms and stochastic approximations[END_REF], [START_REF] Duflo | Algorithmes stochastiques[END_REF]). The basic paradigm in its simplest form is to approximate z * by the following stochastic difference equation:

z * = (m * , λ * ) ∈ R d × R + , which is the unique root of a function h(z) = E[H(X
Z n+1 = Z n + γ n Y n , where
γ n > 0 is the step size that goes to zero as n → ∞, and Y n = H(X n+1 , Z n ), where (X n ) is a sequence of i.i.d random variables with the same law as X. Moreover, if the random variable X is not directly simulatable, [START_REF] Frikha | Multi-level stochastic approximation algorithms[END_REF] extended the scope of multi-level Monte Carlo to the framework of stochastic algorithms. In many cases, the analysis of these algorithms is based on the so-called Ordinary Differential Equation (ODE) method introduced by [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF]. The main idea is to show that, in the long run, the noise is eliminated so that, asymptotically, the behaviour of the algorithm is determined by that of the average ODE:

ż = h(z).
To ensure the convergence of RM algorithms, the function h is assumed to have sub-linear growth. This assumption is often not verified for MSRM, due to the choice of loss functions. In order to avoid reducing drastically the scope of the applications, we introduce a constrained version of the RM algorithm for MSRM. This consists in using the projection into a compact K each time the sequence Z n goes out of K. First introduced by [START_REF] Kushner | Stochastic approximation for constrained systems with observation noise on the system and constraint[END_REF], an excellent survey on projection techniques can be found in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]. In the case of MSRM, the convergence of the constrained algorithm toward the optimal allocation z * is then closely related to the study a projected multivariate ODE, whose solutions remain on K. We provide a rigorous analysis of the asymptotic stability of this ODE, based on the introduction of a Lyapunov function and on the extension of the Lasalle invariant principle to this ODE with discontinuous r.h.s. Under additional assumptions, we also establish the asymptotic normality of the estimator.

The analysis and estimation of the asymptotic covariance matrix is also more difficult due to the multivariate setting. By introducing an averaged "Polyak-Ruppert" version of the stochastic algorithm, we obtain a simpler expression of the asymptotic covariance matrix and propose an estimator allowing the computation of empirical confidence intervals.

Stochastic algorithms also offer a flexible framework that can be conveniently combined with variance reduction techniques such as Importance Sampling (IS). [START_REF] Dunkel | Stochastic root finding and efficient estimation of convex risk measures[END_REF] have adapted the RM for unidimensional SR in order to include IS. However, in order to apply this non parametric approach, it is necessary to find beforehand an appropriate change of measure for each parameter z, which can be very difficult even in the one dimensional case. Inspired by the parametric approach of [START_REF] Bardou | Computing VaR and CVaR using stochastic approximation and adaptive unconstrained importance sampling[END_REF] for the VaR/CVaR, we also propose an adaptive RM algorithm based on recursive IS. This also provides a general variance reduction method for univariate SR, which doesn't require prior knowledge of the change of measure.

The paper is organized as follows. Section 1, is dedicated to MSRM. In particular, we recall the main theorem characterizing the optimal allocations. The main results of consistency and asymptotic normality of the RM algorithm are presented in Section 2, along with the Polyak-Ruppert averaging and adaptative importance sampling. Finally, section 3 is devoted to some numerical experiments of our procedures. We present a first testing example with an exponential loss function, where we have a closed formula for optimal allocations. We also give a second example using a loss function with a mixture of positive part and quadratic functions.

About Multivariate Risk Measures

Let (Ω, F, P ) be a probability space, and denote by L 0 (R d ) the space of F-measurable dvariate random variables on this space with d ≥ 2. For x, y ∈ R d , we say that that x ≥ y We say, for example, for

( x > y resp.) if x k ≥ y k (x k > y k resp.) for every 1 ≤ k ≤ d.
X, Y ∈ L 0 (R d ), that X ≥ Y (or X > Y resp.) if P (X ≥ Y ) = 1 (or P (X > Y ) = 1 resp.).
To simplify the notation, we will simply write L 0 instead of L 0 (R d ).

Now, let X = (X 1 , ..., X d ) ∈ L 0 be a random vector of financial losses, i.e., negative values of X k represents actually profits. We want to assess the systemic risk of the whole system and to determine a monetary risk measure, which will be denoted R(X), as well as a risk

→ (-∞, ∞] is called a loss function if: (A1) l is increasing, that is l(x) ≥ l(y) if x ≥ y;
(A2) l is convex and lower-semicontinuous with inf l < 0;

(A3) l(x) ≥ x k -c for some constant c.

Furthermore, a loss function l is said to be permutation invariant if l(x) = l(π(x)) for every permutation π of its components.

Comment:

The property (A1) expresses the normative fact about the risk, that is, the more losses we have, the riskier is our system. As for (A2), it expresses the desired property of diversification. Finally, (A3) says that the loss function put more weight on high losses than a risk neutral evaluation.

Example 1. Let h : R → R be one dimensional loss function satisfying condition (A1), (A2) and (A3). We could build a multivariate loss function using this one dimensional loss function in the following way:

(C1) l(x) = h( x k ); (C2) l(x) = h(x k ); (C3) l(x) = αh( x k ) + (1 -α) h(x k ) for 0 ≤ α ≤ 1.
More specifically, in (C1), we are aggregating losses before evaluating the risk, whereas in (C2), we evaluate individual risks before aggregating. The loss function in (C3) is a convex 1. About Multivariate Risk Measures combination of those in (C1) and (C2).

One of the main examples we will be studying in this paper are the two following ones:

h(x) = 1 α + 1 ( i e βx i + αe β x i ) - α + d α + 1 , h(x) = i x i + 1 2 i (x + i ) 2 + α i<j x + i x + j ,
where the coefficient α > 0 is called the systemic weight and β > 0 is a risk aversion coefficient.

In the following, we will consider multivariate risk measures defined on Orlicz spaces (see [START_REF] Rao | Theory of orlicz spaces[END_REF] for further details on the theory of Orlicz spaces). This has several advantages. From a mathematical point of view, it is a more general setting than L ∞ , and in the same time, it simplifies the analysis especially for utility maximization problems.

Therefore, we will consider loss vectors in the following multivariate Orlicz heart: [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF] for more details about Orlicz spaces.

M θ = {X ∈ L 0 : E[θ(λX)] < ∞, ∀λ > 0}, where θ(x) = l(|x|), x ∈ R d . See Appendix in
Next, we give the definition of multivariate shortfall risk measures as it was introduced in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF].

Definition 1.2. Let l be a multivariate loss function and X ∈ M θ , we define the acceptance set A(X) by:

A(X) := {m ∈ R d : E[l(X -m)] ≤ 0}.
The multivariate shortfall risk of X ∈ M θ is defined as:

R(X) := inf m k : m ∈ A(X) = inf m k : E[l(X -m)] ≤ 0 . (1.1)
Remark 1.3. When d = 1, the above definition corresponds exactly to the univariate shortfall risk measure in Föllmer and Schied (2002).

The following theorem from [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF] shows that the multivariate shortfall risk measure has the desired properties and admits a dual representation as in the case of univariate shortfall risk measure. We introduce Q θ * the set of measure densities in L θ * , the dual space of M θ :

Q θ * := dQ dP := (Z 1 , ..., Z d ), Z ∈ L θ * , Z k ≥ 0 and E[Z k ] = 1 for every k .
Theorem 1.4. [Theorem 2.10 in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF]] The function

R(X) := inf m k : m ∈ A(X) ,
is real-valued, convex, monotone and translation invariant. Moreover, it admits the dual 1. About Multivariate Risk Measures representation:

R(X) = max Q∈Q θ * {E Q [X] -α(Q)}, X ∈ M θ ,
where the penalty function is given by

α(Q) = inf λ>0 E λl * dQ λdP , Q ∈ Q θ * .
Now, we address the question of existence and uniqueness of a risk allocation which are not straightforward in the multivariate case. [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF] showed that if the loss function is permutation invariant, then risk allocations exist and they are characterized by Kuhn-Tucker conditions. We denote by Z = {m ∈ R d , m i = 0} the zero-sum allocations set.

Definition 1.5. A risk allocation is an acceptable monetary risk allocation m ∈ A(X) such that R(X) = m k . When a risk allocation is uniquely determined, we denote it by RA(X).

We make the following assumption on the loss function l and the vector of losses X ∈ M θ :

(A l ) i. For every m 0 , m → l(X -m) is differentiable at m 0 a.s.;
ii. l is permutation invariant.

Theorem 1.6. [Theorem 3.4 in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF]] Let l be a loss function and X ∈ M θ such that assumption (A l ) holds. Then, risk allocations m * ∈ R d exists and they are characterized by the first order conditions:

1 = λ * E[∇l(X -m * )], E[l(X -m * )] = 0,
where λ * ≥ 0 is a Lagrange multiplier. If moreover l(x + •) is strictly convex along zero sum allocations for every x such that l(x) ≥ 0, the risk allocation is unique.

Comment: Let f 0 (m) = d i=1 m i and f 1 (m) := E[l(X -m)], for m ∈ R d and X ∈ M θ .
The assumption (A l )-i. together with the convexity of the function m → l(X -m), we have that, by Theorem 7.46 in [START_REF] Shapiro | Lectures on stochactic programming: Modeling and theory[END_REF],

f 1 is differentiable at every m ∈ R d and that, ∇f 1 (m) = -E[∇l(X -m)], m ∈ R d .
Therefore, the first order conditions given in the above theorem are equivalent to :

   ∇f 0 (m * ) + λ * ∇f 1 (m * ) = 0, f 1 (m * ) = 0.
Furthermore, we also know, thanks to Theorem 28.3 in Rockafellar (1970), that the above conditions are equivalent to saying that (m * , λ * ) is a saddle point of the Lagrangian associated

Multivariate Systemic Risk Measures and Stochastic Algorithms

to the problem in (1.1),

L(m, λ) := f 0 (m) + λf 1 (m) = d i=1 m i + λE[l(X -m)].
(1.2)

Under the assumptions of the above theorem z * := (m * , λ * ) is the unique solution of h(z) = 0, where:

h(z) :=   λE[∇l(X -m)] -1 E[l(X -m)]   , z = (m, λ) ∈ R d × [0, ∞[. (1.3)
Thus, in order to find the unique risk allocation m * , we can look for the zeros of the function h. We introduce in the section a stochastic algorithm in order to approximate z * , with the advantage of being incremental, less sensitive to dimension, and with offer a flexible framework that can be conveniently combined with features such as importance sampling, as discussed in Section 2.5.

Multivariate Systemic Risk Measures and Stochastic Algorithms

Let l be a loss function and a vector of losses X ∈ M θ . In the remainder of the paper, we assume that the hypotheses of Theorem 1.6 are verified by l, that is:

(A ′ l ) i. For every m 0 , m → l(X -m) is differentiable at m 0 a.s.; ii. l is permutation invariant; iii. m → E[l(X -m)] is strictly convex.
Under (A ′ l ), Theorem 1.6 ensures that there exists a unique risk allocation m * such that z * = (m * , λ * ) is the unique root of the function h(z) := E[H(X, z)], where we set

H(X, z) = λ∇ m l(X -m) -1 l(X -m) , X ∈ M θ . (2.1)
Note that the assumptions (A ′ l ) (i) and (iii) allow us to consider loss functions of interest, which are not necessarily differentiable and/or strictly convex, such as the second example in Section 3.2. The aim of this section is to construct a stochastic algorithm that converges to the root z * = (m * , λ * ) under some suitable assumptions. The standard Robbins-Monro (RM) procedure approximates the function h

(z) = E[H(X, z)] at the n + 1th step by Z n+1 = (m n+1 , λ n+1 ) = Z n + γ n H(X n+1 , Z n ),
with X n+1 a random variable distributed as X. In order to obtain a.s. convergence of the algorithm toward the unique root z * , one crucial condition is the sublinearity of h, which is very constraining and does not allow us to define a framework general enough to cover a wide range of loss functions. Consequently, we propose a constrained version of the RM algorithm. The constrained algorithm requires the sequence (Z n ) n≥1 to remain in a compact set K, for non-explosion reasons. Each time Z n goes out of K, it is replaced by its projection on K leading to the following projected RM algorithm:

Z n+1 = Π K [Z n + γ n H(X n+1 , Z n )], Z 0 = z 0 ∈ K = Π K [Z n + γ n h(Z n ) + γ n δM n ],
(2.2)

where δM n = H(X n+1 , Z n ) -h(Z n ), and (X n ) n≥1 is an i.i.d sequence of random variables
with the same distribution as X, independent of Z 0 . An important point is that the Algorithm (2) approximates simultaneously the optimal allocation m * and the Kuhn-Tucker coefficient λ * , which differs from the univariate case of [START_REF] Dunkel | Stochastic root finding and efficient estimation of convex risk measures[END_REF] where an algorithm which approximates m * directly can be implemented. We refer to Chapter 4 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] for a general introduction on constrained stochastic algorithms.

In the sequel, we denote

F n = σ(Z 0 , X i , i ≤ n).
Note that δM n is a martingale difference sequence with respect to the filtration F = (F n ). (γ n ) n≥1 is assumed to be a deterministic step sequence decreasing to zero and satisfying:

n≥1 γ n = +∞ and n≥1 γ 2 n < +∞. (2.3)
In the following, we will take γ n = c n γ where c is a positive constant and γ ∈]

1 2 , 1]. Finally, K is an hyperrectangle such that z * is in the interior of K: K = {m ∈ R d , a i ≤ m i ≤ b i } × [0, A].
This assumption can be difficult to verify in specific examples. However, the projection on the compact set K can be replaced by an adaptive projection procedure, the so-called projection "à la Chen", introduced in H. [START_REF] Chen | Stochastic approximation procedures with randomly varying truncations[END_REF], H.-F. [START_REF] Chen | Convergence and robustness of the robbinsmonro algorithm truncated at randomly varying bounds[END_REF]. Given an increasing sequence of compact sets (

K n ) n that covers R d , each time Z n leaves the current compact set K σ(n) , the compact set is replaced by K σ(n)+1
before the next step of the algorithm:

Z n+1 = Π K σ(n) [Z n + γ n H(X n+1 , Z n )], Z 0 = z 0 ∈ K 0 , K σ(n+1) =    K σ(n) if Z n ∈ K σ(n) , K σ(n)+1 if Z n / ∈ K σ(n) .
Lelong (2013) proved Central Limit Theorems for this type of projection, that are similar to the ones with the projection over one single compact set (see also [START_REF] Andrieu | Stability of stochastic approximation under verifiable conditions[END_REF]). For the sake of simplicity, we will use here the classical projection over a fixed compact set K.

Main Results

We start by stating the two main results regarding the consistency and asymptotic normality of Algorithm (2). The proofs of Theorem 2.1 and 2.3 can be found in Section 2.2 and 2.3.

Let σ 2 (•), Σ(•) and m 2+p (•), for p > 0, be defined as follows:

       σ 2 (z) = E[||H(X, z) -h(z)|| 2 ]; m 2+p (z) = E[||H(X, z) -h(z)|| 2+p ]; Σ(z) = E[(H(X, z) -h(z))(H(X, z) -h(z)) ⊺ ].
Recall that Assumption (A ′ l ) ensures the existence of a unique risk allocation z * = (m * , λ * ). We introduce the following classical assumptions for stochastic algorithms in order to obtain the almost sure convergence of (Z n ) n to z * :

(A a.s. ) i. h is continuous on K; ii. sup z∈K σ 2 (z) < ∞.
Theorem 2.1 (Almost sure convergence). Assume that the sequence (Z n ) is defined by the algorithm (2) and that assumptions (A ′ l ) and (A a.s. ) hold. Then, Z n → z * P -almost surely as n → ∞.

We introduce the following assumptions for the central limit theorem:

(A a.n. ) i. m → E[∇l(X -m)] is continuously differentiable. Let A := Dh(z * ) (Jacobian matrix of h at z * ); ii. (H(X n+1 , Z n )1 |Zn-z * |≤ρ ) is uniformly integrable for small ρ > 0;
iii. For some p > 0 and ρ > 0, sup

|z-z * |≤ρ m 2+p (z) < ∞; iv. Σ(•) is continuous at z * . Let Σ * := Σ(z * ) .
Remark 2.2. 1. Assumption (A a.s. ) (i) is satisfied for continuously differentiable loss functions l, by the Lebesgue dominating convergence theorem. In the case of univariate SR (see e.g. [START_REF] Dunkel | Stochastic root finding and efficient estimation of convex risk measures[END_REF]), the first assumption is also trivial since h only depends on the continuous loss function l and not on ∇l.

2. In the particular case when X is bounded, the assumptions (A a.n. )(ii), (iii) and (iv) are straightforward.

Theorem 2.3 (Asymptotic Normality). Recall that the step sequence

(γ n ) n≥0 in (2) is de- fined by γ n = c n γ . Assume that γ ∈ ( 1 2 , 1
) and that assumptions (A ′ l ), (A a.s. ) and (A a.n. ) hold. Then,

√ n γ (Z n -z * ) → N 0, c 2 ∞ 0 e cAt Σ * e cA ⊺ t dt .
If furthermore, cA + I 2 is a Hurwitz matrix and cI -P is positive definite with P solution to the Lyapunov's equation:

A ⊺ P + P A = -I, then for γ = 1, √ n(Z n -z * ) → N 0, c 2 ∞ 0 e (cA+ I 2 )t Σ * e (cA ⊺ + I 2 )t dt .
Remark 2.4.

1. Note that, for convex optimization problems, where the matrix A is symmetric negative definite, the two additional conditions reduce to the classical condition cA+ I 2 is negative definite. Indeed, in this case, the solution of the Lyapunov's equation A ⊺ P + P A = -I is simply P = -A -1 /2 and the condition cI -P is positive definite, becomes equivalent to cA + I 2 is negative definite. 2. From a formal point of view, the choice γ = 1 gives the best rate of convergence. The asymptotic variance in this case depends on the constant c. We need to choose it such that cA + I 2 is a Hurwitz matrix and cI -P is positive definite. Setting c too small may lead to no convergence at all, while setting it too large, may lead to slower convergence as the effects of large noises early in the procedure might be hard to overcome in a reasonable period of time.

3. The learning rate c in γ n = c n can be replaced by a conditioning nonsingular matrix Γ. In this case, the results of Theorem 2.3 still hold, where γ n is replaced by a matrix Γ n = Γ/n, and with the modified asymptotic covariance matrix:

∞ 0 e (ΓA+ I 2 )t ΓΣ * Γ ⊺ e (A ⊺ Γ ⊺ + I 2 )t dt.
In particular, if ΓA is a negative definite matrix, then, by the previous remark, the hypothesis of Theorem 2.3 are verified as soon as ΓA + I 2 s a Hurwitz matrix. Actually, the optimal choice of conditioning matrix Γ, which is also called the gain matrix, is the one that will minimize the trace of asymptotic covariance. This is done by taking Γ = -A -1 which yields the asymptomatic optimal covariance: A -1 Σ * (A -1 ) ⊺ . However, the optimal choice of Γ depends on the function h and the equilibrium point z * which are unknown to us. Adaptive procedures that choose the matrix Γ dynamically by estimating Dh(z * ) adaptively have been suggested in the literature (see for example [START_REF] Ruppert | Stochastic approximation[END_REF]), but are generally not as efficient as the Polyak-Ruppert averaging estimators discussed in Section 2.4.

Proof of Theorem 2.1

The proof of the almost sure convergence of (Z n ) to z * rely on the so-called ODE method for the study of constrained RM algorithms. Introduced by Ljung and Kushner, the main idea of the method consists in replacing the analysis of the sequence (Z n ) by the analysis of continuous-time interpolations of the sequence of stochastic approximations. As detailed in Theorem 4.8 of Appendix 4.2, these interpolations converge under suitable assumptions to a solution of the projected ODE defined by :

ż = h(z) + C(z), C(z) ∈ -C(z),
(2.4)

where C(z) is the convex cone determined by the outer normals to the faces that need to be truncated at z, and C(z) is the minimum force needed to bring back z to K. In particular,

C(z(t)) ̸ = 0 only for z(t) ∈ ∂K. When z(t) ∈ ∂K with h(z(t)) directed outside of K, then C(z(t)) take the "smallest" value in -C(z) to maintain z(•) in K. More precisely, for z ∈ R d+1 and 1 ≤ i ≤ d: C i (z) = -h i (z i ) 1 {z i =a i }∪{h i (z)<0} + 1 {z i =b i }∪{h i (z)>0} .
Note that the existence of the projected ODE is obtained as a byproduct of the convergence of the SA algorithm (see Theorem 4.8 in Appendix). We refer to Section 4.3 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] for more details on the projected ODE, or to Section 2.2 of [START_REF] Nagurney | Projected dynamical systems and variational inequalities with applications[END_REF] for a more general presentation of projected ODEs.

Theorem 5.2.1 of [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] (Theorem 4.8) yields that for fixed ω ∈ Ω, the constrained sequence (Z n (ω)) defined by ( 2) converges to a point in the limit set of the projected ODE (2.4). In order to prove Theorem 2.1, the main point is thus to show that the limit set of projected ODE is reduced to z * , or in other words, that z * is globally asymptotically stable in the sens of Lyapunov. This point, which is straigthforward in the univariate case, is more complicated for MSRM due to the introduction of the Lagrangian function L from which h derives.

Since z * is in the interior of K and h(z * ) = 0, z * is an equilibrium point for the projected ODE (2.4). In order to study the asymptotic stability of the equilibrium z * , we define in the following proposition some convenient Lyapunov function V . A natural and classical choice for this type of problems is V (z) = ||z -z * || 2 . Some recalls on Lyapunov stability and Lyapunov functions are presented in Section 4.1 (see also e.g. [START_REF] Slotine | Applied nonlinear control[END_REF] for more details). The following proposition shows that its derivative along any state trajectory is negative semi-definite on K.

Proposition 2.5. The function

V (z) = ||z -z * || 2 is such that z → V (z) = ⟨∇V (z), h(z) + C(z)⟩ is negative semi-definite on K with the respect to the ODE in (2.4). Furthermore, V (z) ≤ -W (z) := 2⟨z -z * , h(z)⟩ ≤ 0, ∀z ∈ K. (2.5) Proof. First, let z = (m, λ) ∈ int(K) so that V (z) = ⟨∇V (z), h(z)⟩ = 2⟨z -z * , h (
z)⟩, and define L the Lagrangian as defined in (1.2). We have:

⟨z -z * , h(z)⟩ = ⟨m -m * , λE[∇l(X -m)] -1⟩ + (λ -λ * )E[l(X -m)] = -⟨m -m * , ∂ m L(m, λ)⟩ + (λ -λ * )∂ λ L(m, λ).
Now, thanks to the convexity of L with respect to m, we have:

L(m * , λ) ≥ L(m, λ) + ⟨m * - m, ∂ m L(m, λ)⟩. This in turn implies that L(m * , λ) -L(m, λ) + (λ -λ * )∂ λ L(m, λ) ≥ ⟨m * -m, ∂ m L(m, λ)⟩ + (λ -λ * )∂ λ L(m, λ).
But, we also have,

L(m, λ) = m i + λE[l(X -m)] = m i + λ * E[l(X -m)] + (λ -λ * )E[l(X -m)] = L(m, λ * ) + (λ -λ * )∂ λ L(m, λ).
The previous inequality becomes then

⟨m -m * , ∂ m L(m, λ)⟩ -(λ -λ * )∂ λ L(m, λ) ≥ L(m, λ * ) -L(m * , λ)
The RHS of the last inequality is non-negative, because, (m * , λ * ) is a saddle point, that is

L(m * , λ) ≤ L(m * , λ * ) ≤ L(m, λ * )
. Therefore, we get that, ⟨z -z * , h(z)⟩ ≤ 0.

(2.6)

Moreover, since L is strictly convex with respect to m, the inequality (2.6) is strict if m ̸ = m * .

Finally, let z ∈ ∂K. If z i = b i (resp. a i ) and h i (z) > 0 (resp. h i (z) < 0) for some i, then

C i (z) = -h i (z). Since z * ∈ int(K), a i < z * i < b i and thus (z i -z * i )C i (z) ≤ 0. Hence, V (z) = 2⟨z -z * , h(z) + C(z)⟩ ≤ ⟨z -z * , h(z)⟩ ≤ 0, which achieves the proof. Remark 2.6. We cannot conclude that V is negative definite on K, because z ̸ = z * does not imply that m ̸ = m * . Besides, if z = (m * , λ) such that λ ̸ = λ * , we have V (z) = 0 and z ̸ = z * .
As stated in the previous remark, V is not negative definite, and thus the standard Liapunov theorem (see e.g. Theorem 3.3 in [START_REF] Slotine | Applied nonlinear control[END_REF]) cannot be applied. In order to prove the asymptotic stability of z * , we apply LaSalle's invariant principle. The invariant set theorem characterizes the limit set of the ODE as the largest invariant set in {z; V (z) = 0}. Note that, due to projection on K, the r.h.s of the ODE 2.4 is discontinuous.

Versions of Lasalle's invariance principle for very general nonsmooth ODEs can be found in [START_REF] Fischer | Lasalle-yoshizawa corollaries for nonsmooth systems[END_REF] (see also e.g. [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF], [START_REF] Shevitz | Lyapunov stability theory of nonsmooth systems[END_REF]). In the case of (2.4), the direct proof presented below is facilitated by the fact that K is a compact set.

Proposition 2.7. The equilibrium point z * of the ODE (2.4) is globally asymptotically stable.

Multivariate Systemic Risk Measures and Stochastic Algorithms

Proof. Let us first show that limit points of the ODE (2.4) are in the set R = {z ∈ K; W (z) = 0}. This can viewed as special case of Corollary 1 in [START_REF] Fischer | Lasalle-yoshizawa corollaries for nonsmooth systems[END_REF]. Let (z(t)) be a solution to the projected ODE (2.4) with initial condition z(0) ∈ K.

The first step of the proof is proving that W (z(t)) → 0 when t → ∞, with W the scalar function defined by W (z) = 2⟨z -z * , h(z)⟩, by applying the Barbalat's Lemma 4.6. By Proposition 2.5, V ≤ 0. Hence, t → V (z(t)) is decreasing and

t 0 W (z(s))ds ≤ - t 0 V (z(s))ds = V (z(0)) -V (z(t)).
(2.7)

Since V (z(t)) ≥ 0, the l.h.s of the previous inequality is increasing and bounded by V (z(0)), and thus lim t→∞ t 0 W (z(s))ds < ∞. Since z is absolutely continuous and h is bounded on K, z is uniformly continuous. Finally, since W is a continuous function, it is uniformly continuous on the compact K and thus t → W (z(t)) is uniformly continuous. By Barbalat's lemma, we can conclude that lim t→∞ W (z(t)) = 0. Thus, limit points of the ODE (2.4) are in the set R = {z ∈ K; W (z) = 0}.

Since solutions of (2.4) are in the bounded domain K, the limit set of the projected ODE is an invariant set (see e.g. Section 4.3 of [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] or Lemma 12.5 in [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF]). It thus remains to show that {z * } is the only invariant set in R. Let M be the largest invariant set in R, and z = (m, λ) ∈ M . As discussed in the proof of Proposition 2.5, if (m, λ) ∈ R, then necessarily m = m * . Thus, since

M is an invariant set in R, if z(0) = (m, λ) ∈ M , then z(t) = (m * , λ(t)) ∈ R for all t ≥ 0, with z(•) solution of the following ODE,      0 = λ(t)E[∇l(X -m * )] -1 + C(z(t)) m , t ≥ 0, dλ(t) dt = E[l(X -m * )] + C(z(t)) λ , t ≥ 0. (2.8) Now, since z * = (m * , λ * ) ∈ int(K), we get that, C(z(t)) m = C((m * , λ(t)) m = 0. We obtain that λ(t)E[∇ m l(X -m * )
] -1 = 0, and consequently t → λ(t) is a constant function equal to λ * , since λ * is the only solution of the first order condition by Theorem 1.6.

We have then showed that the limit set is simply {z * } and therefore z * is globally asymptotically stable equilibrium for the ODE (2.4).

Proof of Theorem 2.1. By Proposition 2.7, z * is the only point in the limit set of ODE (2.4).

Thus, if the assumptions 4.7 (1)-( 3) are verified, we can conclude that Z n → z * as n → ∞ by Theorem 4.8 in the Appendix. The first assumption (1) of the theorem is exactly the Assumption (A a.s. )-(ii). Assumption (2) corresponds to the continuity of h, i.e. Assumption (A a.s. )-(i). Finally, Assumptions (3) corresponds to the hypothesis (2.3) on the step sequence (γ n ), which concludes the proof.

Proof of the asymptotic normality

The proof of Theorem 2.3 is obtained as an application of Theorem 4.10 in the Appendix (Theorem 10.2.1 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]), with g = ch and Y n = cH(X n+1 , Z n ) and

ϵ n = 1 n γ .
Proof of Theorem 2.3. By the proof of Theorem 2.1, Assumption 4.7 is satisfied. Let us prove that Assumption 4.9 holds.

First, let us start with the case γ ∈ ( 1 2 , 1). Then Assumption 4.9.1 is true since

µ n = 1 + γ 2n + o( 1 n 2
). Assumption 4.9.2 corresponds to the assumption (A a.n. )-ii.. Let (Z n (•)) be the sequence of continuous time processes defined by (4.2). For all t ≥ 0 and n ∈ N, Z n (t) ∈ {Z p ; p ≥ n} ⊂ K. Hence, the sequence (Z n (•)) is uniformly bounded. Furthermore, since (Z n ) converges a.s. to z * by Theorem 2.1, Aldous' criterion for tightness (see e.g. Theorem VI.4.5 in Jacod and Shiryaev ( 2013)) is satisfied and hence the family processes (Z n (•)) is tight. Finally, for all t ≥ 0, Z n (t) = Z 0 (t n + t) converges a.s. to the constant z * by definition of the interpolation process Z 0 . Thus, (Z n (•)) converges weakly to z * , and Assumption 4.9.3 holds. Assumption 4.9.4 corresponds to (A a.n. )-i.. The first part of 4.9.6 is (A a.n. )-iii.. By (A a.n. )-iv.

(continuity of Σ at z * ) and since Z n → z * a.s.,

E[(Y n -g(Z n ))(Y n -g(Z n )) ⊺ 1 |Zn-z * |≤ρ ] → Σ 1 := c 2 Σ * .
Finally, 4.9.7 holds for ϵ n = 1 n γ (see e.g. p.341 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]). It remains to show that 4.9.5 holds, that is the matrix A ′ = cA is a Hurwitz matrix. We have:

A =       λ * DE[∇l(X -m * )] E[∇l(X -m * )] -E[∇l(X -m * )] 0       = -       Â -1 λ * 1 λ * 0      
, where  := -λ * DE[∇l(X -m * )] corresponds to the second derivative of the Lagrangian L with respect to m. Note that L is strictly convex with respect to m due to the strict

convexity of m → E[l(X -m)].
This implies that  is a positive definite matrix. Thanks to Theorem 3.6 in [START_REF] Benzi | Numerical solution of saddle point problems[END_REF], we deduce that A (and thus A ′ ) is a Hurwitz matrix.

For the case γ = 1, it remains to prove Assumption 4.9.7, ie that λ > 1 where λ is the largest real number such that I ≥ λP ′ , with P ′ the solution of the Lyapunov equation

(A ′ ) ⊺ P ′ + P ′ A ′ = -I. Recalling A ′ = cA, P ′ is also solution of A ⊺ cP ′ + cP ′ A = -I. By
uniqueness of the Lyapunov equation's solution, we obtain that P = cP ′ , with P as defined in Theorem 2.3. Since cI -P = c(I -P ′ ) is positive definite by assumption, I -P ′ is positive definite and thus λ > 1, which achieves the proof.

Polyak-Ruppert Averaging

In order to ease the tuning of the step parameter which is known to monitor the numerical efficiency of RM algorithms, we are led to modify our algorithm and to use an averaging procedure. Averaging algorithms were introduced by Ruppert (see [START_REF] Ruppert | Stochastic approximation[END_REF]) and Polyak (see [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF]) and then widely investigated by many authors. [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] and [START_REF] Kushner | Stochastic approximation with averaging and feedback: Faster convergence[END_REF] studied these algorithms in combination with projection and proved a Central Limit Theorem (CLT) for averaging constrained algorithms.

The following theorem describes the Polyak-Ruppert algorithm for MSRM and states its asymptotic normality. It is a direct consequence of theorem 11.1.1 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF].

Theorem 2.8. Assume γ ∈ ( 1 2 , 1) and that assumptions (A ′ l ), (A a.s. ) and (A a.n. ) hold. For any arbitrary t > 0, we define

Zn = γ n t n+t/γn-1 i=n Z i , (2.9)
where any upper summation index u ∈ R + is interpreted as its integer part. If Σ * is positive definite, then we have the following CLT:

t γ n Zn -z * -→ n→∞ N 0, V + O 1 t , (2.10)
where

V = A -1 Σ * (A -1 ) ⊺ .
Remark 2.9.

1. In (2.9), the window of averaging is t/γ n for any arbitrary real t > 0. Equivalently, γ n × (size of window) does not go to infinity as n → ∞, hence the name "minimal window" of averaging. In contrast, the "maximal window" of averaging corresponds to a window size q n t such that γ n q n → ∞. A natural and a classical choice is taking γ n = c/n γ and q n = n. In the case of maximal window of averaging, under some extra conditions, we are able to achieve the optimal asymptotic variance without an extra term O(1/t)(see Theorem 11.3.1 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]).

2. Two sided averages can also be used instead of the one-sided average in (2.9).

Estimator of asymptotic variance

The previous CLT theorems assert that, under some suitable conditions, our RM and PR algorithms converge to the root z * with a corresponding rate. More specifically, in Theorem 2.8, the asymptotic variance V depends on Σ * and A. In practice, these two quantities are unknown and need to be approximated in order to derive confidence intervals for our estimators. In Theorem 2.3, in both cases, γ = 1 and γ ∈ ( 1 2 , 1), the asymptotic variance is expressed as an infinite integral that involves Σ * and A. The numerical evaluation of these integrals is a non-trivial exercise even when Σ * and A are known. In [START_REF] Hsieh | Confidence regions for stochastic approximation algorithms[END_REF], they described an approach that produces confidence 2. Multivariate Systemic Risk Measures and Stochastic Algorithms regions and that avoids the necessity of having to explicitly estimate these integrals.

In the following proposition, we provide consistent estimators of these two quantities. The proof relies mainly on the Martingale Convergence Theorem.

Proposition 2.10. Assume (A ′ l ), (A a.s. ) and (A a.n. ) hold. If z → E[||H(X, z)|| 4 ] is locally bounded around z * , then,

Σ n := 1 n n k=1 H(X k , Z k-1 )H(X k , Z k-1 ) ⊺ → Σ * a.s.
(2.11)

Let A ϵ n the matrix whose elements A ϵ n (i, j) for i, j ∈ {1, ..., d + 1} are defined as follows:

A ϵ n (i, j) := 1 ϵn n k=1 H i (X k , Z k-1 + ϵe j ) -H i (X k , Z k-1 ), then, lim ϵ→0 lim n→∞ A ϵ n = A a.s.
(2.12)

Proof. Let (S n ) n∈N * be the sequence defined as:

S n = H(X n , Z n-1 )H(X n , Z n-1 ) ⊺ -Σ(Z n-1 ) -h(Z n-1 )h(Z n-1 ) ⊺ , n ≥ 1 (S n ) n∈N * is
a martingale difference sequence adapted to F and consequently the following sequence (M n ) n∈N * defined as:

M n = n i=1 S i i , n ≥ 1,
is a F-martingale. Moreover, since Z n → z * a.s., h and Σ are continuous in z * (Assumptions (A a.s. )-i. and (A a.n. )-iv.), and by the boundedness of z → E[||H(X, z)|| 4 ] around z * , we have that:

sup n≥1 E[||S n || 2 |F n-1 ] < ∞ a.s.
Thus, the martingale convergence theorem ensures the existence of a finite random variable

M ∞ such that M n → M ∞ a.s. Kronecker's lemma then guarantees that 1 n n i=1 S i → 0. Now, since, Σ n = 1 n n i=1 S i + 1 n n i=1 Σ(Z i-1 ) + 1 n n i=1 h(Z i-1 )h(Z i-1 ) ⊺ , we deduce that Σ n → Σ * .
The proof of (2.12) follows using the same arguments above.

Remark 2.11.

1. Instead of averaging on all observations, one could modify the estimators above and average only on recent ones. This might improve the behaviour of these estimators.

If we denote

V n := (A ϵ n ) -1 Σ n ((A ϵ n ) -1
) ⊺ , then we obtain an approximate confidence interval for PR estimator with a confidence of 1 -α in the following form:

Zj,n - V jj,n tn γ q α , Zj,n +
V jj,n tn γ q α , j ∈ {1...d}, γ ∈ (0, 1), (2.13)

where q α is the 1 -α 2 quantile of a standard normal random variable. Note that this confidence interval has the advantage of being obtained with one simulation run. For RM estimators, confidence intervals could be estimated empirically, however a naive approach would multiple runs of the SA algorithm.

Variance reduction with adaptive importance sampling (IS)

The convergence of the stochastic algorithms presented in this section can be improved by implementing variance reduction techniques such as important sampling (IS). The original idea of IS for the Monte-Carlo computation of an expectation E[G(X)], is to find a change of probability measure

P ′ , such that E[G(X)] = E ′ [G ′ (X)] and Var P ′ (G ′ (X)) < Var(G(X)).
The main issue is then to find the right change of measure in order to reduce the variance. A possible parametric approach is to look for the optimal parameter θ * minimizing the variance, among a family of distribution (P θ ). Under some assumptions, the minimization problem can be reinterpreted as finding the root of a function written as an expectation. This is the main idea of recursive IS by stochastic algorithm, where the parameter θ * is estimated by a RM algorithm.

The Esscher transform (or exponential tilting) is a natural tool for parametric IS, used extensively in actuarial science and for variance reduction. This transfom is a natural change of measure since it corresponds to a translation in the gaussian case, and a change of intensity for Poisson processes. Let us first recall the main idea behind the Esscher transform and recursive IS in the standard case. For more details, see e.g. see e.g. [START_REF] Arouna | Adaptative monte carlo method, a variance reduction technique[END_REF], [START_REF] Kawai | Optimal importance sampling parameter search for lévy processes via stochastic approximation[END_REF], or [START_REF] Lemaire | Unconstrained recursive importance sampling[END_REF].

Esscher transform and parametric IS In the following, we assume that the R dvalued random variable X has a density f . For θ ∈ R d , let ψ be the log-Laplace transform of X, assumed to be finite:

ψ(θ) = log E[exp(⟨θ, X⟩] < ∞.
Then, ψ is convex and infinitely differentiable. The Esscher transform is defined by the following change of measure introducing the density g θ :

g θ (x) = exp (⟨θ, x⟩ -ψ(θ)) f (x).
(2.14)

In the following we denote by X (θ) the r.v. of density g θ . The expectation can be replaced by

E[G(X)] = E[G(X (θ) ) exp -⟨X (θ) , θ⟩ + ψ(θ) ].
The main idea of the recursive IS procedure is to choose θ in order to minimize the variance of G(X (θ) ) exp -⟨X (θ) , θ⟩ + ψ(θ) , which boils down to minimizing the function

E[G 2 (X (θ) ) exp -2⟨X (θ) , θ⟩ + 2ψ(θ) ] = E[G 2 (X) exp (-⟨X, θ⟩ + ψ(θ))].
In other words, the optimal parameter is

θ * = argmin E[G 2 (X, z) exp (-⟨X, θ⟩ + ψ(θ))].
(2.15)

Under some regularity assumptions, the minimization problem can be solved using stochastic algorithms.

Adaptative RM algorithm with IS Our setting of multivariate systemic risk measures is more complex. For 1

≤ i ≤ d + 1, let H i be the i component of the function H defined by (2.1). Then, for z = (m, λ) ∈ R d × [0, ∞[ and θ ∈ R d , h i (z) = E[H i (X, z)] = E[H i (X (θ) , z) exp -⟨X (θ) , θ⟩ + ψ(θ) ],
(2.16)

Let V = (V i ) 1≤i≤d , with V i the function defined for z ∈ R d × [0, ∞[ and θ ∈ R d by V i (θ, z) = E[H 2 i (X (θ) , z) exp -2⟨X (θ) , θ⟩ + 2ψ(θ) ] = E[H 2 i (X, z) exp (-⟨X, θ⟩ + ψ(θ))].
(2.17)

In order to reduce the asymptotic variance obtained in Theorem 2.8, it is natural to take

θ i = θ * i in (2.16) such that θ * i = argmin V i (θ, z * ) = argmin E[H 2 i (X, z * ) exp (-⟨X, θ⟩ + ψ(θ))]. (2.18) Assume that for each 1 ≤ i ≤ d, E[|X|H 2 i (X, z) exp (-⟨X, θ⟩ + ψ(θ))] < ∞, ∀z ∈ R d × [0, ∞[ , θ ∈ R d . Then, V(•, z) is differentiable and for 1 ≤ i ≤ d, ∇V i (θ, z * ) = E[(∇ψ(θ) -X)H 2 i (X, z * ) exp (-⟨θ, X⟩ + ψ(θ))]. (2.19)
Lemaire and Pagès (2010) (Proposition 3, see also [START_REF] Bardou | Computing VaR and CVaR using stochastic approximation and adaptive unconstrained importance sampling[END_REF]) give sufficient conditions for V i (•, z * ) to be convex and for the set {θ ∈ R d ; ∇V i (θ, z * ) = 0} to be non empty.

the VaR/CVaR, the expression of the matrix A -1 is explicit, and this might be the case for MSRM for specific loss functions. Another potential drawback of the adaptative RM with IS (2.21) is the number d(d + 1) of additional parameters introduced by the variance reduction procedure. Once again, the number of parameters could be significantly reduced depending on the loss function l.

Remark 2.12. In dimension 1, it is straightforward to show that θ * minimizes V (θ). Hence the adaptative algorithm (2.21) also provides a general variance reduction method for univariate shortfall risk measures. This is very useful when the optimal change of measure is not known beforehand, which is the case most of the time.

Numerical examples

In this section, we test the performance of the proposed stochastic algorithms schemes for MSRM. We provide numerical examples showing the consistency properties of the different estimators and then their normal asymptotic behaviour, with (PR algorithm) and without averaging (RM algorithm). Two examples are considered. In the first one, we consider a loss function of an exponential type coupled with a normal distribution. This example is relevant for our numerical analysis as we can explicitly express the optimal allocations in a closed form. In the second example, we consider a loss function that involves positive part function with a Gaussian and a compound Poisson distributions. All the algorithms have been implemented on a standard computer using the programming language Python, and are publicly available at https://github.com/AchrafTamtalini/SRM/.

In the following, n will denote the number of steps in one simulation run and N the number of simulations. We introduce the following sequences:

Dn := √ tn γ ( Zn -z * ) γ ∈ ( 1 2 , 1), (3.1)
D n := √ n γ (Z n -z * ), γ ∈ ( 1 2 , 1]. (3.2)

Toy example

As a first simple example, we will consider a exponential loss function of the following form:

l(x 1 , ..., x d ) = 1 1 + α d i=1 e βx i + αe β d i=1 x i - α + d α + 1 (3.3)
We will set d = 2 and consider a bivariate normal vector X = (X 1 , X 2 ) ∼ N(0, M ) with

M = σ 2 1 ρσ 1 σ 2 ρσ 1 σ 2 σ 2 2
. α is a systemic weight parameter taken to be non negative and β > 0 is the risk aversion coefficient. In this case, we can explicitly solve the first order conditions and derive closed formulas for optimal allocations. This will be useful to test our algorithms. The following Lemma extends Example 3.12 in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF]:

Lemma 3.1. Let l be the loss function defined by (3.3), and X the vector of random losses as above. Then, the optimal risk allocation is:

m * i =        βσ 2 i 2 , if α = 0, βσ 2 i 2 + 1 β SRC(ρ, σ 1 , σ 2 , α, β), if α > 0. (3.4) With SRC(ρ, σ 1 , σ 2 , α, β) = ln αe ρβ 2 σ 1 σ 2 -1 + 1 + α(α + 2)e ρβ 2 σ 1 σ 2 .
This shows that, in the case α > 0, the risk allocations are disentangled into two components: an individual contribution βσ 2 i 2 and the Systemic Risk Contribution term (SRC). Note that taking α → 0 makes the SRC null since the systemic weight α is responsible for the systemic contribution in the loss function l. One can also show, by easy calculations, that the SRC is increasing with respect to ρ: the higher the correlation is, the more costly the acceptable monetary allocations are. This could be explained by the fact that, with a higher correlation between the two components, the losses of one will induce the loss of the other and consequently the system will become riskier. Note also that we could also express in a closed form the Jacobian matrix A and Σ * .

Proof. The optimal risk allocations m * i are characterized by the first order conditions given in Theorem 1.6, i.e.

     λ * 1 + α E e β(X i -m * i ) + αe β(X 1 +X 2 -m * 1 -m * 2 ) = 1, i = 1, 2, 1 1 + α E e β(X 1 -m * 1 ) + e β(X 2 -m * 2 ) + αe β(X 1 +X 2 -m * 1 -m * 2 ) - α + 2 α + 1 = 0. The two first equations imply that E[e β(X 1 -m * 1 ) ] = E[e β(X 2 -m * 2 )
], which in turn gives that,

βσ 2 1 2 -m * 1 = βσ 2 2 2 -m * 2 . (3.5)
The third equation gives e

β 2 σ 2 1 2 -βm * 1 + e β 2 σ 2 2 2 -βm * 2 + αe ρβ 2 σ 1 σ 2 + β 2 2 (σ 2 1 +σ 2 2 )-β(m * 1 +m * 2 ) = 2 + α.
Thanks to (3.5) and denoting

Q = e β 2 σ 2 i 2 -βm * i , we get that, αe ρβ 2 σ 1 σ 2 Q 2 + 2Q -(2 + α) = 0.
Taking the positive solution of the last equation gives Q

= -1+ 1+α(α+2)e ρβ 2 σ 1 σ 2 αe ρβ 2 σ 1 σ 2 . Now, denoting by SRC = -log(Q), we obtain that m * i = βσ 2 i 2 + 1 β SRC.
In all this example, we fix α = 1, β = 1 and σ 1 = σ 2 = 1. With ρ ∈ {-0.5, 0, 0.5}, we obtain the exact values in the table below. Note that since we have X 1 ∼ X 2 ∼ N(0, 1) and

l is permutation invariant, it follows that m * 1 = m * 2 .
Table 1: Exact optimal risk allocations. For RM/PR algorithms, we used a number of steps n = 10 5 . As for the compact K,

we took K = [0, 2] 3 and Z 0 was taken uniformly on K. We run the different algorithms for γ = 1 and γ = 0.7. We chose an averaging parameter t = 10 and we set c = 2 in a first step.

Figure 1 shows that, for different values of ρ ∈ {-0.5, 0, 0.5}, our RM algorithm with γ = 1 converges relatively quickly to the optimal allocations, whereas when γ = 0.7, noise is still persisting. This is due to the step parameter c as discussed in the previous section. In order to get a smoother numerical behaviour, two solutions are available to us: either we use PR averaging (c.f. Figure 1), or we reduce the value of the parameter c (Figure 2). Figure 3 shows the increase of computational time when the error become smaller. We have also compared in Table 2 the computational time (CT) of the RM algorithms (with c = 2 and γ = 1) with the Monte Carlo method used in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF]. As it is often the case, the RM algorithm is slower than the Monte Carlo method in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF].

However, as stated above, the latter does not provide any convergence results, and might not converge in practical applications. Note that we can easily verify that all conditions in (A a.s. ) and (A a.n. ) hold. We can also verify thanks to the exact formula of Σ * , that this matrix is positive definite for the different values of ρ used. This is a condition needed in Theorem 2.8.

For any random estimator, constructing confidence intervals is important to assess the error in the estimation. For PR estimator, confidence interval can be obtained in one simulation run after estimating matrices Σ * and A and hence the asymptotic variance matrix V . Figure 4 shows the convergence, in the case ρ = 0, of the estimator of

V n = A -1 n Σ n (A -1 n ) ⊺
where A n and Σ n are as introduced in Proposition 2.10. In the following table, we give the estimated confidence interval for PR estimator with a confidence coefficient of 95%: As for RM estimators, it is difficult to estimate the asymptotic covariance matrix due to its complexity. In order to visualize the normal behaviour of these estimators, we give the

Numerical examples

empirical probability density function (EPDF) in both cases γ = 1 and γ = 0.7. To this end, we use again a number of steps n = 100000 and we repeat the procedure N = 10000 times. We restrict our attention to the case ρ = 0. Figure 5 shows that D n,i = √ n γ (m n,i -m * i ), i ∈ {1, 2} are very close to a normal distribution. From the two figures above, the width of the 90% confidence interval of the RM estimator for the case γ = 0.7 is approximately 8% and for the case γ = 1 is roughly 2%.

Second example

As a second example, we will consider consider the following loss function used in [START_REF] Armenti | Multivariate shortfall risk allocation and systemic risk[END_REF]:

l(x 1 , ..., x d ) = d i=1 x i + 1 2 d i=1 (x + i ) 2 + α i<j x + i x + j .

First case: Gaussian distribution and d = 2

We start by a simple case where we fix d = 2 and use standard two dimensional Gaussian distribution for the loss vector X. We take K = [0, 2] 3 , n = 10 5 , t = 10, α = 1 and c = 6. Again, we compare RM and PR estimators for different values of ρ. Note that all assumptions can be easily verified except the one that requires the matrix cA + I/2 to be a Hurwitz matrix. However, the PR estimator solve this issue.

The following figure 7 allows us to draw the same conclusions as in the previous example:

RM estimator with γ = 1 and PR estimator are better than RM estimator with γ = 0.7. RM estimator with γ = 0.7 is noisy and one can remediate to this by choosing a smaller value of c as we did in the first example. In order to assess the accuracy of our PR estimator, we give the confidence interval with a 95% confidence coefficient (CI), using the estimators of Proposition 2.10. Since m * 1 = m * 2 in this example, we only present in Table 4 the results for the estimator m * 1 of m * 1 . For RM estimator, we plotted the EPDF of D n,i , i ∈ {1, 2} as well as the ECDF of the error m n -m * for the case ρ = 0 and γ = 1. These figures shows that the length of the confidence interval of 90% in the case γ = 0.7 is much higher that in the case γ = 1 (approx 0.2 against 0.04).

Second case: Compound Poisson Distribution and higher dimensions

In this section, we propose to use compound Poisson processes to model the loss vector X. The scope of application of compound Poisson processes is very wide. It ranges from statistical physics and biology to financial mathematics. In biology, they are used to study dynamics of populations. In the modern financial modeling, compound Poisson processes are used to describe dynamics of risk factors such as interest rates (see for instance [START_REF] Li | Stochastic interest model based on compound poisson process and applications in actuarial science[END_REF]), foreign exchange rates and option pricing (see [START_REF] Jaimungal | Catastrophe options with stochastic interest rates and compound poisson losses[END_REF]). In actuarial science, compound processes are extensively used to model claims sizes and to compute the ruin probability, i.e. the probability that the initial reserves increased by premiums received from clients and decreased by their claims, drops below zero.

More precisely, given a final time T , we consider a multivariate Poisson random vector N T = (N T 1 , ..., N T d ), where each N T i ∼ P(λ i T ) and the loss corresponding to the i th component is

X i = N T i k=1 G k i and (G k i )
k is an i.i.d sequence representing the jump sizes and independent of N T i . We will take two examples for the distribution of the jumps sizes: One with a Gaussian distribution and another one with an exponential one. The correlation between the different components of X will be done through the correlations between components of N T . In what follows, we detail the method of generating a multivariate Poisson random vector, N = (N 1 , ..., N d ) with a vector of corresponding intensities (λ 1 , ..., λ d ). To do so, we will use a method that is based on the Gaussian vectors. More precisely, denote η = (η 1 , ..., η d ) to be a Gaussian random vector having a centered normal distribution with correlation matrix R = (ρ kl ) and Φ to be the standard normal cdf. Then, the random vector ξ = (Φ(η 1 ), ..., Φ(η d ))

has a multivariate distribution with standard uniform marginal distributions. Let P λ (x) =

ρ * kl = E(ζ k ζ l ) -E(ζ k )E(ζ l ) σ(ζ k )σ(ζ l ) = E(ζ k ζ l ) -λ k λ l √ λ k λ l .
We need to express the expectation E(ζ k ζ l ) as a function of ρ kl . We have:

E(ζ k ζ l ) = E P -1 λ k (Φ(η k ))P -1 λ l (Φ(η l )) = ∞ m=1 ∞ n=1 mn P (ζ k = m, ζ l = n) = ∞ m=1 ∞ n=1 mn P (u k m-1 ≤ Φ(η k ) ≤ u k m , u l n-1 ≤ Φ(η l ) ≤ u l n ),
where u i j = P λ i (j), i, j = 1, ..., d. It remains to explicit the probabilities in the last equality. If we denote Φ 2 (•, •, ρ kl ) the bivariate Normal distribution function, we get finally that,

Z mn (ρ kl ) := P (u k m-1 ≤ Φ(η k ) ≤ u k m , u l n-1 ≤ Φ(η l ) ≤ u l n ) = Φ 2 (A k m , B l n , ρ kl ) -Φ 2 (A k m-1 , B l n , ρ kl ) -Φ 2 (A k m , B l n-1 , ρ kl ) + Φ 2 (A k m-1 , B l n-1 , ρ kl )
where A k m = Φ -1 (P λ k (m)) and B l n = Φ -1 (P λ l (n)). As a conclusion, we obtain,

∞ m=1 ∞ n=1
mnZ mn (ρ kl ) = λ k λ l + ρ * kl λ k λ l .

(3.6)

The equation (3.6) gives an implicit relation between ρ * kl and ρ kl . It also involves two infinite sums which makes it hard to solve. In practice, one needs to truncate this sum and choose some appropriate upper-limits M * and N * . We are then able to compute the elements of the correlation matrix ρ of the Gaussian vector given the correlation matrix ρ * of the vector N .

However, there is a problem of sufficient conditions for a given positive semi-definite matrix to be a correlation matrix of a multivariate Poisson random vector. This issue is tackled in [START_REF] Griffiths | Aspects of correlation in bivariate poisson distributions and processes[END_REF] where it is shown that each ρ * kl has to be in a certain range, The following figure shows the covariance matrix for the loss vector X of dimension d = 10 obtained by generating a random correlation matrix (ρ kl ) and using a Gaussian distribution for the jump sizes, i.e. G k i ∼ N(1, 1). The intensity vector was taken uniformly in [1, 3] 10 .

-1 < ρ min kl ≤ ρ * kl ≤ ρ max

Numerical examples

the same optimal allocations. This is the case for components 4 -5 and components 9 -10.

However, once α is taken non null, we see that the same components have no longer the same optimal allocations. For instance, component 10 has higher optimal allocation than 9 when α = 1. This could be explained by the fact that component 10 is more correlated with other components that have high variances, such as components 4 and 5, than component 9. We now consider an exponential distribution for jump sizes as a second example, i.e. G k i ∼ E(a i ). The parameters a i were generated randomly in [0.

Appendix

well-chosen function (most of the time W = V ). We recall below Barbalat's Lemma (see e.g. Lemma 3 in [START_REF] Fischer | Lasalle-yoshizawa corollaries for nonsmooth systems[END_REF]), with is central to proving LaSalle's invariance theorem.

Lemma 4.6 (Barbalat's Lemma). Let ϕ : R + → R be a uniformly continuous function.

Suppose that lim t→∞ t 0 ϕ(s)ds < ∞. Then, lim t→∞ ϕ(t) = 0.

4.2

Strong and weak convergence of constrained stochastic algorithms Let g : R n → R n be a continuous function, and let z * ∈ R n be an equilibrium point of h, i.e.

h(z * ) = 0.

The constrained stochastic approximation algorithm is defined by Continuous time interpolation let t 0 = 0 and for n ≥ 1, t n = n-1 i=0 ϵ i . The piecewise constant continuous-time interpolation of the sequence (Z n ) is defined by in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]. Note that the framework presented here is restricted to the case where E[Y n |F n ] = g(Z n ) (g does not depend on n). Hence, some assumptions are either trivial and not rewritten here or simplified.

Z n+1 = Π K [Z n + ϵY n ], ( 4
Z 0 (t) = Z n , ∀ t

Strong convergence

The assumption below can be found in Section 5.2 of [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF].

  We denote by || • || the Euclidean norm, and x • y = x k y k . For a function f : R d → [-∞, ∞], we denote by f * (y) = sup x {x • y -f (x)} the convex conjugate of f . The space L 0 (R d ) inherits the lattice structure of R d and therefore, we can use the classical notations in R d in a P -almost-sure sens.

Figure 1 :Figure 2 :

 12 Figure 1: Consistency of RM/PR estimators with for different values of ρ.

Figure 3 :

 3 Figure 3: Computational time (in seconds) as a function of the (average) absolute error (log log scale), for the RM algorithm (γ = 1, c = 2 and ρ = 0.5).

  Comparison of computation times for RM algorithms (γ = 1 and c = 2) and Monte Carlo Method in Armenti et al. (2018) (m * 1 = m * 2 = 0.6364).

Figure 4 :

 4 Figure 4: Convergence of the estimator V n .

Figure 5 :Figure 6 :

 56 Figure 5: Empirical cumulative density function of m n -m * .

Figure 7 :

 7 Figure 7: Consistency of RM/PR estimators with for different values of ρ.

Figure 8 :Figure 9 :

 89 Figure 8: Empirical cumulative density function of m n -m * .

  j /j!)e -λ be the cdf of the Poisson distribution with parameter λ. Now, consider the vector ζ = (ζ 1 , ..., ζ d ) where ζ k = P -1 λ k (Φ(η k )), k = 1, ..., d. ζ has therefore Poisson marginal distributions with intensities (λ 1 , ..., λ d ). We can express the correlation coefficient ρ * kl = corr(ζ k , ζ l ) as a function of ρ kl = corr(η k , η l ):

56

  Algorithm for generating a sample of X with Compound Poisson Distribution Input: (λ i ) i=1,...,d , (σ i ) i=1,...,d , T , and ρ * correlation matrix of N T ; Ensure: For each k > l, ρ * kl verifies the inequality in (3.7); 1 Solve the equation (3.6) to find ρ the correlation of the Gaussian vector η = (η 1 , ..., η d );2 Generate a sample of Gaussian vector η with correlation matrix ρ and for i = 1, ..., d;3 for i = 1, ..., d do 4 Set N T i = P -1 λ i T (Φ(η i )); Generate a i.i.d sample of G k i of size N T i ; Set X i = k=1 G k i . Output: X;

Figure 12 :

 12 Figure 12: Correlation matrix of the vector loss X with exponential distribution for the law of jump sizes.

  .1) with (Y n ) a sequence of R n random variables such that E[Y n |F n ] = g(Z n ), Π K the projection operator on the constraint set K = {z; a i ≤ z i ≤ b i } ⊂ R n and (ϵ n ) a sequence decreasing to zero, such that ∞ n=0 ϵ n = ∞.

  Assumption 4.7. 1. (A2.1) sup E[|Y 2 n |] < ∞.

Table 3 :

 3 Confidence intervals for PR estimators.

	ρ	CI for m * 1	CI for m * 2
	-0.5	[0.3772, 0.4047]	[0.3679, 0.3949]
	0	[0.4962, 0.5259]	[0.4912, 0.5213]
	0.5	[0.6194, 0.6629]	[0.6203, 0.6665]

Table 4 :

 4 Confidence intervals for PR estimators.

	ρ	m * 1	CI
	-0.5	0.188	[0.1790, 0.2089]
	0	0.21	[0.1963, 0.2303]
	0.5	0.25	[0.2415, 0.2769]

  ∈ [t n , t n+1 [.The shifted process is defined byZ n (t) = Z 0 (t n + t), ∀t ≥ 0. √ ϵ n , ∀ t ∈ [t n , t n+1 [,andU n (t) = U 0 (t n + t), ∀t ≥ 0. (4.3)We recall below the assumption and statement of Theorem 5.2.1 (p.127) and 10.2.1 (p.329)

		(4.2)
	Similarly, let
	U 0 (t) =	Z

n -z *

The authors research is part of the ANR project DREAMeS (ANR-21-CE46-0002) and benefited from the support of the "Chair Risques Emergents en Assurance" under the aegis of Fondation du Risque, a joint initiative by Le Mans University and Covéa.

We also refer to [START_REF] Kawai | Optimal importance sampling parameter search for lévy processes via stochastic approximation[END_REF] in the special case when X = ξ T with ξ a Levy process. A RM algorithm as in (2) can be implemented to find the root of each gradient ∇V i , replacing H i with F i (X, z * , θ) = (∇ψ(θ) -X)H 2 i (X, z * ) exp (-⟨θ, X⟩ + ψ(θ)) .

(2.20)

However such an algorithm cannot be implemented in practice, since z * is unknown. The same issue appears in the case of the VaR/CVaR estimation in [START_REF] Bardou | Computing VaR and CVaR using stochastic approximation and adaptive unconstrained importance sampling[END_REF]. Hence, the authors introduced an adaptive procedure combining the RM algorithm and recursive Importance Sampling, by substituting the unknown parameters at the nth step the by their estimation at the previous step. This idea can be adapted to our framework of multivariate systemic risk measures.

Let us first introduce some notations. For Θ = (θ 1 , . . . , d+1) , let H be the function replacing H in (2.16), defined by

where F i has been defined in (2.20). Finally, we denote by

valued random vector, where each R d component X (θ i ) has the density g θ i and is independent of the others. Finally, let Σ(z, Θ) be the covariance matrix of H(X, z, Θ).

The new stochastic algorithm is defined as follow:

(2.21)

where K ′ is an hyper rectangle containing 0 and Θ * . In the specific case of the VaR-CVaR approximation, [START_REF] Bardou | Computing VaR and CVaR using stochastic approximation and adaptive unconstrained importance sampling[END_REF] propose an unconstrained version of the algorithm, based on [START_REF] Lemaire | Unconstrained recursive importance sampling[END_REF]. However, the procedure cannot be applied here since H does not verify the sublinearity condition, as discussed in Section 2. However, as for K, K ′ can also be replaced by an increasing sequence of compact (K ′ σ(n) ) n∈N . The proof of Theorem 2.1 and 2.3 can be adapted to obtain the a.s. convergence and asymptotic normality of the approximations (Z n , Θ n ) n introduced in (2.21), as well as its Polyak-Ruppert counterpart. In particular, the covariance matrix Σ * = Σ(z * , 0) will be replaced by Σ(z * , Θ * ), minimizing the variance of the vector's H(X, z * , •) components.

The asymptotic variance for the Polyak-Ruppert estimator with adaptative IS is actually

) ⊺ , with A the Jacobian matrix of h at z * . The matrix A does not depend on the Θ * , but the diagonal elements of the covariance matrix V might not minimize the function V (Θ) = A -1 Σ(z * , Θ * )(A -1 ) ⊺ . In the case of the joint computation of The above figure shows that there are components with the same optimal allocations for the case α = 0. This is something we expect to see, since with α = 0, correlations between components are not involved, so components with the same variance should have

Appendix

ODE method and related concepts

In this section, we recall some key concepts of the stability of an ODE ż = h(z).

Concepts of stability of an ODE

A state z * is an equilibrium of the ODE if h(z * ) = 0. To describe the behaviour of the system around the equilibrium, a number of stability concepts are needed. Let us first introduce the basic concepts of stability. To alleviate the notations, we will take 0 as an equilibrium state.

Definition 4.1. The equilibrium z * = 0 is said to be stable, if for any R > 0, there exists r > 0 such that if ||z(0)|| < r, then ||z(t)|| < R for all t ≥ 0. Otherwise, the equilibrium is unstable.

Appendix

Essentially, this means, the system can be kept arbitrarily close to the origin by starting sufficiently close to it. This is also know as Lyapunov stability. In some applications, Lyapunov stability is not enough and the concept of asymptotic stability is needed.

Definition 4.2. An equilibrium point z * = 0 is asymptotically stable if it is stable, and if in addition, there exists some r > 0 such that ||z(0)|| < r implies that z(t) → 0 as t → ∞. The ball B r is called a domain of attraction of the equilibrium point.

If asymptotic stability holds for all initial states, the equilibrium is said to be globally asymptotically stable.

Lyapunov function and invariant set

The first property that need to be verified by a Lyapunov function is positive definiteness.

Definition 4.3. A scalar continuous function V (z) is said to be locally positive definite if V (0) = 0 and in around 0, we have, z ̸ = 0 ⇒ V (z) > 0.

If the above property holds over the whole state space, then V (z) is said to be globally positive definite.

The above definition implies that the function V has a unique minimum at the origin 0. Next, we define the "derivative of V" with respect to time along the system trajectory.

Assuming that V is differentiable, this derivative is defined as:

Definition 4.4 (Lyapunov function). Let V be a positive definite function and continuously differentiable. If its time derivative along any state trajectory is negative semi-definite, i.e.,

then V is said to be a Lyapunov function for the system.

It is possible to draw conclusions on asymptotic stability, with the help of the invariant set theorems introduced by LaSalle. The central concept in these theorems is that of invariant set, a generalization of the concept of equilibrium point.

Definition 4.5. Let z(•) be a solution of some ODE. A set G is said to be an invariant set

For instance, the singleton {z * } where z * is an equilibrium point is an invariant set. Its domain of attraction is also an invariant set. One other trivial invariant set is the whole state-space, ∪

Informally, the Lasalle invariance theorem allows the identification of an ODE's limit set as a subset of the largest invariant set contains in the set {z; W (z) = 0}, with W some 4. Appendix 2. (A2.3) g is continuous.

3. (A2.4) i ϵ 2 i < ∞.

Theorem 4.8 (Theorem 5.2.1 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] (partial)). Let {Z n (•)} be the sequence of continuous-time processes defined by (4.2) and let Assumption (4.7) holds.

There is a set N of probability zero such that for ω / ∈ N , the family of function {Z n (ω, •)} is equicontinuous. Furthermore, if z(•) is a limit of a convergent subsequence, then the function is solution of the projected ODE

and {Z n (ω)} converge to some limit set of the ODE in K.

Weak convergence

The assumptions below can be found in Section 10.2 and 10.4 of [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]:

Assumption 4.9.

1. (A2.0) Either

) is uniformly integrable for some ρ > 0.

3. (A2.2) Let (Z n (•)) be the sequence of continuous time processes defined by (4.2). Then (Z n (•)) converges weakly to z * . 4. (A2.4) g is continuously differentiable. 5. (A2.6 and A4.4) The Jacobian matrix A ′ of g at z * is a Hurwitz matrix. If ϵ n = 1 n , then A ′ + I 2 is also a Hurwitz matrix. 6. (A2.7) For some p > 0 and ρ > 0, sup

in probability, with Σ 1 some non-negative definite matrix.

(A4.5)

• Case ϵ n = 1 n : let P ′ be the unique symmetric and positive definite solution of the Lyapunov equation (A ′ ) ⊺ P ′ + P ′ A ′ = -I, and λ the largest positive number such that I ≥ λP ′ . Then λ > 1.

• Case µ n = o(ϵ n ): Let m(t) be the smallest integer n such that t n ≤ t < t n+1 , with m(t) = 0 for t ≤ 0. Then for all T > 0, lim inf n min m(tn-T )≤i≤n ϵn ϵ i = 1. Theorem 4.10 (Theorem 10.2.1 in Kushner and Yin (2003) (partial)). Assume that the assumptions 4.7 and 4.9 hold, and let U n (•) be the process defined by (4.3).

If µ n = o(ϵ n ), then the sequence of processes (U n (•)) converges weakly in the Skorohod space D n (R + ) to the stationary process U defined for all t ≥ 0 by:

e A ′ (t-s) dW s .

(4.4)

If ϵ n = 1 n , then the sequence of processes (U n (•)) converges weakly to the stationary process U defined for all t ≥ 0 by:

2 )(t-s) dW s , (4.5)

with in both cases W : R → R n a Wiener process of covariance matrix Σ 1 .

Proof. By Assumption 4.7 and 4.9, the assumption (A2.0)-(A2.7) p.329 of Theorem 10.2.1 are straightforward, at the exception of assumption (A2.3). The latter condition is also satisfied as a consequence of Theorem 10.4.1 p.341 in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] since the assumptions (A4.1)-(A4.5) of this Theorem are also satisfied by Assumption 4.7 (1), Theorem 4.8, and Assumption 4.9 ( 5) and ( 7).