Peggy Cenac
email: peggy.cenac@u-bourgogne.fr

Véronique Maume-Deschamps
email: veronique.maume@univ-lyon1.fr

Clémentine Prieur
email: clementine.prieur@imag.fr

P C Énac

SOME MULTIVARIATE RISK INDICATORS: MINIMIZATION BY USING A KIEFER-WOLFOWITZ APPROACH TO THE MIRROR STOCHASTIC ALGORITHM

published or not. The documents may come

Some multivariate risk indicators: minimization by using

a Kiefer-Wolfowitz approach to the mirror stochastic algorithm

Introduction

The new regulation rules for insurance industry, namely Solvency 2 in Europe, lead companies to adjust their solvency margins to the underlying risks. Once the overall company capital requirement has been computed, it must be splited into solvency capitals for each line of business. In other words, given an initial capital u we assume that it is allocated to each line of business: u k is the initial capital of the k th line of business, then u 1 + • • • + u d = u. We aim at optimizing the capital allocation with respect to some risk indicators. In this context u is fixed and we search for an optimal choice of the u k 's. The ruin probability is a quite standard risk indicator and has been widely studied in dimension 1 (see [START_REF] Gerber | An introduction to Mathematical Risk Theory[END_REF] for a review). The ruin probability is related to the much used Value at Risk (VaR) which is simply a quantile and do not take into account the amount or severity of ruin, moreover it is not coherent in Artzner's et al. sense (see [START_REF] Artzner | Heath Coherent Measures of Risk[END_REF]). Risk measures incounting severity of ruin or some penalties functions have also been studied in dimension 1 (see for example [START_REF] Dufresne | The surpluses immediately before and at ruin, and the amount of the claim causing ruin[END_REF][START_REF] Picard | On some measures of the severity of ruin in the classical Poisson model[END_REF]) and they are related to the Tail Value at Risk (TVaR) or Tail conditional Expectation (see [START_REF] Artzner | Heath Coherent Measures of Risk[END_REF]). Generalizations of these notions in higher dimension are not straightforward. Several attempts exist. The simplest way to construct risk measures or risk indicators for vectors could be to consider risk measures or indicators associated to the sum of the coordinates (this is a way to aggregate risks). But, in a ERM (Enterprise Risk Management) point of view, one should wonder about the risk induced by each business line. Generalizations of the ruin probability or the tail conditional expectation may be found in [START_REF] Cai | Dependence properties and bounds for ruin probabilities in multivariate compound risk models[END_REF][START_REF] Hult | Heavy-tailed insurance protfolios: buffer capital and ruin probabilities[END_REF][START_REF] Di Bernadino | Plugin estimation of level sets in a non-compact setting with applications in multivariate risk theory[END_REF].

Key words and phrases. Multivariate risk processes, risk indicators, stochastic algorithms, reserve allocation.

We would like to thank the referees for their careful reading of our manuscript and their comments. We are grateful to Anatoli Juditsky for his valuable advices. This work has been partially supported by the French Research National Agency (ANR) under the reference ANR-08-BLAN-0314-01.

In a multidimensional framework, risk indicators have been introduced by Loisel [START_REF] Loisel | Differentiation of some functionals of risk processes and optimal reserve allocation[END_REF], which do not take into account the dependence structure or the ruin amount. We introduce some new risk indicators, both involving dependence and local insolvency severity. Considering that the main risk drivers for the overall company have been identified and that the global solvency capital requirement has been computed, they reveal the marginal solvency capitals for each line of business. A way to avoid as far as possible that some lines of business become insolvent too often could be to minimize these risk indicators, under a fixed total capital constraint. This might be achieved if some capital fungibility between lines of business or between entities is possible. One possible way to define optimality of the global reserve allocation is to minimize the expected sum of the penalties that each line of business would have to pay due to its temporary potential insolvency.

We consider a vectorial risk process X i = (X 1 i , . . . , X d i ) t , where X k i corresponds to the gains of the k th business line during the i th period. That is, X k i = G k i -L k i where G k i denotes the incomes and L k i denotes the losses. We are interested in the cumulative gain:

Y k i = i p=1 X k p .
We shall say that the k th line of business makes default during the n first periods if there exists i = 1, . . . , n such that Y k i + u k < 0, where we recall that u k is the initial capital of the k th line of business. The ruin probability in a multivariate setting may be defined as the probability that one of the d lines of business makes default during the n first periods that is: R(u 1 , . . . , u d ) = P(∃k = 1, . . . , d ∃i = 1, . . . , n / Y k i + u k < 0). Other definitions of ruin probability and some bounds on them are provided in [START_REF] Cai | Dependence properties and bounds for ruin probabilities in multivariate compound risk models[END_REF].

Let R k p denotes the capital of the k th line of business at time p, that is

R k p = u k + Y k p .
We may also be interested in the total cost of ruin :

A(u 1 , . . . , u d ) = - d k=1 E   n p=1 R k p 1 {R k p <0}  
The asymptotic (as u → ∞) properties of A have been studied in [START_REF] Biard | Veraverbeke Asymptotic behavior of the finitetime expected time-integrated negative part of some risk processes and optimal reserve allocation[END_REF]. More precisely, in the case d = 1, they study the asymptotics of A(u). In higher dimension, they provide the asymptotic behavior of the vector u ∈ R d realizing the minimum of A under the constraint

u 1 + • • • + u d = u, u i ≥ 0, i = 1, . . . , d.
The risk indicator A does not take into account the dependence structure between the lines of business. A multivariate risk indicator that takes into account the dependence structure may be :

B(u 1 , . . . , u d ) = d k=1 E   n p=1 1 {R k p <0} 1 { d j=1 R j p >0}   .
This risk indicator is a measure of the time to ruin, it has been introduced in [START_REF] Loisel | Differentiation of some functionals of risk processes and optimal reserve allocation[END_REF]. It does not take into account the severity (or the cost) of the insolvency. We introduce the following risk indicators, related to the ruin cost and dependending on penalty functions g k . Given convex functions g k : R → R, k = 1, . . . , d satisfying g k (x) ≥ 0 for x ≤ 0, and

E(|g k (R k p )|) < ∞, we consider the risk indicator (1.1) I(u 1 , . . . , u d ) = d k=1 E   n p=1 g k (R k p )1 {R k p <0} 1 { d j=1 R j p >0}   .
The function g k represents the cost that the k th business branch has to pay when it becomes insolvent while the whole company is solvent. The risk indicator I is constructed with the idea that the solvency requirement for the whole company has been determined (namely u is determined). Minimizing I allows to control a partial insolvency, that is, insolvency measured for some branches while the whole company is solvent.

A natural choice for the cost could be g k (x) = -x. The following picture represents the risk indicator I (light area) for g k (x) = -x and d = 2.

time 1 2 3 Local insolvency Global insolvency u1 R 2 1 u2 R 1 2 R 1 1 + R 2 1 u = u1 + u2 d k=1 R k p R 1 p R 2 p
We aim at finding a minimum under constraint on the initial capital u. Formally, we are looking for u

∈ (R + ) d such that (1.2) I(u ) = inf v 1 +...+v d =u I(v), v ∈ (R + ) d .
Unless for very specific examples, no explicit solution is available. We propose to solve this problem by using stochastic algorithms (see [START_REF] Duflo | Random Iterative Models[END_REF][START_REF] Bartoli | Del Moral Simulation & algorithmes stochastiques[END_REF][START_REF] Hj Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] for a review). The classic algorithm introduced by Robbins and Monro [START_REF]Robbins & Monro A stochastic approximation method[END_REF] consists in applying a gradient descent method with a noisy gradient. In case such a noisy gradient is not observable, Kiefer and Wolfowitz (see [START_REF] Kiefer | Wolfowitz Stochastic estimation of the maximum of a regression function[END_REF]) developed a finite difference version of the Robbins-Monro method. The inconvenient of such algorithms is their very high sensitivity to the choice of the stepsize, which are very difficult to calibrate in practice. Moreover, we are here in the context of constraint optimization problem. To handle the constraints, we first turned to the Lagrange multipliers approach. It was not tractable in practice in our setting, it did not converge on our simulations. Improvements of Robbins-Monro and Kiefer-Wolfowitz algorithms have been developed more recently using an averaging tool (see for instance [START_REF] Polyak | Juditsky Acceleration of stochastic approximation by averaging SIAM[END_REF][START_REF] Pelletier | Asymptotic almost sure efficiency of averaged stochastic algorithms[END_REF]) or with mirror versions of the stochastic descent algorithms. Deterministic form of the mirror descent algorithm is introduced by Nemirovski and Yudin [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF].

The stochastic one has been proposed by Nesterov [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] and used in [START_REF] Juditsky | Recursive aggregation of estimators via the mirror descent algorithm with average[END_REF][START_REF] Juditsky | Learning by mirror aver aging[END_REF] to solve non-smooth stochastic convex optimization problems. Recently new stochastic approximation algorithms, with stronger convergence properties, based on Nesterov's method have been proposed (see [START_REF] Ghdimi | Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization[END_REF][START_REF] Lan | Efficient methods for stochastic composite optimization[END_REF][START_REF] Lan | An optimal method for stochastic composite optimization[END_REF][START_REF] Nemirovski | Robust stochastic approximation to stochastic programming[END_REF]). All these algorithms rely on the existence of a stochastic oracle, that is a mechanism to generate noisy versions of the gradient. In our problem, we do not have access to a noisy gradient (or to a noisy subgradient). We are thus led to use approximations of ∇I in stochastic algorithms. That is the reason why we have chosen to follow a Kiefer-Wolfowitz approach, in order to take into account an estimation of the gradient online. Our proof for the convergence of the algorithm is inspired by the work of Tauvel [START_REF] Tauvel | Optimisation stochastique à grande échelle[END_REF].

We construct a stochastic algorithm to estimate u , the optimal allocation. Under a moment condition (order > 2) we prove the almost sure convergence of the estimator. Notice that our approach do not require any hypothesis, apart from a moment assumption, on the distribution of X i and the existence of a density for some parts of our work. Also there may be temporal dependence on a period of length n. Nevertheless, in order to achieve the algorithm, we shall need N independent copies of the distribution of (X 1 , . . . , X n ) ∈ (R d ) n , that is of the gains over a period of length n.

The paper is organized as follows: in section 2 we prove some convexity properties on I. In section 3 we make explicit the algorithm. In section 4 we prove the convergence results. In section 5 we show how to apply these results to our risk indicator and provide some simulations. While it is well known that standard Robbins-Monro and Kiefer-Wolfowitz algorithms are very sensitive to the choice of the stepsize, it is remarkable that our Kiefer-Wolfowitz version of the mirror algorithm is very stable. Namely, it seems, from our simulations, that any choice of the sequences in the theoretical convergence domain is convenient. So that our algorithm is very tractable while more standard ones are very hard to use in practice. For our purpose, more classical algorithms were not tractable, either because of the instability of the algorithm or because we do not have access to the gradient or even to an oracle of this gradient. Also, we have performed high dimensional simulations (10 and 40) and still the results are very satisfactory.

Convexity property

Let us first prove that, under suitable assumptions on g k , the risk indicator I is convex.

Let v (-k) = j =k v j and S p = d k=1 Y k p .
Proposition 2.1. Assume that the functions g k are convex, that g k : R → R with g k (0) = 0, g k (x) ≥ 0 for x < 0, and for all k = 1, . . . , d and p = 1, . . . , n, g k (R k p ) admits an order 1 moment. Then the risk indicator I is convex on the convex set

U u = {(v 1 , . . . , v d ) ∈ (R + ) d / v 1 + • • • + v d = u}. If moreover, for some k = 1, . . . , d, g k is strictly convex, g k (x) > 0 for x < 0 and the joint distribution of (Y k p , S p ) contains ] -u, +∞[×] -∞, u[ in its support then I admits a unique minimum in U u .
Proof. We rewrite I as

I(u) = I(u 1 , . . . , u d ) = d k=1 E   n p=1 g k (Y k p + u k )1 {Y k p <-u k } 1 {Sp>-u}   . Since g k (0) = 0, (u 1 , . . . , u d ) → g k (Y k p + u k )1 {Y k p <-u k } 1 {Sp>-u}
is a convex function on U u . So that I is a convex function as the sum of expected values of convex functions. The unicity of the minimum follows from the fact that I is strictly convex, due to the strict convexity of the g k 's.

Remark 1. If, in addition to the hypothesis of Proposition 2.1, the functions g k are differentiable and such that for all k = 1, . . . , d and p = 1, . . . , n, g k (Y k p + u k ) admits an order 1 moment and for all i = 1, . . . , d, (Y i p , S p ) has a joint density distribution denoted by f Y i p ,Sp , then the risk indicator I is differentiable. Indeed, it writes:

I(u) = I(u 1 , . . . , u d ) = d k=1 E   n p=1 g k (Y k p + u k )1 {Y k p <-u k } 1 {Sp>-u}   = d k=1 n p=1 -u k -∞ ∞ -u g k (y + u k )f Y k p ,Sp (y, s)dyds.
Hence, the assumptions on g k imply that I is differentiable. Moreover, we have:

(∇I(v)) i = d k=1 n p=1 -v k -∞ g k (y + v k )f Y k p ,Sp (y, -u)dy + n p=1 E(g i (Y i p + v i )1 {Y i p <-v i } 1 {Sp>-u} ) - n p=1 g i (0) +∞ -u f Y i p ,Sp (-v i , s)ds = d k=1 n p=1 -v k -∞ g k (y + v k )f Y k p ,Sp (y, -u)dy + n p=1 E(g i (Y i p + v i )1 {Y i p <-v i } 1 {Sp>-u} ).
Thus ∇I(v) is a sum of two terms. The second term appears as an expectation of an observable function which could be an oracle for a stochastic algorithm design. For the first part of the sum, such an oracle does not seem simply reachable, this is why we turned to a Kiefer-Wolfowitz approach.

The fact that I is a convex function is used to prove the convergence of our stochastic algorithm. The convexity assumption is also standard for classical stochastic algorithms (as Robbins-Monro or Kiefer-Wolfowitz) and also in more recent extensions ([9, 15, 16, 20]).

Estimation by use of stochastic algorithms

In this section we present our algorithm which is a Kiefer-Wolfowitz version of the mirror algorithm. In section 4 we shall prove its convergence under general assumptions (see Assumption 1) that are satisfied in particular for our specific risk indicator I.

3.1.

Auxiliary function for the mirror algorithm. The stochastic mirror algorithm requires the use of an auxiliary function that will be used to push the trajectory into the set of constraints. We consider that R d is endowed with the L 1 norm and the dual space (R d ) is endowed with the dual norm :

x = d i=1 |x i |, ξ = sup i=1,...,d |ξ i |.
Let C be a compact convex subset of R d . Recall that a convex function δ : C -→ R is strongly convex with parameter α > 0 with respect to the norm if for any λ ∈ [0, 1] and any x, y belonging to C,

δ(λx + (1 -λ)y) ≤ λδ(x) + (1 -λ)δ(y) - α 2 λ(1 -λ) x -y 2 .
The mirror algorithm uses an α-convex function δ differentiable on a point x 0 ∈ C to define the auxiliary function V . Following Definition A.3 of [START_REF] Tauvel | Optimisation stochastique à grande échelle[END_REF], we consider:

V (x) = δ(x) -δ(x 0 )-< ∇δ(x 0 ), x -x 0 > .
It is then easy to see that the function V is also α-convex. Setting a positive β, W β denotes the Legendre-Fenchel transform of βV :

W β (ξ) = sup x∈C { ξ, x -βV (x)} .
Due to the properties of Legendre-Fenchel transform of α-convex functions, we know that W β is continuously differentiable and that its gradient takes its values in C. Moreover, the gradient of W β has a Lipschitz constant equal to (αβ) -1 . We consider a function Ψ which is an approximation of the gradient of I.

The mirror algorithm uses two positive sequences (β n ) n and (γ n ) n and a sequence of i.i.d. random vectors (Y n ) n in the following way (the choices for these sequences will be made explicit a bit later).

Algorithm 1

Initialization :

ξ 0 = 0 ∈ (R d ) χ 0 ∈ C Update : for i = 1, . . . , N do ξ i = ξ i-1 -γ i Ψ(χ i-1 , Y i ) χ i = ∇W β i (ξ i ) Output : S N = N i=1 γ i χ i-1 N i=1 γ i
Under assumptions on the function Ψ, we shall prove that S N goes a.s. to the unique minimum x of I.

Approximate gradient.

In our case, we do not have access directly to the gradient of I, ∇I, but we may approximate it online.

Recall that

I(u 1 , . . . , u d ) = d k=1 E   n p=1 g k (R k p )1 {R k p <0} 1 { d k=1 R k p >0}   .
Let us recall that d denotes the number of lines of business, and n the number of periods of interest. The problem is to solve (1.2). We will consider an approximation Ψ of ∇I with the following decomposition Ψ = ∇I + η + r, where η is a martingale difference and r is negligible with respect to ∇I and η. Arguing as for Kiefer-Wolfowitz, we first notice that

I(u 1 , . . . , u d ) = E (I(u 1 , . . . , u d , Y))
where

Y =   Y 1 1 • • • Y 1 n • • • • • • • • • Y d 1 • • • Y d n   and
I(u 1 , . . . , u d , y) = d k=1 n p=1 g k (y k p + u k )1 {y k p +u k <0} 1 { d k=1 y k p +u k >0} .
Let us denote a realization of Y by

y =   y 1 1 • • • y 1 n • • • • • • • • • y d 1 • • • y d n   .
For sake of shortness, we introduce the following notations :

I k c + i , Y = I(χ 1 i-1 , . . . , χ k-1 i-1 , χ k i-1 + c i , χ k+1 i-1 , . . . , χ d i-1 , Y), I k c - i , Y = I(χ 1 i-1 , . . . , χ k-1 i-1 , χ k i-1 -c i , χ k+1 i-1 , . . . , χ d i-1 , Y), I k c + i = I(χ 1 i-1 , . . . , χ k-1 i-1 , χ k i-1 + c i , χ k+1 i-1 , . . . , χ d i-1 ), I k c - i = I(χ 1 i-1 , . . . , χ k-1 i-1 , χ k i-1 -c i , χ k+1 i-1 , . . . , χ d i-1 ). We then consider D c i I (resp. D c i I) the random vector whose k th coordi- nate D k c i I(u 1 , . . . , u d , Y) (resp. D k c i I(u 1 , . . . , u d , Y)) is defined by I k c + i , Y -I k c - i , Y 2c i resp. I k c + i , Y -I k c - i , Y 2c i . Let Ψ c i (χ i-1 , Y i ) = D c i I(χ i-1 , Y i ).
Considering independent copies of the random matrix Y, denoted by Y 1 , . . . , Y N , we perform Algorithm 2 below :

Algorithm 2

Initialization :

ξ 0 = 0 ∈ (R m ) χ 0 ∈ C Update : for i = 1, . . . , N do ξ i = ξ i-1 -γ i Ψ c i (χ i-1 , Y i ) χ i = ∇W β i (ξ i ) Output : S N = N i=1 γ i χ i-1 N i=1 γ i
Then defining

η c i (χ i-1 , Y i ) = D c i I(χ i-1 , Y i ) -D c i I(χ i-1 ), r c i (χ i-1 ) = D c i I(χ i-1 ) -∇I(χ i-1 ),
it enhances the decomposition

(3.1) Ψ c i (χ i-1 , Y i ) = ∇I(χ i-1 ) + η c i (χ i-1 , Y i ) + r c i (χ i-1 ).
In the next section, we give conditions on (γ i ) i , (β i ) i and (c i ) i for the convergence of the algorithm.

Convergence of the algorithm

Let us consider C a compact convex subset of R d . We make also the following assumptions on I :

Assumption 1.
(1)

I is a convex function from R d to R, (2)
I is of C 2 class, (3)
I admits a unique minimum x on C, Moreover, we have for any x in R d and any c > 0 E(η c (x, Y)) = 0 and r c (x) converges to zero as c goes to 0. We perform Algorithm 2 with χ 0 being any random vector taking values in C, and F i being the σ-field generated by χ 0 , Y 1 , . . . , Y i . Then ξ i is F i measurable and Y i is independent of F i-1 . We get therefore

E(η c i (χ i-1 , Y i )|F i-1 ) = 0.
We shall need the following additional assumption on the function I.

Assumption 2. There exists some nonnegative real number σ such that for all vector (v 1 , . . . , v d ) ∈ R d , the random variable I(x, Y) satisfies the moment condition ∀i

Var(I(v 1 , . . . , v d , Y i )|F i-1 ) ≤ σ 2 .
Remark 2. Assumption 2 yields the existence of a real number σ > 0 such that for any i = 1, . . . , N , (4.1)

E η c i (χ i-1 , Y i ) 2 * |F i-1 ≤ d σ 2 c 2 i .
Indeed, one has

E η c i (χ i-1 , Y i ) 2 * |F i-1 ≤ d k=1 Var I k (c + i , Yi)|F i-1 4c 2 i + Var I k (c - i , Y i )|F i-1 4c 2 i - 1 2c 2 i E I k (c + i , Y i ) -I(c + i ) I k (c - i , Y i ) -I(c - i ) .
From Assumption 2 together with the Cauchy-Schwartz inequality, one gets

E η c i (χ i-1 , Y i ) 2 * |F i-1 ≤ d σ 2 c 2 i ,
which concludes the proof of (4.1). Condition (4.1) is sufficient to obtain the L 1 -convergence. In order to prove the almost sure convergence, we shall need a higher order moment condition.

Assumption 3. There exists p > 2 such that almost surely

sup i>0 E (|I(v 1 , . . . , v d , Y)| p |F i-1 ) < ∞.
Remark 3. Assumption 3 yields the existence of a real p > 2 such that almost surely

sup i>0 E( c i η c i (χ i-1 , Y i ) p |F i-1 ) < ∞.
Our main result is the following.

Theorem 4.1. Let (β i ) i≥0 , (γ i ) i≥0 and (c i ) i≥0 be sequences in R * + N . Assume moreover that (β i ) i≥0 is non decreasing and that the following assumptions are satisfied :

(i) β N / N i=1 γ i -----→ N →+∞ 0, (ii) N i=1 γ i c i / N i=1 γ i -----→ N →+∞ 0, (iii) N i=1 γ 2 i c 2 i β i-1 / N i=1 γ i -----→ N →+∞ 0, (iv) +∞ i=1 γ i c i 2 < ∞.
Then, provided that Assumptions 1 and 2 are satisfied, we have :

S N L 1 -→ x .
If moreover, Assumption 3 is satisfied, we have :

S N a.s. -→ x .
Remark 4. Let 0 < c < a -1 2 < 1 2 . One can choose for (β i ) i≥0 the constant sequence equal to 1, for (c i ) i≥0 the sequence (i -c ) i≥0 and for (γ i ) i≥0 the sequence (i -a ) i≥0 .

Remark 5. For the "classical" Kiefer-Wolfowitz algorithm, the assumptions on the gains are particular cases of our assumption with a constant sequence (β i ) i≥0 :

γ n = ∞, γ n c n < ∞, and γ 2 n c -2 n < ∞.
The proof of Theorem 4.1 uses an estimation of ε N = I(S N ) -I(x ) inspired from [START_REF] Tauvel | Optimisation stochastique à grande échelle[END_REF].

Lemma 4.2. Assume that (β i ) i∈N is a non decreasing sequence, then 0 ≤ N i=1 γ i ε N ≤ β N V (x ) - N i=1 γ i η c i (χ i-1 , Y i ), χ i-1 -x - N i=1 γ i r c i (χ i-1 ), χ i-1 -x + N i=1 γ 2 i 2αβ i-1 Ψ c i (χ i-1 , Y i ) 2 .
Proof. We remark that because I (x ) = inf x∈C I(x) then ε N ≥ 0 and for any i = 1, . . . , N , I(χ i-1 ) -I(x ) ≥ 0. Because I is a convex and differentiable function,

N i=1 γ i ε N ≤ N i=1 γ i (I(χ i-1 ) -I(x )) ≤ N i=1 γ i ∇I(χ i-1 ), χ i-1 -x . Thus, N i=1 γ i ε N ≤ N i=1 γ i Ψ c i (χ i-1 , Y i ), χ i-1 -x - N i=1 γ i η c i (χ i-1 , Y i ), χ i-1 -x - N i=1 γ i r c i (χ i-1 ), χ i-1 -x = ξ N , x + N i=1 γ i Ψ c i (χ i-1 , Y i ), χ i-1 - N i=1 γ i η c i (χ i-1 , Y i ), χ i-1 -x - N i=1 γ i r c i (χ i-1 ), χ i-1 -x
Now, because the sequence (β i ) i∈N is non decreasing and using the convexity and the differentiability of W β i-1 (as the Legendre transform of a strictly convex function), Lemma A.1. of Tauvel [START_REF] Tauvel | Optimisation stochastique à grande échelle[END_REF] yields :

W β i (ξ i ) -W β i-1 (ξ i-1 ) ≤ W β i-1 (ξ i ) -W β i-1 (ξ i-1 ) ≤ ∇W β i-1 (ξ i-1 ), ξ i -ξ i-1 + γ 2 i 2αβ i-1 Ψ c i (χ i-1 , Y i ) 2 * = -γ i χ i-1 , Ψ c i (χ i-1 , Y i ) + γ 2 i 2αβ i-1 Ψ c i (χ i-1 , Y i ) 2 * .
As a consequence, we get:

N i=1 γ i Ψ c i (χ i-1 , Y i ), χ i-1 ≤ W β 0 (ξ 0 ) -W β N (ξ N ) + N i=1 γ 2 i 2αβ i-1 Ψ c i (χ i-1 , Y i ) 2 * = -W β N (ξ N ) + N i=1 γ 2 i 2αβ i-1 Ψ c i (χ i-1 , Y i ) 2 * .
The result finally follows from:

x , ξ N -W β N (ξ N ) ≤ β N V (x ).
Remark 6. With our assumptions, we have that

E(< η c i (χ i-1 , Y i ), χ i-1 -x >) = 0, indeed, (< η c i (χ i-1 , Y i ), χ i-1 -x >) i is a martingale difference sequence.
The following lemma may be deduced from Taylor inequality, recalling that χ i ∈ C ∀i. Lemma 4.3. There exists κ > 0 such that, for any i,

r c i (χ i-1 ) ≤ κc i .
Proof. It follows from the fact that the Hessian matrix of I on C is bounded.

Remark 7. As a consequence of Lemmas 4.2 and 4.3, we have that (recall that ε N ≥ 0)

0 ≤ E(|ε N |) = E(ε N ) ≤ β N V (x ) N i=1 γ i + 2κD N i=1 γ i c i N i=1 γ i + N i=1 Bγ 2 i 2c 2 i αβ i-1 N i=1 γ i , where D is the diameter of C and B c 2 i denotes a bound of E Ψ c i (χ i-1 , Y i ) 2 . Indeed, because of Decomposition (3.1) Ψ c i 2 ≤ ∇I 2 + r c i 2 + η c i 2 (4.2) + 2 ( ∇I + r c i η c i + ∇I . r c i ).
Recall that ∇I is bounded (because C is compact and ∇I is continuous), also following Lemma 4.3, r c i ≤ κc i . Finally, using (4.1),

E( η c i 2 ) ≤ dσ 2 c 2 i .
So that, sufficient conditions to have the L 1 -convergence of (ε N ) are Assumptions (i), (ii) and (iii) of Theorem 4.1.

We shall apply some standard results for martingales. Let us fix some notations : Definition 1. Let M N be a square integrable martingale written as

M N = N i=1 e i
with (e i ) i a martingale difference sequence. Let

M N = N i=1 E(e 2 i |F i-1 ),
with F i-1 = σ (e 1 , . . . , e i-1 ).

The following results are well known.

Theorem (distribution of large numbers for martingales, see for instance [START_REF] Neveu | Martingales à temps discret[END_REF]). Let M N be a square integrable martingale if

M ∞ = lim N →∞ M N < ∞ a.s.
then (M N ) N ∈N converges a.s. to a square integrable M ∞ .

Theorem (Chow Lemma, see for instance [START_REF] Duflo | Random Iterative Models[END_REF] p.22). Suppose (a N ) N ∈N is a bounded sequence of positive random variables, suppose that 1 < p ≤ 2. For

N ∈ N, let A N = 1 + N k=0 a k and A ∞ = lim N →∞ A N . Suppose that (Z N ) N ∈N is
a positive sequence of random variables adapted to F N and K is a constant such that

E(Z N +1 |F N ) ≤ K and sup N E(Z p N +1 |F N ) < ∞
then we have the following properties almost surely :

(4.3) on {A ∞ < ∞} ∞ k=1 a k Z k+1 converges (4.4) on {A ∞ = ∞} A -1 N N k=1 a k Z k+1 ≤ K.
Proof of Theorem 4.1. We have proved the convergence in L 1 of (ε N ) N to 0 (see Remark 7). The convergence of S N to x follows: I(S N ) converges in L 1 to I(x ), so it converges also in probability to I(x ). The unicity of x (as a minimum of I) and the continuity of I gives that for any δ > 0, there exists > 0 such that:

|s -x | > δ =⇒ |I(s) -I(x )| > .
Now, since S N and x belong to the compact set C, there exists K > 0 such that for any N ∈ N and any δ,

E|S N -x | ≤ δ + |S N -x |P (|S N -x | > δ) ≤ δ + KP (|S N -x | > δ) ≤ δ + KP (|I(S N ) -I(x )| > )
Thus, for any δ, lim sup

N →∞ E|S N -x | ≤ δ (using the convergence in probability)
and finally we get the convergence in L 1 of S N to x . Let us prove the almost sure convergence. We apply the distribution of large numbers for martingales to

M N = N i=1 γ i < η c i (χ i-1 , Y i ), χ i-1 -x >,
M N is a martingale. Under Assumption (iv), M N is bounded and hence converges because it is increasing, thus M N converges almost surely. Using Lemma 4.2 we get:

0 ≤ ε N ≤ β N V (x ) N i=1 γ i - M N N i=1 γ i + 2Dκ N i=1 c i γ i N i=1 γ i (4.5) + N i=1 γ 2 i 2αβ i-1 Ψ c i (χ i-1 , Y i ) 2 N i=1 γ i . (4.6)
The three terms in (4.5) go to zero a.s. by using the convergence of M N , (i) and (ii), remark that (i) implies that γ i = +∞. In order to prove that the term in (4.6) converges to zero a.s., we use the decomposition (4.2) and we apply the Chow Lemma twice to

Z N = c 2 N -1 η c N -1 2 and a N = γ 2 N β N -1 c 2 N (resp. Z N = c N η c N and a N = γ 2 N β N -1 c N ).
Remark that in that case, we apply the Chow Lemma to constant random variables (a n ) n∈N . With our hypothesis on the sequences (γ n ) n∈N , (c n ) n∈N and (β n ) n∈N , the two considered sequences (a n ) n∈N are summable. Theorem 4.1 is proved. Indeed, we have proved that

I(S N ) a.s. -→ I(x ).
Now, one gets S N a.s.

-→ x * because S N belongs to the compact set C. If s is an accumulation point then because I(S N ) -I(x ) goes to zero, s is also a minimum of I. So that s = x by uniqueness of the minimum. Remark 8. If some higher order moment assumption holds we could obtain an estimation of the rate of convergence, by using Markov inequality and the Chow Lemma.

Application to risk indicators and simulations

We go back to our initial problem: we aim at estimating the optimal allocation of the function I (defined by (1.1)) on the simplex

U u = {(v 1 , . . . , v d ) ∈ R d + , v 1 + • • • + v d = u}.
In order to run our mirror algorithm, we have to find an auxiliary function V , whose Legendre-Fenchel transform is computable. A natural strongly convex function on U u is the entropy function :

δ u (x) = d k=1 x k u ln x k u ,
it is 1 u 2 -convex and admits a unique minimum x 0 = ( u d , . . . , u d ). As in [START_REF] Tauvel | Optimisation stochastique à grande échelle[END_REF], we consider the function

V (x) = δ u (x) -δ u (x 0 ) -∇δ u (x 0 ), x -x 0 = ln d + d k=1 x k u ln x k u .
This function is also 1 u 2 -convex and its Legendre-Fenchel transform is easily computable :

∇W β (ξ) = β ln 1 d d k=1 exp[ξ k u β ] .
In Proposition 2.1 above, we have proved that our risk indicator I is convex and admits a unique minimum provided that g k (x) > 0 for x < 0 and that the joint distribution of (

Y k p , S p ) contains ]-u, +∞[×]-∞, u[ in its support.
The present section is devoted to simulations, whose aim is more to discuss the performances of our algorithm on simple test cases than to investigate more complex models. A more complete study on realistic models (heavy tailed distributions, temporal and spatial dependencies) will be the purpose of a future work, which should be written as a useful tool for risk managers. Here we focus on what happens when dealing with Gaussian distributions, and also when handling the so-called common shocks model. In all the examples below, the g k 's are chosen equal to -Id : g k (x) = -x for all k = 1, . . . , d. Moreover the cumulative gain is considered only on one period. In that case, n = 1 thus Y = X ∈ R d . In Section 5.1, the components of X are independent. We first consider that they also are identically distributed (see Subsection 5.1.1), then one knows that the minimum is achieved for u * = ( u d , . . . , u d ). We can therefore study the performances of our algorithm for various tuning of the parameters, (γ n ) n∈N and (c n ) n∈N , on this simple test case. The sequence (β n ) n∈N is fixed equal to 1 in the following. Various choices for (γ n ) n∈N and (c n ) n∈N are tested, provided that they satisfy conditions (i)-(iv) of Theorem 4.1. It is well known that for standard algorithms such as Robbins-Monro or Kiefer-Wolfowitz it is hard to calibrate these sequences to get the convergence of the algorithm in practice. In the following, the sequences (γ n ) n∈N and (c n ) n∈N are chosen as follows :

• γ n = 1 (n + 1) α , • c n = 1 (n + 1) δ , with 0 < δ < α -1/2 < 1/2.
On our test cases, it seems that choosing α too close to one can deteriorate the mean squared error, but the overall feeling is that our algorithm is not too sensitive to the parameters (see Subsection 5.1.1). Thus from Subsection 5.1.2, we fixed α = 3 4 + 1 10 and δ = 1 4 . In Subsection 5.1.2 we still consider independent normal coordinates for X but the marginals may have different means or variances. In Section 5.2, but some dependencies on the coordinates of X are allowed.

In section 5.3, some higher dimensional problems are considered (dimensions 10 and 40).

Also, we have chosen to perform the simulations with the initial capital u = d (except when another value for u is specified). The initialization of the algorithm (χ 0 ) is done at random uniformly in the simplex U u . 5.1. Independent models. The aim of this section is to provide some benchmark for the use of the algorithm. So we consider : d = 2, n = 1, the coordinates X 1 and X 2 are independent. At first, we assume that X 1 and X 2 have the same normal distribution, then we consider different normal distributions.

5.1.1. Same normal distribution. We have chosen an independent normal distribution with mean equal to 0.3 and standard deviation equal to 1. For N = 10000 independent simulations we obtain : S N = (0.996; 1.004) and the following graphs. Remark that, as expected, it seems that the minimum of our risk indicator is reached for u

* 1 = u * 2 = u 2 = 1.
In order to illustrate the convergence, we have plotted for each of the two coordinates the trajectory χ i (black line) and S i (gray line); we have plotted the whole trajectories (i = 1, . . . , 10000) and a zoom on the end of the trajectories (i = 2000, . . . , 10000). We see on these graphs the importance of the averaging step which provides a solution which is much more smooth. We conclude this subsection by a study of k = 30 simulations of length N = 1000, with various sequences (γ n ) n∈N and (c n ) n∈N . These simulations allow to see the influence of these sequences on the results of the algorithms. The performance is measured with the mean squared error :

mse = 1 k k j=1 u -u j 2 ,
where u j is the estimated value of the vector u on the j th simulation. The table below gives for each of the two coordinates, the mean of the estimation ( u 1 , u 2 ) of the minimum (u 1 , u 2 ), the standard error and the mean squared error mse for sequences (γ n ) n∈N and (c n ) n∈N :

γ n = 1 (n+1) α and c n = 1 (n+1)
δ for various values of α and δ satisfying the hypothesis of Theorem 4.1. The parameter δ is used for finite differences and it seems that it is better to choose it not too close to 1, as it deteriorates the mean squared error. The parameter α (used to define the descent step (γ n ) n∈N ) does not seem to have a great influence on our test case. We observe that there is a significant difference between the two coordinates. This is intuitive that in order to minimize the risk indicator, one should affect more initial capital to the branch whose expected income is the less (which is the more risky branch in that case). We achieve this section with independent normal distributions with different variance (and same mean). As above, we have performed 50 simulations of length N = 1000, with two independent normal distributions. For the first coordinate we consider N (0.3, 1) and for the second one N (0.3, 4). The table below gives for each of the two coordinates, the mean of the estimation ( u 1 , u 2 ) of the minimum (u 1 , u 2 ) and the standard error.

α = 1 α = 1 α = 3 4 + 1 10 α = 3 4 + 1 10 α = 3 4 α = 3 4 δ = 1 2 δ = 1 4 δ = 1 2 δ = 1 4 δ = 1 2 δ = 1 4 1 st coord
first coord. second coord. mean 0.787 1.213 standard error 0.067 0.067

In that case also, we observe that there is a significant difference between the two coordinates. Here, the more risky branch is the second one (which has a higher variance).

5.2. Models with spatial dependencies. We propose two simple models of spatial dependence. First of all we consider random vectors in R 3 with common shocks. A second example of spatial dependence is given by Gaussian vectors.

5.2.1.

Model with common shocks. Models with common shocks are very popular in insurance because they allow to model dependencies due to some common factor. We consider a dimension three vector: We have also performed simulations for a Gaussian vector in R 3 with covariance matrix Σ =   1 0 0 0 1 0.9 0 0.9 1   and expectation m = (0.3, 0.3, 0.3). As above, we have performed 50 simulations of length N = 1000, of the Gaussian vector X. The table below gives for each of the three coordinates, the mean of the estimation ( u 1 , u 2 , u 3 ) of the minimum (u 1 , u 2 , u 3 ) and the standard error. first coord. second coord. third coord. mean 0.785 0.604 0.612 standard error 0.045 0.03 0.028

X 1 = I × Z 1 + (1 -I) × W , X 2 = I × Z 2 + (1 -I) × W , X 3 = Z 3 with Z 1 , Z 2 , Z 3 independent normal distributions N (0.
As in the example with common shocks, the strategy is to put more capital on the entity constituted with the positively correlated branches.

5.3. Simulations in higher dimension. We have performed simulations in higher dimension, still considering one period for the cumulative gain (n = 1 ). First for Gaussian random vectors in dimension 10 then in dimension 40, in order to see if our algorithm gives good results in high dimension. In order to measure the performance of the algorithm, we have performed k = 30 simulations of length N = 1000, with independent normal Gaussian vectors of R 10 then of R 40 , each coordinate following a normal N (0.3, 1). We have chosen for the sequences γ n = 1 (n+1) and c n = 1 (n+1) .

5.3.1.

Results in dimension 10. The table below contains for each coordinate, the mean -over the 30 simulations -of the estimated value of u i and the standard-error. Then we have computed the mse. Remark that as before, since the coordinates of the vector are i.i.d. then the minimum of the indicator will be realized for u i = u 10 = 1 for i = 1, . . . , The mean squared error is 0.02. Even if the mean squared is greater (around 5 times greater as in dimension 2), it seems that the algorithm is still performant in dimension 10.

We have also considered in dimension 10, a Gaussian vector with mean 0.3 for all the coordinates and with the covariance matrix:

               
1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5

               
The table below contains for each coordinate, the mean -over the 30 simulations -of the estimated value of u i and the standard-error.

coordinate 0.811 0.034 Clearly, the algorithm recognizes two blocks in the random vector and puts more capital on the more risky (with higher variance) one.

5.3.2.

Results in dimension 40. We conclude this simulation session by a simple example in dimension 40. We have performed k = 30 simulations of length N = 1000, with independent normal Gaussian vectors each coordinate following a normal distribution N (0.3, 1). We have chosen for the sequences γ n = 1 (n+1) and c n = 1 (n+1) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 10 20 30 40 0.99 1.00 1.01 1.02

Estimation for an i.i.d. vector in dimension 40

Index estimated values

All these results show that our algorithm is tractable even for relative high dimensions. Of course the mean squared error increases with the dimension but even in dimension 40 we obtain a quite good precision for our estimator (it seem that the mse increases linearly with the dimension). Note that for all these simulations, we did not encounter any difficulty for fitting our various parameters. More realistic models, including temporal dependencies and heavy tailed distributions, have to be investigated from a theoretical and simulation point of view, in order to propose a useful tool for risk managers. This will be the aim of a forthcoming work.

Figure 1 .Figure 2 .

 12 Figure 1. Convergence of the algorithm, first coordinate

1 4 .

 14 The mean squared error is 0.09.The picture below shows the mean (over the 30 simulations) of the estimated values u 1 , . . . , u 40 ofu 1 = • • • = u 40 = u 40 = 1.

Figure 3 .

 3 Figure 3. Estimated values of u 1 , . . . , u 40 .

  The table below gives for each of the two coordinates, the mean of the estimation ( u 1 , u 2 ) of the minimum (u 1 , u 2 ).

	However, even with a naive fitting of the parameters we observe the conver-
	gence of the algorithm in practice.
	5.1.2. Different normal distributions. For this subsection, we have chosen to
	present the results of 50 simulations of length N = 1000, with two indepen-
	dent normal distributions. First of all, we consider for the first coordinate a
	normal distribution N (0.3, 1) and for the second one a normal distribution
	N (0.8, 1).
			first coord. second coord.
		mean		1.226	0.774
	standard error	0.051	0.051
	. mean	1.02	1.01		1.01	1.02	1.01	1.01
	1 st standard error 0.084	0.071	0.04	0.04	0.03	0.03
	2 nd coord. mean	0.98	0.99		0.99	0.98	0.99	0.99
	2 nd standard error 0.084	0.071	0.04	0.04	0.03	0.03
	mse	0.015	0.01		0.003	0.003	0.002	0.002

  As above, we have performed 50 simulations of length N = 1000, of the Gaussian vector X. The table below gives for each of the two coordinates, the mean of the estimation ( u 1 , u 2 ) of the minimum (u 1 , u 2 ) and the standard error.

	first coord. second coord. third coord.
	mean	0.61	0.61	0.78
	standard error	0.024	0.033	0.042
	The first two coordinates are positively correlated. If one fails then the
	probability that the other one fails is more important. Putting more capital
	for the entity constituted with these first two branches should reduce the
	insolvency risk of the branches.
	5.2.2. Gaussian vectors. We conclude this paragraph on vectorial depen-
	dence with an example of a Gaussian vector in dimension 2 and an example
	of a Gaussian vector in dimension 3. First, we consider the Gaussian vector
	with covariance matrix
		Σ =	1 0.8 0.8 1
	and expectation m = (0.3, 0.3). first coord. second coord.
	mean		0.999	1.001
	standard error	0.034	0.034
	It seems that the result is the same as in the non correlated case but with
	a smaller standard deviation.
	We have also performed the simulations for a Gaussian vector with expec-
	tation m = (0.3, 0.8) and the covariance matrix Σ. The table below gives
	for each of the two coordinates, the mean of the estimation ( u 1 , u 2 ) of the
	minimum (u 1 , u 2 ) and the standard error.
		first coord. second coord.
	mean		1.21	0.79
	standard error	0.0.34	0.034
	As in the above example where the marginal distributions where the same, it
	seems that the result is similar to the non correlated case but with a smaller
	standard error.
				3, 1), W = 0.3 + T ,
	T a Student distribution with 5 degrees of freedom and I a Bernoulli distri-
	bution with parameter 1 5 .
	The Z's, W , I are independent.
	We have simulated 30 times 1000 independent periods of length n = 1 with
	u = 2.

Acknowledgments. We are very grateful to Anatoli Juditsky for his valuable advices. Moreover, we would like to thank the referees for their careful reading of our manuscript and their comments.


C Berthon

S Cordier

O Delestre

M H Le

An analytical solution of Shallow Water system coupled to Exner equation

In this paper, an exact smooth solution for the equations modeling the bedload transport of sediment in Shallow Water is presented. This solution is valid for a large family of sedimentation laws which are widely used in erosion modeling such as the Grass model or those of Meyer-Peter & Müller. One of the main interest of this solution is the derivation of numerical benchmarks to valid the approximation methods.

Introduction

Soil erosion is a consequence of the movements of sediments due to mechanical actions of flows. In the context of bedload transport, a mass conservation law, also called Exner equation [START_REF] Exner | Über die wechselwirkung zwischen wasser und geschiebe in flüssen, Sitzungsber[END_REF], is used to update Email addresses: christophe.berthon@math.univ-nantes.fr (C. Berthon), stephane.cordier@math.cnrs.fr (S. Cordier), Delestre@unice.fr (O. Delestre), mh.le@brgm.fr (M.H. Le).

the bed elevation. This equation is coupled with the shallow water equations describing the overland flows (see [START_REF] Castro Díaz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF] and references therein) as follows:

∂ t h + ∂ x (hu) = 0, (1)
∂ t (hu) + ∂ x (hu 2 + gh 2 /2) + gh∂ x z b = 0, ( 2
) ∂ t z b + ∂ x q b = 0, ( 3
)
where h is the water depth, u the flow velocity, z b the thickness of sediment layer which can be modified by the fluid and g the acceleration due to gravity. The variable hu is also called water discharge and noted by q. Finally, q b is the volumetric bedload sediment transport rate. Its expressions are usually proposed for granular non-cohesive sediments and quantified empirically [START_REF] Grass | Sediment transport by waves and currents[END_REF][START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF][START_REF] Cordier | Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help[END_REF]. Many numerical schemes have been developed to solve system (1-3) (see [START_REF] Cordier | Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help[END_REF] and references therein). The validation of such schemes by an analytical solution is a simple way to ensure their working. Nevertheless, analytical solutions are not proposed in the literature. Up to our knowledge, asymptotic solutions, derived by Hudson in [START_REF] Hudson | Numerical technics for morphodynamic modelling[END_REF], are in general adopted to perform some comparisons with approximated solutions. The solutions are derived for Grass model [START_REF] Grass | Sediment transport by waves and currents[END_REF], i.e q b = A g u 3 , when the interaction constant A g is smaller than 10 -2 . In this paper, we propose a non obvious analytical solution in the steady state condition of flow.

Solution of the equations

In order to obtain an analytic solution, we consider q b as a function of the dimensionless bottom shear stress τ * b (see [START_REF] Cordier | Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help[END_REF]). By using the friction law of Darcy-Weisbach, τ * b is given by

τ * b = f u 2 8(s -1)gd s ,
where f is the friction coefficient, s = ρ s /ρ the relative density of sediment in water and d s the diameter of sediment. The formulae of q b is usually expressed under the form

q b = κ(τ * b -τ * cr ) p (s -1)gd 3 s , (4)
where τ * cr is the threshold for erosion, κ an empirical coefficient and p an exponent which is often fixed to 3/2 in many applications. The expression (4) can be written in the simple form

q b = Au 2p e , (5)
where the effective velocity u e and the interaction coefficient A are defined by

           u 2 e = u 2 -u 2 cr , u 2 cr = τ * cr f 8(s -1)gd -1 , A = κ f 8(s -1)gd p (s -1)gd 3 s . (6)
Remark. The Grass model [START_REF] Grass | Sediment transport by waves and currents[END_REF] is one of the simplest case by using p = 3/2, τ * cr = 0 and an empirical coefficient A g instead of A. The Meyer-Peter & Müler model [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF] is one of the most applied by using p = 3/2, κ = 8, τ * cr = 0.047. The following result is valid for all models rewriting in form (5-6). Proposition 2.1 Assume that q b is defined by (4). For a given uniform discharge q such that τ * b > τ * cr , system (1-3) has the following analytical unsteady solution

                 u 2 e = αx + β A 1/p , u = u 2 e + u 2 cr , h = q/u, z 0 b = - u 3 + 2gq 2gu + C, z b = -αt + z 0 b . (7)
where α, β, C are constants and A, u cr are defined by [START_REF] Hudson | Numerical technics for morphodynamic modelling[END_REF].

Proof. We are here concerned by the smooth solution. In view of the assumption hu = q = cst, equations (1-3) reduce to

∂ t h = 0, ∂ x (q 2 /h) + gh∂ x H = 0, ( 8
) ∂ t H + ∂ x q b = 0, ( 9
)
where H = h + z is the free surface elevation. Differentiating equation ( 8) with respect to t and then equation ( 9) with respect to x, we obtain

∂ xt H = 0, ∂ 2 x q b = 0. ( 10
)
Note that we can write q b = q b (h, q) to have

∂ t q b = ∂ h q b ∂ t h + ∂ q q b ∂ t q = 0, so q b is not time-depending.
Thank to (10), the expression of q b is obtained under the form

q b = αx + β, (11)
where α and β are constant. From (3), we obtain

∂ t z b = -∂ x q b = -α to write z b = -αt + z 0 b (x). (12)
Moreover, from (5) we deduce the effective velocity as follows:

u 2 e = αx + β A 1/p .
Plugging (12) into the momentum equation ( 8) and using a direct calculation, we have

∂ x z 0 b = q u 2 - u g ∂ x u ⇒ z 0 b = - u 3 + 2gq 2gu + C
which concludes the proof.

Remark. As h and u are stationary, the initial condition of ( 7) is (h, u, z 0 b ). Moreover, the solution (h, u) applied to the Grass model is also an analytical solution of the Shallow Water Equations with the variable topography z 0 b . Concerning the shallow-water model, other solutions can be found in [START_REF] Delestre | SWASHES: a library of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies[END_REF].

Numerical experiments

In this section, we consider the analytical solution [START_REF] Delestre | SWASHES: a library of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies[END_REF] applied to the Grass model with q = 1, A g = α = β = 0.005 and C = 1. A relaxation solver is applied to approximate the solution of the model. We will not give here the details of the relaxation solver (for details see [START_REF] Audusse | Sediment transport modelling : Three layer models and relaxation schemes[END_REF]), but just the relaxation model for the equations (1-3). Thus, we solve the following relaxation system:

                       ∂ t h + ∂ x (hu) = 0, ∂ t (hu) + ∂ x (hu 2 + π) + gh∂ x z b = 0, ∂ t π + u∂ x π + a 2 h ∂ x u = 0, ∂ t z b + ∂ x q r = 0, ∂ t q r + b 2 h 2 -u 2 ∂ x z b + 2u∂ x q r = 0,
that is completed with π = gh 2 /2 and q r = q b at the equilibrium. Figure 1 presents the numerical result with J = 500 space cells, a CFL fix condition of 1 and T = 7s. We only notice little difference on the velocity, near the inflow boundary.

Figure 1 .

 1 Figure 1. Comparison between the exact solution and the relaxation method for : the water height and the topography (left) and the velocity (right).

Preprint submitted to the Académie des sciencesDecember 7, 2011