Manon Costa
email: manon.costa@math.univ-toulouse.fr

Sébastien Gadat
email: sebastien.gadat@tse-fr.eu

Lorick Huang
email: lorick.huang@insa-toulouse.fr.

CV@R penalized portfolio optimization with biased stochastic mirror descent

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Given a portfolio of nancial assets, what is the best way to allocate resources under a risk management constraint? Without any xed constraint, it is likely that the asset with the highest expected return would be favored, leading to a trivial optimization problem. However, ethical issues arise when this latter return possesses high variability, since even though in the long run some prots are expected, large losses can occur in between wins. This problematic is especially relevant for institutions with some large capital investments such as banks and insurances.

Hence, the necessity for regulation appears. As a matter of fact, this is precisely the point of the Basel III banking regulation that has been implemented in response to the various past nancial crises. Risks can be quantied with various metrics, and policies developed using one or another metric could lead to dierent investment optimization problems and strategies. In this work, we consider a penalty on the conditional value at risk denoted by CV@R (also referred to as the expected shortfall in the literature) of the portfolio, which stands for a coherent measure of risk of an investment [START_REF] Artzner | Coherent measures of risk[END_REF][START_REF] Rockafellar | Conditional value-at-risk for general loss distributions[END_REF]. We refer to Equation [START_REF] Alfonsi | On the discretization schemes for the cir (and bessel squared) processes[END_REF] in the next section for a formal denition. Heuristically, the CV@R measures how big the average loss is, at a given rate for worst scenarii (say 10% for example). Our problem is therefore to optimize the expected portfolio returns while keeping the risk (measured via CV@R) under a xed threshold.

Portfolio optimization under CV@R constraint has already been investigated in previous works, we can notably cite Rockafellar and Uryasev [START_REF] Rockafellar | Optimization of conditional value-at-risk[END_REF] who develop a Monte-Carlo strategy to approximate the CV@R of the portfolio. These results have been rened in subsequent works of the authors, notably in Krokhmal et al. [START_REF] Krokhmal | Portfolio optimization with conditional valueat-risk objective and constraints[END_REF], Rockafellar et Uryasev [START_REF] Rockafellar | Conditional value-at-risk for general loss distributions[END_REF] and more recently in [START_REF] Rockafellar | Random variables, monotone relations, and convex analysis[END_REF] , but basically the authors of these former works use the approximation of the V@R and CV@R with a batch empirical sum to feed an optimization algorithm via classical tools of convex analysis.

The approximation with a Monte-Carlo empirical mean to handle the V@R and the CV@R can be traced back to Ben Tal and Teboulle [START_REF] Ben Tal | Expected utility, penalty functions and duality in stochastic nonlinear programming[END_REF] and Pluf [START_REF] Ch | Some remarks on the value-at-risk and the conditional value-at-risk[END_REF], who link these two quantities by a min and arg min of a convex function. This article spun a series of articles dealing with the Robbins-Monro approximation of V@R and CV@R. We can cite a series of papers by Bardou, Frikha and Pagès [START_REF] Bardou | Cvar hedging using quantization-based stochastic approximation algorithm[END_REF], [START_REF] Bardou | Computing var and cvar using stochastic approximation and adaptive unconstrained importance sampling[END_REF], and the book [START_REF] Bardou | Recursive computation of value-at-risk and conditional value-at-risk using mc and qmc[END_REF], as well as [START_REF] Frikha | A multi-step richardsonromberg extrapolation method for stochastic approximation[END_REF] all dealing with Robbins Monro stochastic approximation. However, as in [START_REF] Frikha | A multi-step richardsonromberg extrapolation method for stochastic approximation[END_REF], a question still remains when the underlying asset's distribution is unobtainable. A natural solution to that problem would be to approximate the underlying distribution, introducing a bias. This trend of biased stochastic approximation is very natural but presents some challenges in itself, which can be addressed in multiple ways. We can cite the works of [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF] or [START_REF] Laruelle | Stochastic approximation with averaging innovation applied to nance[END_REF] (see also the references therein) for an overview of the techniques used for biased stochastic approximation algorithms.

In mathematical nance, assets are commonly modeled using stochastic dierential equations, and in general, these equations do not have any explicit solutions. A natural approach to this issue is simply to replace the SDE with a discretization scheme. However, a such discretization naturally introduces a bias in the stochastic approximation procedure, and the eect of this bias should be addressed to nely tune the optimization algorithm. Discretization for SDEs has been studied since the pioneering work of Talay and Tubaro [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic dierential equations[END_REF], where strong regularity is assumed on the coecients. As results developed, regularity assumptions have been dropped, and sharp bounds have been derived for both weak and strong error (see e.g. [START_REF] Kloeden | Numerical solution of SDE[END_REF]) on specic stochastic models for specic discretization schemes (see e.g. [START_REF] Alfonsi | On the discretization schemes for the cir (and bessel squared) processes[END_REF][START_REF] Alfonsi | High order discretization schemes for the cir process: application to ane term structure and heston models[END_REF] among others).

In our paper, we assume in addition that one of the asset in our portfolio is risk-less, which allows to model debt obligations, or treasury bonds. In mathematical terms, we consider a stochastic rate, modeled by a Cox-Ingersoll-Ross (shortened as CIR below) [START_REF] Cox | A theory of the term structure of interest rates[END_REF] process and we use this rate for the growth ratio of one component of our portfolio (see Subsection 3.1 below the denition). It can be shown that CIR processes can be calibrated to always stay positive (which is desirable for modeling rates -but sometimes not necessary). Unfortunately, positivity is not preserved with the usual explicit Euler scheme. This fact drives us, when discretizing the CIR, to use a drift-implicit Euler scheme, introduced in [START_REF] Alfonsi | On the discretization schemes for the cir (and bessel squared) processes[END_REF], that is known to preserve positivity (see Section 3.2 and Proposition C.2 below for more details). Thankfully, rate of convergence for the approximation of the CIR process has been widely studied, see e.g. [START_REF] Dereich | An Euler-type method for the strong approximation of the cox ingersoll ross process[END_REF][START_REF] Alfonsi | Strong order one convergence of a drift implicit euler scheme: Application to the cir process[END_REF] and these quantitative results may be exploited to control the bias and its propagation in the specic terms induced by some functional related to the CV@R variational formulation.

Finally, we also impose a reinvestment condition, saying that the strategy must reallocate all funds available. This means that if we have an initial capital of V and n elements in our portfolio, we are looking for a vector (v 1 , . . . , v n ) such that v 1 + • • • + v n = V . Equivalently, we consider u = (u 1 , . . . , u n ) in the simplex: u i ∈ [0, 1] and u 1 + • • • + u n = 1 and look for an optimal set of weights u. Note that this is asking each weights to be positive. This prevents short selling, which is the practice of selling a stock one does not own in hopes of buying it back in the future at a lower price, and is yet another behavior that is regulated by nancial institutions.

In order to determine a strategy u = (u 1 , . . . , u n ) in the simplex, one natural idea would be to derive an algorithm without the constraint u 1 + • • • + u n = 1, and project back the result over the simplex.

However, projection over the simplex tends to favor extremal points, which slows down the convergence of algorithms while projection is itself an additional step that may be numerically costly when some high-dimensional portfolio are considered. To overcome this issue, we set up a stochastic mirror descent, as pioneered by Nemirovskii & Yudin [START_REF] Nemirovskij | Problem complexity and method eciency in optimization[END_REF]. Mirror descent (and its variants) stands for a widely used class of optimization algorithms, whose main interest is to avoid a costly projection into the set of constraints by replacing the Euclidean distance by a Bregman divergence in the second order Taylor expansions. It may be adapted to the stochastic framework with the help of carefully tuned step-size sequences, we can cite among others [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] that initiates the use of a noisy stochastic gradient in the mirror descent, [START_REF] Lan | Validation analysis of mirror descent stochastic approximation method[END_REF] for specic derivations on convex optimization problems on the simplex and more recently Zhou et al [START_REF] Zhou | Stochastic mirror descent in variationally coherent optimization problems[END_REF] that derives some results on the a.s. convergence in the general case with asymptotic pseudo-trajectories, that are common tools of random dynamical systems (see e.g. [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent ows, with applications[END_REF]).

In summary, to nd an investment strategy while managing risk, we will use a stochastic mirror descent with an approximation of the portfolio, combined with a ne approximation for the CV@R.

The discretization of the Portfolio is handled by an implicit Euler scheme. This discretization feeds the approximation of the V@R and the CV@R, which is done through a biased Robbins-Monro algorithm, which in turn feeds a stochastic mirror descent to determine the optimal portfolio allocation. Even if every individual steps have a known (sometimes involved) solution, we show that plugging-in these solutions together yields a strategy that solves the constrained optimization problem. Our main results are as follows:

First, using some sophisticated tools of stochastic approximation (ODE method) and of mirror descent (Fenchel conjugate among others), we establish the convergence of our biased procedure towards the minimum of our penalized criterion.

Second, we establish some new sharp probabilistic upper bounds dedicated to the simulation of the CIR process that are necessary to assess the inuence of the number of assets m on our stochastic algorithm. These bounds are at the cornerstone of our work and of crucial importance for the ne tuning of the numerical implicit Euler scheme.

In what follows, we carefully study the inuence of the number of considered assets and of the discretization step-size on the numerical cost of the overall procedure. We then use our method to construct a Markowitz portfolio envelope [START_REF] Markowitz | Portfolio selection[END_REF] that crosses the risk (built with the CV@R) and the expected return to identify the optimal investment strategy. This article is organized as follows. In Section 2 we detail the optimization problem and the stochastic mirror descent that we propose. We state the almost sure convergence result as well as the non asymptotic bounds that are proved in Appendix B. Then in Section 3 we describe the structure of a portfolio of interest and the simulation strategy. In particular we highlight the link between the step size of the discretization and the bias that will be induced in the SMD algorithm. Finally, we illustrate our results in Section 4 both on simulated and on real data.

2 CV@R α constrained mean value optimization

Financial model and constrained optimization

We consider the relative return of a portfolio, described by a random vector Z = (Z 1 , . . . , Z m ), that is composed of m assets whose capital at the initialization time t 0 = 0 is scaled to 1. We address the question of the optimal investment strategy in order to maximize the mean relative return at a xed time horizon t 0 + T where T = 1, i.e. each random variable Z i corresponds to Z i = Ai(T ) Ai(0) -1 where A i (t) is the price of asset i at time t. An investment strategy corresponds to an allocation of the initial capital, modeled by a collection of m positive weights u = (u 1 , . . . , u m ) that belongs to the m -1 dimensional simplex denoted by ∆ m and dened as:

∆ m := u ∈ R m + : m i=1 u i = 1 .
We are interested in a constrained optimization of the mean of the random variable ⟨Z, u⟩ dened by:

⟨Z, u⟩ = m i=1 u i Z i .
The constraint we are using involves a risk measure, named the conditional value at risk, which quanties the mean value of a loss given that a loss has occurred. More precisely, if α refers to a chosen risk level, V@R α (u) corresponds to the classical statistical α-quantile dened as:

V@R α (u) := sup{q ∈ R : P[⟨Z, u⟩ ≤ q] ≤ α}, while CV@R α , that denes the active constraint we will use, is the mean value of the loss when ⟨Z, u⟩ is below V@R α , namely:

CV@R α (u) = E[-⟨Z, u⟩ | ⟨Z, u⟩ ≤ V@R α (u)]. (1)
Note that here, we implicitly assume that the risk level α is chosen such that the associated quantile V@R α (u) is negative, which implies that CV@R α is a positive quantity. Then, for any xed positive level M , the optimization problem of u -→ E[⟨Z, u⟩] with CV@R α constraints we are interested in, is the problem dened by:

P M := arg max u∈∆m m i=1 u i E[Z i ] : CV@R α (u) ≤ M = arg min u∈∆m - m i=1 u i E[Z i ] : CV@R α (u) ≤ M
For our purpose, we also introduce the unconstrained penalized optimization problem: we consider for any λ > 0 the solution v λ of the following (convex) optimization problem:

Q λ := arg min u∈∆m - m i=1 u i E[Z i ] + λCV@R α (u) .
The (standard) following result is key for our forthcoming analysis.

Proposition 2.1. The collection of the two convex problems (P M ) M >0 and (Q λ ) λ>0 are equivalent. More precisely, for any M > 0 that denes a feasible constraint, a solution u ⋆ M exists such that:

∃λ ⋆ M > 0 u ⋆ M = arg min u∈∆m - m i=1 u i E[Z i ] + λ ⋆ M CV@R α (u) . Moreover, λ ⋆ M is a decreasing function of M . Oppositely, any solution v λ of Q λ solves P M with M = CV@R α (v λ ).
We emphasize that this latter result is the keystone of our forthcoming analysis: instead of solving directly P M , we will instead dene a collection of penalized problems Q λ for several values of λ. Proposition 2.1 then establishes that this approach leads to optimal portfolio allocation under CV@R constraints.

We also highlight that the relationship between M and λ ⋆ M remains unknown and challenging to obtain.

The proof of this Lagrangian formulation is deferred to Appendix A. Using the result of [START_REF] Rockafellar | Optimization of conditional value-at-risk[END_REF][START_REF] Krokhmal | Portfolio optimization with conditional valueat-risk objective and constraints[END_REF] and in particular the convex representation of the CV@R α , we observe that:

CV@R α (u) = min θ∈R ψ α (u, θ),
where ψ is the convex coercive Lipschitz continuous and dierentiable function dened by:

ψ α (u, θ) = θ + 1 1 -α E (⟨Z, u⟩ -θ) + ,
where x + = max(0, x). We emphasize that despite the denition with () + , when Z has an absolutely continuous distribution with respect to the Lebesgue measure, ψ α is dierentiable (we refer to [START_REF] Bardou | Cvar hedging using quantization-based stochastic approximation algorithm[END_REF] for further details on this remark).

Then, we deduce from Proposition 2.1 and the above representation that the collection of the optimization problems (P M ) M >0 are then described equivalently by the convex unconstrained problem:

Q λ = arg min (u,θ)∈∆m×R {p λ (u, θ)}, (2)
where the key function p λ is dened by:

p λ (u, θ) = - m i=1 u i E[Z i ] + λψ α (θ, u).
We emphasize that (p λ ) λ>0 forms a collection of convex functions dened on ∆ m × R, coercive w.r.t. θ and that Q λ enables to avoid the CV@R α constraint with a suitable reparametrization of the problem.

Biased stochastic Mirror Descent

To solve the minimization problem [START_REF] Alfonsi | High order discretization schemes for the cir process: application to ane term structure and heston models[END_REF] we are led to use stochastic approximation theory, that originates in the pioneering works of [START_REF] Robbins | A stochastic approximation method[END_REF] and [START_REF] Robbins | A convergence theorem for non negative almost supermartingales and some applications[END_REF] since the (convex) function p λ is written as an expectation.

However, we encounter here specic diculties. First we need to handle the minimization over the simplex ∆ m . Second, the random variables Z involved in p λ cannot generally be simulated exactly, and it will therefore be necessary to control the bias coming from the stochastic simulation. Lastly, note that even though p λ is dierentiable, it is the expectation of a non-dierentiable function of (u, θ) which will require some specic attention in the following algorithms.

Deterministic case

Mirror Descent (MD below) originates from the pioneering work of [START_REF] Nemirovskij | Problem complexity and method eciency in optimization[END_REF] and permits to naturally handle constrained optimization problems especially when the mirror/proximal mapping is explicit, which is indeed the case for a convex problem constrained on ∆ m (see e.g. [START_REF] Lan | Validation analysis of mirror descent stochastic approximation method[END_REF], [START_REF] Bubeck | Convex optimization: Algorithms and complexity[END_REF]). The MD approach has the nice feature to dene a smooth evolution that lives inside the constrained set without adding some supplementary projection step and pushes the frontiers of the simplex at an innite distance from any point strictly inside ∆ m . Though we will need to encompass the stochastic framework, since p λ is dened through an expectation which is not explicit, for the sake of clarity we describe rst the deterministic version.

For this purpose, we rst introduce the strongly convex negative entropy function over the simplex ∆ m of probability distributions, dened by:

∀u ∈ ∆ m , φ(u) = m i=1 u i log u i .
If ⟨•, •⟩ refers to the standard Euclidean inner product, the Bregman divergence D φ is dened by:

∀(u, v) ∈ ∆ 2 m , D φ (u, v) = φ(u) -φ(v) -⟨∇φ(v), u -v⟩.
This Bregman divergence will induce the natural metric associated to the component u involved in the problem Q λ . In the same way, we also introduce the standard Euclidean square distance over R 2 , which may be viewed as a Bregman divergence. It then leads to a Bregman divergence over pairs (θ, u) and (θ ′ , v) dened as:

∀(θ, θ ′ ) ∈ R 2 , ∀(u, v) ∈ ∆ 2 m , D Φ ((u, θ), (v, θ ′ )) = (θ -θ ′ ) 2 2 + D φ (u, v),
which is the Bregman divergence associated to the strongly convex function:

Φ(u, θ) = φ(u) + θ 2 2 .
We emphasize that the strong convexity of Φ induces the following lower bound:

∀(θ, θ ′ ) ∈ R 2 , ∀(u, v) ∈ ∆ 2 m , D Φ ((u, θ), (v, θ ′ )) ≥ (θ -θ ′ ) 2 2 + 1 2 ∥u -v∥ 2 .
(

) 3
The (deterministic) MD with a step-size sequence (η k ) k≥1 consists in minimizing from k to k + 1 the rst order Taylor approximation of p λ penalized by the Bregman divergence. To alleviate the notations , we will denote by X k = (U k , θ k ) the position of the algorithm and generally x = (u, θ) an element of ∆ m × R. We observe that the MD iterative step corresponds to:

X k+1 = arg min x∈∆m×R ⟨∇p λ (X k ), x -X k ⟩ + 1 η k+1 D Φ (x, X k ) .
A remarkable feature is that this latter minimization can be made explicit :

X k+1 = U k+1 θ k+1 with    U k+1 = U k e -η k+1 ∂u p λ (U k ,θ k ) ∥U k e -η k+1 ∂u p λ (U k ,θ k ) ∥1 θ k+1 = θ k -η k+1 ∂ θ p λ (U k , θ k )
where the rst equation has to be understood within a m dimensional vector structure. For a xed λ > 0, we dene x ⋆ λ = (u ⋆ λ , θ ⋆ λ ) the minimizer of p λ . Following the work of [START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Nemirovskij | Problem complexity and method eciency in optimization[END_REF], it can be shown that an averaged version of the sequence (X k ) k≥1 dened by

χn = n k=0 η k -1 n k=0 η k X k ,
satises the next error bound: Theorem 1 (Proposition 1 in [START_REF] Lan | Validation analysis of mirror descent stochastic approximation method[END_REF]). For any λ > 0, we note by

L = arg max x∈∆m×R {∥∇p λ (x)∥}, then p λ ( χn ) -p λ (x ⋆ λ ) ≤ {∆ 0 Φ } 2 + L 2 n k=1 η 2 k 2 n k=1 η k , where {∆ 0 Φ } 2 = (θ0-V@Rα(u ⋆ λ )) 2 2 + log m.
It is possible to produce several variations around this previous result. It is also possible to assess an upper bound on L 2 , using in particular Equations ( 4) and (5) stated below, which implies

L 2 ≤ m(1 -α + λ) 2 E[∥Z∥ 2 ] + λ 2 α 2 (1 -α) 2 .
Nevertheless, we emphasize that we will not use this upper bound as our study will be more involved.

We should only keep in mind that it is possible to nely tune the step-size sequence to obtain a O(n -1/2 ) upper bound. We refer to [START_REF] Lan | Validation analysis of mirror descent stochastic approximation method[END_REF] for further details.

Stochastic Mirror Descent (SMD) using biased simulations

In our problem, p λ and ∇p λ involves the computation of several expectations, which are not explicit and for which the computational time of approximation may not be reasonable. If we denote by ∂ u and ∂ θ the partial derivatives with respect to u and θ, and using that p λ is expectation of a convex function, we then verify that:

∂ u p λ (u, θ) = -E[Z] + λE[Z1 ⟨Z,u⟩≥θ ] 1 -α , (4)
and

∂ θ p λ (u, θ) = λ 1 - 1 1 -α E[1 ⟨Z,u⟩≥θ ] . (5)
Remark 1. We formulated the above convergence and rate estimation results for ∇p λ . However, these results would stay true as soon as the key function p λ is convex when we can access to a sub-gradient, which has been shown to be compatible with mirror descent.

We assume that we observe a sequence of mutually independent random variables ( Ẑk ) k≥0 that are also sampled independently from the previous positions of the algorithm. The expressions (4) and ( 5) lead to a natural (possibly biased) stochastic approximation of ∇p λ with the help of the sequence ( Ẑk ) k≥0 .

Assuming that the algorithm is at step k at position (U k , θ k ), we introduce the stochastic approximation of the sub-gradients:

     ĝk+1,1 = -Ẑk+1 + λ 1-α Ẑk+1 1 ⟨ Ẑk+1 ,U k ⟩≥θ k ĝk+1,2 = λ 1 -1 1-α 1 ⟨ Ẑk+1 ,U k ⟩≥θ k . ( 6
)
We now describe the necessary assumptions on the sequence ( Ẑk ) k≥0 to build a consistent SMD algorithm.

Assumption 1. Biased simulations Let L( Ẑk+1 ) and L(Z) denote the distributions of Ẑk+1 and Z respectively. We assume that the sequence ( Ẑk ) k≥0 satises both:

W 1 (L( Ẑk+1 ), L(Z)) ≤ δ k+1 , (7)
and

∀u ∈ ∆ m , ∀θ ∈ R, E ⟨Z, u⟩1 ⟨Z,u⟩≥θ -⟨ Ẑk+1 , u⟩1 ⟨ Ẑk+1 ,u⟩≥θ | F k ≤ υ k+1 . (8)
where W 1 stands for the Wasserstein-1 distance.

The sequences (δ k+1 ) k≥0 and (υ k+1 ) k≥0 translate the fact that we may observe some biased realizations of the assets Z at step k, the perfect simulation framework being translated by δ k+1 = υ k+1 = 0 for all integer k. The Wasserstein distance used in δ k+1 is indeed an easiest way to upper bound the bias naturally involved in our approximation procedure. Indeed, Equation [START_REF] Frikha | A multi-step richardsonromberg extrapolation method for stochastic approximation[END_REF] (see below) will involve the Kolmogorov distance instead, denoted by d Kol below, which quanties the dierence between cumulative distribution function. As the discretization of S.D.E. with (implicit) Euler scheme are well understood in terms of Wasserstein-1 distance, and because of the straightfoward inequality d Kol ≤ 2 √ W 1 , we prefered to introduce Assumption 1 in terms of W 1 .

We emphasize that both (δ k+1 ) k≥0 and (υ k+1 ) k≥0 might heavily depend on m the dimension of the vector Z. A specic example will be detailed below in our work.

The method we propose is then dened with Algorithm 1. The next result states the asymptotic almost sure convergence of the sequence (X k ) k≥0 constructed in Algorithm 1 towards the target point

x ⋆ λ = (u ⋆ λ , V@R α (u ⋆ λ )).
     ĝk+1,1 = -Ẑk+1 + λ 1-α Ẑk+1 1 ⟨ Ẑk+1 ,U k ⟩≥θ k ĝk+1,2 = λ 1 -1 1-α 1 ⟨ Ẑk+1 ,U k ⟩≥θ k .
Update the algorithm X k+1 = arg min x∈x∈∆m×R ⟨ĝ k+1 , x -X k ⟩ + 1 η k+1 D Φ (x, X k ) using:

X k+1 = (U k+1 , θ k+1 ),    U k+1 = U k e -η k+1 ĝk+1,1 ∥U k e -η k+1 ĝk+1,1 ∥ 1 θ k+1 = θ k -η k+1 ĝk+1,2 . end Algorithm 1: Biased SMD i) k≥0 η k+1 (p λ (X k ) -min(p λ )) ≤ k≥0 η k+1 ⟨∇p λ (X k ), X k -x ⋆ λ ⟩ < +∞ a.s.
ii) The Cesaro average Xη k dened by

Xη k := k i=0 η i -1 k i=0 η i X i (9)
is almost surely convergent and

p λ ( Xη k ) -→ min(p λ ) a.s.
iii) Assume that:

k≥0 δ k+1 + υ k+1 < +∞,
then the sequence (X k ) k≥1 almost surely converges and veries that:

lim k -→+∞ p λ (X k ) = min(p λ ) a.s.
The results obtained in i) and ii) are rather classical in the stochastic optimization community. The diculty here is to derive a Lyapunov function and then to quantify the eect of the bias induced by

Ẑk+1 satisfying [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF] and [START_REF] Bally | The law of the Euler scheme for stochastic dierential equations[END_REF]. To do so, we will follow the arguments developed in [START_REF] Lan | Validation analysis of mirror descent stochastic approximation method[END_REF] in order to obtain a recursion expression but we will need to use some ad-hoc supplementary steps to handle the bias and to use the Robbins Siegmund theorem. The last item is much more dicult to obtain, roughly speaking we use the O.D.E. method introduced in [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent ows, with applications[END_REF]. Even though it is a standard tool in stochastic optimization, applying this method to our framework is challenging because of the mirror descent O.D.E. and it deserves some sophisticated additional developments that are far from being a straightforward application of [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent ows, with applications[END_REF].

In particular, we will need to use and adapt some ne results of [START_REF] Mertikopoulos | On the convergence of gradient-like ows with noisy gradient input[END_REF] to our framework. The proof is given in Appendix B.

Asymptotic results such as the one stated in Theorem 2 are commonly criticized because algorithms are always ran within a nite number of iterations. In particular, Theorem 2 does not provide any insight on the rate of convergence of the method, which may crucially inuence the practical usefulness of the algorithm. Nevertheless, it shows that Algorithm 1 converges and minimizes the function p λ up to a Cesaro averaging procedure (see i) and ii) of Theorem 2). It may be shown indeed that the algorithm converges without Cesaro averaging at the price of a more stringent condition on the terms involved in the bias of the algorithm quantied by δ k+1 and υ k+1 .

To cope with the legitimate asymptotic criticisms, we also state non asymptotic theoretical guarantees for the value of the objective function.

Theorem 3 (Finite-time guarantees). Let us consider (X k ) k≥0 dened in Algorithm 1 and its Cesaro averaging counterpart introduced in ( 9), then for any n > 1,

E[p λ ( Xη n )] -p λ (x ⋆ λ ) ≤ n-1 j=0 η j+1 -1 D 0 Φ + n-1 k=0 a k+1 D k Φ + b k+1 ,
where

D k Φ = E[D Φ (x ⋆ λ , X k )]
and

a k+1 = 2η k+1 2 δ k+1 + δ k+1 + λυ k+1 1 -α , b k+1 = C η 2 k+1 + a k+1 .
The constant C > 0 appaering in b k+1 is precised in [START_REF] Markowitz | Portfolio selection[END_REF].

Finally, we aim at applying these guarantees in the case of a nite horizon of computation n and in this view, D k Φ can be controlled using a recursion inequality (35). We consider the specic case of the SMD with a constant step-size sequence stopped at iteration n:

η k+1 = η > 0, ∀ 0 ≤ k ≤ n.
We will also assume a constant upper bound of the bias in the simulation of the random variables Ẑk : we denote by ω its resulting impact in the SMD. This is legitimate since we can reduce this bias with the use of an arbitrarily small discretization step-size, which of course harms the computational cost. More precisely we consider xed values of δ k+1 and υ k+1 such that:

2 δ k+1 + δ k+1 + λυ k+1 1 -α = ω > 0, ∀1 ≤ k ≤ n. Corollary 2.1. For a given n ∈ N, if (η, ω) are chosen such that η = ∆ 0 Φ 2 √ n + 1 and ω = 1 √ n + 1∆ 0 Φ with {∆ 0 Φ } 2 = (θ0-V@Rα(u ⋆ λ )) 2 2
+ log m, then there exists C > 0 large enough such that:

E[p λ ( Xη n )] -p λ (x ⋆ λ ) ≤ C |θ 0 -V@R α (u ⋆ λ )| + √ log m √ n + 1 .
This previous corollary is built using the optimal tuning of the parameters η and ω derived from our proof that is detailed in Appendix B.2. Therefore, these values may be seen as purely theoretical as we do not exactly know the value of ∆ 0 Φ . Nevertheless, any upper bound of ∆ 0 Φ may be used to derive a strategy and an upper bound of the excess risk (see Equation (37) in Appendix (B.2)). In the meantime, the choice:

η = 1 2 √ n + 1 and ω = 1 √ n + 1 , yields: E[p λ ( Xη n )] -p λ (x ⋆ λ ) ≤ C |θ 0 -V@R α (u ⋆ λ )| 2 + log m √ n + 1 ,
which is slightly worse than the upper bound stated in Corollary (2.1).

3 Portfolio hedging under CV@R α constraint

We aim to apply our optimization strategy to the specic situation of portfolio hedging. As a consequence, we precise the structure of Z, a portfolio of m nancial assets, and detail the strategy of biased simulation that can be used. We will illustrate this setting using numerical simulations in Section 4. We are exactly in the eld of application described in Section 2.

Description of the portfolio dynamics

In what follows, we consider the situation where Z contains m = m ′ + 1 assets: a return Y obtained as the baseline short-term interest rate Cox-Ingersoll-Ross process or CIR, (r t ) t≥0 with no drift, and a family S = (S 1 , . . . S m ′ ) of m ′ = m -1 geometric Brownian motions that encode some risky assets in the portfolio.

Recall that the CIR has been introduced in [START_REF] Cox | A theory of the term structure of interest rates[END_REF] as a diusion process and that it is commonly used for the description of the dynamics of interest rates. The CIR short rate model induces a trajectory t -→ r t , whose stochastic dierential equation depends on a triple (a, b, σ 0 ) and is given by:

dr t = a(b -r t )dt + σ 0 √ r t dB 0 (t), (10)
where (B 0 (t)) t≥0 stands for a standard real Brownian motion. The parameter b stands for the long-time mean of the short rate while a quanties the strength of the mean-reversion eect. The volatility σ 0 is multiplied by √ r t , and the condition 2ab > σ 2 0 guarantees that the interest rate remains positive with probability 1. We refer to [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF] for further details. This almost sure positivity motivated our decision to use the CIR process instead of the Vasicek one. Moreover, the CIR model also incorporates both the mean reversion and the conditional heteroscedasticity since the volatility of the short rate process is increased when the short rate increases.

The assets Z t = (Y t , S 1 t , . . . , S m ′ t ) are then described by the following system of stochastic dierential equations:

∀t ≥ 0 dY t = r t Y t dt, dS i t = µ i S i t dt + σ i S i t dB i (t), ∀i ∈ {1, . . . , m ′ }, (11)
where B = (B 0 , B 1 , . . . , B m ′ ) refers to a multivariate Brownian motion with correlated components. For the sake of simplicity, these components are assumed to satisfy:

E[B i (t)B j (t)] = ρ i,j t,
but more general correlation structures could be handled in our framework with further eorts. The correlation matrix is the symmetric positive matrix, denoted by Σ = (ρ i,j ) 1≤(i,j)≤m .

The rst process (Y t ) t≤1 corresponds to a neutral risk process whereas each geometric Brownian motion (S i t ) t≤1 corresponds to a specic risky asset parametrized by known parameters µ i and σ i . We refer to Pitman and Yor [START_REF] Pitman | A decomposition of bessel bridges[END_REF] and to Gulisashvili and Stein [START_REF] Gulisashvili | Asymptotic behavior of the stock price distribution density and implied volatility in stochastic volatility models[END_REF] for several details on these classical processes used (among others) for portfolio modeling. Below we also assume that all the parameters that describe the CIR evolution and the portfolio dynamics (see Equations ( 10) and ( 11)) are known and we are simply interested in the optimal hedging strategy, i.e. we are looking for the optimal u ⋆ λ associated to the opti- mal mean return penalized by λCV@R α dened in Equation ( 2) at time T = 1. As indicated in Section 2, we will use a Stochastic Mirror Descent strategy and we then need to sample some realizations of Z at time T = 1. Before detailing our sampling strategy, we summarize our assumptions on the parameters.

Assumption 2. Assumptions on the portfolio parameters

We assume that: i) the CIR parameters satisfy ab > σ 2 0 and a > 2 √ 2σ 0 . ii) the correlation matrix of the Brownian motions Σ in invertible.

These assumptions on the coecients dening the CIR ensure a control of the L 2 moment of the weak error rate as well as exponential integrability of the integral of the CIR for all time t (see Appendix C).

Simulation of the portfolio

We are interested in an ecient simulation method that satises Assumption 1: both a good approximation of the law of Z 1 (7) and a bias upper bound of the form [START_REF] Bally | The law of the Euler scheme for stochastic dierential equations[END_REF].

For this purpose, we emphasize that each G.B.M. used in Equation ( 11) may be simulated exactly since an exact representation of S i t is available:

∀t ≥ 0 S i t = S i 0 exp µ i t + σ i B i (t) .
In the meantime, we can also observe that some exact simulation method exists for the CIR process sampled at a given xed time. Nevertheless, the CIR process induces some numerical diculties.

First, the SDE ( 10) is not explicitly solvable at any time t ∈ [0, 1] so that the integral of the CIR between 0 and 1 needs to be approximated. Second the correlations between the dierent Brownian components of the portfolio are hardly compatible with the existing exact simulation methods of the CIR: the presence of correlated noise, which is commonly accepted in nancial modeling, is hardly tractable with the use of an exact simulation of the CIR (with the help of χ 2 square distributions, see, e.g. [START_REF] Alfonsi | Ane diusions and related processes: simulation, theory and applications[END_REF]).

Third, one viable strategy to obtain an approximation could be the use of a discretization scheme as the implicit Euler or Milstein schemes with a xed step-size h > 0 for example, thanks to a recursion of the form:

(r (k+1)h , Y (k+1)h ) = F (r kh , Y kh , h, ξ k+1 ),
where F corresponds to an iterative update that uses the position (r kh , Y kh ) and ξ k+1 is a two dimensional Gaussian innovation with a specic covariance. Nevertheless, such an attractive point of view induces some theoretical diculties because of the form of the function Z1 ⟨Z,u⟩≥θ we need to approximate. This function is non-smooth and unbounded, which generates some technical diculties when trying to use or adapt classical approaches on weak discretization errors.

We are led to develop our own ad-hoc simulation of the portfolio and prove its associated properties on the approximation.

Correlated Brownian motion For this purpose, we introduce some independent standard Brownian paths W (t) = (W 0 (t), W 1 (t), . . . , W m ′ (t)) and use the Cholesky decomposition to encode the correlations and recover correlated Brownian motions B(t) = (B 0 (t), B 1 (t), . . . , B m ′ (t)).

More precisely, the Cholesky decomposition used on Σ yields a matrix L such that:

LL T = Σ. (12)
We now dene B(t) = LW (t).

There are two important consequences of constructing (B i (t)) t≥0 for i = 0, . . . , m ′ in this way. The rst one is that the Cholesky method yields a lower triangular matrix L whose rst line is 1, 0, . . . , 0:

L =       1 0 . . . 0 ℓ 01 ℓ 11 . . . . . . . . . . . . 0 ℓ 0m ′ . . . . . . ℓ m ′ m ′       =      1 0 . . . 0 ℓ 01 . . . L ℓ 0m ′      .
In particular, we see that we can take B 0 (t) = W 0 (t) and W (t) = (W 1 (t), . . . , W m ′ (t)) is independent of W 0 (t). The m -1 components may be simply written as:

           B 1 (t) = ℓ 01 W 0 (t) + ℓ 11 W 1 (t) B 2 (t) = ℓ 02 W 0 (t) + ℓ 12 W 1 (t) + ℓ 22 W 2 (t) . . . B m ′ (t) = ℓ 0m ′ W 0 (t) + ℓ 1m ′ W 1 (t) + • • • + ℓ m ′ m ′ W m ′ (t).
The second important and standard consequence is that these stochastic processes (B i (t)) t≥0 for i = 0, . . . , m ′ are indeed Brownian motions with the desired covariance matrix.

Geometric Brownian Motion simulations From now, we consider t = 1 and we alleviate the notations by noting S i 1 = S i , B i (1) = B i and W i (1) = W i . The asset i at time 1 for all i ∈ {1, . . . m ′ } is a Geometric Brownian motion driven by (B i (t)) t≥0 , and we have:

S i = S i 0 exp µ i - (σ i ) 2 2 + σ i B i = S i 0 exp µ i - (σ i ) 2 2 + σ i ℓ 0i W 0 + ℓ 1i W 1 + • • • + ℓ ii W i = S i 0 exp µ i - (σ i ) 2 2 + σ i ( L W ) i independent of W 0 ×e σ0ℓ0iW0 , ( 13
) since ℓ 1i W 1 + • • • + ℓ ii W i = ( L W ) i when W = (W 1 , . . . , W m ′ ).
Interest rate simulation Finally we focus on a discretization method for the CIR process. We mention that discretizing the CIR process leads to some theoretical issues, as the coecients in the SDE are not uniformly elliptic and bounded, as assumed in the seminal works of Bally and Talay [START_REF] Bally | The law of the Euler scheme for stochastic dierential equations[END_REF]. Besides, a classical explicit Euler scheme generates positivity issues (because of the square root). However, many authors, notably Alfonsi [START_REF] Alfonsi | On the discretization schemes for the cir (and bessel squared) processes[END_REF][START_REF] Alfonsi | High order discretization schemes for the cir process: application to ane term structure and heston models[END_REF] proposed implicit Euler schemes and provided weak and strong error rates in the previously mentioned works. We choose to use an drift-implicit Euler scheme which was introduced by [START_REF] Alfonsi | On the discretization schemes for the cir (and bessel squared) processes[END_REF] and studied in numerous articles [START_REF] Dereich | An Euler-type method for the strong approximation of the cox ingersoll ross process[END_REF][START_REF] Alfonsi | Strong order one convergence of a drift implicit euler scheme: Application to the cir process[END_REF] . The drift-implicit Euler scheme on a discrete time grid (kh) 0≤k≤N can be written by considering the SDE satised by y t = √ r t which leads to

ŷ(k+1)h = ŷkh + 4ab -σ 2 0 8ŷ (k+1)h - a 2 ŷ(k+1)h h + σ 0 2 ∆B (k) 0 ,
where ∆B (k) 0 = B 0 ((k + 1)h) -B 0 (kh). This implicit scheme can be solved explicitely on R + from iteration k to iteration k + 1, which ensures the positivity of the scheme (r kh ) 0≤k≤N . In particular, the update from kh to (k + 1)h is given by:

r(k+1)h = √ rkh + σ0 2 ∆B (k) 0 2(1 + ah 2 ) + √ rkh + σ0 2 ∆B (k) 0 2 4(1 + ah 2 ) 2 + (4ab -σ 2 0 )h 8(1 + ah 2 ) 2 . ( 14
)
Convergence results for this discretization scheme will be recalled in Appendix C.

We then use the integral representation Y t = Y 0 exp t 0 r s ds and propose to approximate Y 1 with a Riemann integral approximation between 0 and 1:

Îh := 1 N N k=1
rkh .

This leads to the denition of our approximation:

Ŷ (h) 1 := Y 0 exp( Îh ) = Y 0 exp 1 N N k=1 rkh . (15)
Overall simulation algorithm Combining all these steps, we are led to dene the following algorithm to compute a biased simulation of Z. for k = 1, . . . , N do Compute the recursive approximation rkh using Equation (14). Algorithm 2: Approximation of the portfolio Ẑ1

In order to propose a convergent algorithm derived from Theorem 2, we need to derive an upper bound of the bias induced by Algorithm 2 involving the CIR approximation and the Riemann approximation of its integral computed in Equation [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent ows, with applications[END_REF]. The accuracy of our simulations depends on the number of discretization points N = h -1 sampled between 0 and 1 in Algorithm 2 used to build the Riemann approximation [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent ows, with applications[END_REF]. It is necessary to assess the accuracy of our method in terms of the value of h to make [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF] and ( 8) explicit. The simulation given in Algorithm 2 satises the next statements.

Proposition 3.1. Assume that the portfolio parameters satisfy Assumption 2 then using the above discretization scheme presented in Algorithm 2, a constant C exists (dependent on the CIR parameters) such that:

W 1 (L( Ẑ1 ), L(Z 1 )) = W 1 (L( Ŷ (h) 1 ), L(Y 1 )) ≤ C √ h.
Proposition 3.2. Assume that the portfolio parameters satisfy Assumption 2 then for any e > 0, there exists a constant K e independent of h and m such that:

E ⟨Z 1 , w⟩1 ⟨Z1,w⟩≥θ -⟨ Ẑ1 , w⟩1 ⟨ Ẑ1,w⟩≥θ 2 ≤ K e √ me {σ + } 2 m 2 4e 2 h 1 6 -e .
where σ + = sup 1≤i≤m σ i < +∞.

The proofs of these two results is postponed to Appendix C. We emphasize that these two satellite results contain some sharp new estimations of numerical probability nature on the Euler scheme used for the simulation of our portfolio, and may be of independent interest for applications in mathematical nance. They are not so easy to obtain in particular regarding the inuence of the number of assets m and need some technical computations to obtain this last dependency that is crucial in our stochastic optimization procedure to design an ecient bias reduction strategy.

Stochastic Mirror Descent with sampling approximation

We nally aggregate the optimization procedure described in Algorithm 1 with our sampling scheme in Algorithm 2 and we address the problem of adapting the step size of the discretization scheme to the step size of the SMD approximation in light of Theorems 2, 3 and Propositions 3.

1
h 1/4 + h 1 6 -e ∼ n -1/2 ,
which entails that we could choose a discretization step-size close to n -3 .

SMD with decreasing step-size sequence Let us now choose adequate step-size sequences in order to obtain an a.s. convergent optimization algorithm (Theorem 2). Recall that the sequence (η k ) k≥1 scales the step-size of stochastic gradient descent and assume that

η k = k -α , α ∈ ( 1 2 , 1].
We now chose the step sequence for approximating the CIR process as:

h k = h (m) 0 k -β , β > 0.
According to Propositions 3.1 and 3.2, using the notations of Theorem 2, we deduce that:

δ k = k -β 2 and v k = K e √ me {σ + } 2 m 2 4e 2 h -β( 1 6 -e) .
Assumption in Theorem 2 k≥0 η k+1 ( δ k+1 + υ k+1 ) < +∞ now reads:

k≥1 k -α k -β 4 + k -β( 1 6 -e) < ∞,
which is equivalent to:

α + β 6 > 1, with α ∈ ( 1 2 , 1].
With these conditions on α and β, we can now derive the non-asymptotic bound Theorem 3. We set ϵ = α + β 6 -1 > 0, using the notations of Theorem 3, we deduce that:

a k+1 ∼ k -1-ϵ and b k+1 ∼ k -(2α∧(1+ϵ)) .
Using comparisons between series and integral, we nd that:

D k ϕ ≤ D 0 Φ exp(k -ϵ )) + k -(2α-ϵ+1) .
If we use the notation ≲ that denotes an inequality up to a constant C m which heavily depends on the number of assets in the portfolio, we then obtain that:

E[p λ ( Xη n )]-p λ (x ⋆ λ ) ≤ 1 n-1 j=0 η j+1 D 0 Φ + n-1 k=0 a k+1 D k Φ + b k+1 ≲ 1 n 1-α D 0 Φ + n-1 k=0 k -1-ϵ [D 0 Φ exp(k -ϵ )) + k -(2α+ϵ-1) ] + k -(2α∧(1+ϵ)) ≲ n α-1 D 0 Φ 1 + exp(n -ϵ ) + n -α-2 + n -α∧(1+ϵ-α) ≲ n α-1 D 0 Φ + n -α∧ β 6 .
We point out that choosing α = 1/2 and β > 3 allows for a convergence rate O(n -1/2 ) which is similar to the previous case.

Simulations

In this brief numerical paragraph, we illustrate the behavior of our algorithm on both synthetic and real datasets. The simulations have been driven using Python 3.7.13 on a standard computer.

Description of the simulation study

Bregman divergence vs Euclidean projection In order to assess the specicity and eciency of our algorithm, we also consider the projected stochastic gradient optimization that is obtained after a stochastic gradient step-size + a Euclidean projection on the simplex of probability distribution. The projection has been exactly computed with the help of the standard recursive method derived from the Lagrangian formulation. Below, we will refer to SMD for the stochastic mirror descent and PSGD for the projected stochastic gradient descent (whose denition is straightforward (see [START_REF] Duchi | Ecient projections onto the l 1 -ball for learning in high dimensions[END_REF][START_REF] Fujino | Large-scale multiclass support vector machine training via euclidean projection onto the simplex[END_REF][START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF]). In what follows, we have used the recent algorithm introduced in [START_REF] Duchi | Ecient projections onto the l 1 -ball for learning in high dimensions[END_REF] to project any points into the simplex of probability measure for which a fast implementation is available in [START_REF] Fujino | Large-scale multiclass support vector machine training via euclidean projection onto the simplex[END_REF] and for which a complexity bound of O(m 2 ) is established in [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF]. That being said, it is also shown in [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF] that O(m 2 ) is the worst-case complexity of the method while in practice the observed complexity is only O(m). A such complexity is also assessed in expectation when using a random pivot strategy as proposed in [START_REF] Kiwiel | Breakpoint searching algorithms for the continuous quadratic knapsack problem[END_REF]. As reported in Table 1 of [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF], the observed complexity of several projection methods, including the one of [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF], is O(m).

We have chosen to compare SMD with PSGD with the projection method of [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF] as indicated in Section 4 of [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF], their state of the art algorithm generally outperforms the other methods for reasonable vectors to be projected into the simplex.

SMD vs. MC-MD We also compare our algorithm with a batch Monte-Carlo Mirror Descent, e.g.

we use a mini-batch simulation instead of a single one at each step of the algorithm. Therefore, in the Monte-Carlo Mirror Descent (MC-MD), an additional loop is included in order to perform a Monte-Carlo estimation of the gradients (4) and [START_REF] Artzner | Thinking coherently[END_REF]. Of course, a such strategy improves the quality of the estimation but has a supplementary cost in terms of simulation times. We present below a comparison of the two methods on synthetic data where we will compare the results of SMD and MC-MD in terms of numerical costs in seconds.

Diversication and inuence of λ

To undertake the inuence of the penalty parameter λ, and in particular its eect on the composition of the optimal portfolio, we vary its value in a regular grid

λ min + s λmax-λmin Λ 0≤s≤Λ
dened with the help of a lower and upper value (λ min , λ max ) and a number of points Λ+1 in the grid. According to Proposition 2.1, it is expected that large values of λ essentially favor low risk portfolio composition whereas small values of λ allow for risky portfolio with higher expected returns.

Ecient frontier construction Since the value of λ has to be determined from a practical point of view, we decided to use the Markowitz ecient frontier approach [START_REF] Markowitz | Portfolio selection[END_REF]. The ecient frontier was initially introduced with the standard deviation as a natural measure of the risk of the portfolio. This curve was initially obtained by representing the risk in the X-axis and the Expected Return in the Y-axis.

Below, we replace the volatility by the CV@R as a natural risk measure as proposed in [START_REF] Krokhmal | Portfolio optimization with conditional valueat-risk objective and constraints[END_REF] to build our For convenience, we denote

p λ (u, θ) = -J(u) + λψ α (θ, u) where J(u) = - m i=1 u i E[Z i ].
We then draw the frontier obtained (CV @R α (û λs ), J(û λs )) that may be computed on-line all along the iterations of the SMD algorithm with a Cesaro averaging strategy: if (X k ) k≥0 = (u k , θ k ) k≥0 corresponds to a random sequence built with Algorithm 1 stopped at iteration n and if ( Ẑk ) k≥0 is the biased porfolio simulation sequence used all along Algorithm 1, then we may approximate the expected return J(û λs ) and the risk CV @R α (û λs ) by:

1 n n k=1 ⟨u k , Ẑk ⟩ and 1 n n k=1 θ k + 1 1 -α ⟨ Ẑk , u k ⟩ -θ k + .
We emphasize that since this construction requires the optimization of the portfolio for each value of λ, it is mandatory to design an ecient and fast procedure to solve each optimization problem, which legitimates the use of our stochastic algorithm.

Finally, we propose to use the Sharpe ratio criterion (see e.g. [START_REF] Sharpe | Mutual fund performance[END_REF]) to select the optimal portfolio among all the portfolio found on the Markowitz frontier. More precisely, if 1 CIR refers to the risk-free portfolio hedging 1 CIR = (1, . . . , 0), the Sharpe ratio criterion selects the portfolio that maximizes the ratio between the dierence J(û λs ) -J(1 CIR ) and the risk measured by CV @R α (û λs ).

Synthetic dataset

In this rst set of simulations, we implement our SMD with purely articial assets that are driven with a CIR and several GBM, whose parameters are chosen such that the mean reward / Expected Return (ER for short below) and the individual CV@R monotonically vary together. Figure 1 represents one simulation run of the assets we used. At this stage, we do not use any correlation between assets, which is of course a very simplifying assumption that is typically false for nancial application purposes. In a second stage in Section 4.3, we will consider correlated assets.

Comparison between SMD and PSGD Our rst experiment consists in comparing the behaviour of the SMD and the PSGD. The value of λ used here is λ = 0.9.

Figure 2 shows that there is no real evidence that SMD and PSGD lead to very dierent results in terms of the rate of convergence of p λ (X k ) -→ min(p λ ), even though as indicated in Figure 2 the limiting weights reached by the two algorithms may not be the same. As indicated in Figure 3, the performances in terms of ER and CV@R are equivalent, which indicates that the convex function p λ may have an Even though both algorithms achieve almost the same results in terms of weights, CV@R and ER, the numerical cost associated to the Euclidean projection is important and prevents from the use of the PSGD with a large number of assets. Table 1 provides a brief comparison of the the computational time ratio between PSGD and SMD, which demonstrates that SMD with several dozen of assets requires generally at least 5 times less computations when compared to the update with PSGD. This phenomenon is of course amplied when m increases since the projection algorithm requires more and more operations in larger dimensions. Table 1: Ratio of computation with n iterations on m assets between PSGD and SMD (rounded at the closest integer).

All the more, as indicated in Figure 3 (the same conclusion can be drawn from Figure 2), the speed of convergence of SMD seems fastest than the one of PSGD from a numerical point of view. Of course, this remark is purely based on empirical observations from numerical simulations and not on a mathematical results as the theoretical rates of convergence of the methods are equivalent (both rates evolve as n -1/2 where n is the number of iterations of the algorithms and as n -1 with the supplementary Ruppert-Polyak averaging strategy). This is certainly due to the fact that the averaged SMD algorithm evolves smoothly over the simplex while the averaged PSGD approach is much more irregular due to the sequence of gradient step outside + projection inside the simplex. Finally, this last remark associated with the supplementary computational time due to the projection over the simplex (instead of the instantaneous mirror update) amplies the benets of SMD when compared to PSGD regarding the speed of computation. This last remark is inline with some recent observations in the machine learning community (see e.g. [START_REF] Zhou | Stochastic mirror descent in variationally coherent optimization problems[END_REF]).

Comparison between SMD and MC-MD We now present the result of the comparison of SMD and MC-MD. To take into account the computational cost of the MC step, we represent the estimation of the algorithm with respect to the computational time (in second). In these simulation, we draw 10 independent realization of the portfolio in the MC step. We observe that the supplementary accuracy of the estimation of the gradient, obtained at each iteration, is not sucient to overcome the additional supplementary cost (in terms of seconds) generated by the Monte-Carlo step. As indicated by the right hand side of Figure 4, it appears that both SMD and MC-MD converge to the same portfolio when we push the number of iterations (e.g. the computational time), which lead to the same value of the V@R of the portfolio. We have compared in Figure 4 the learning speed of SMD and of MC-MD of the optimal weights and of the optimal V@R. Our simulations seem to recommend the use of a single stochastic simulation per iteration instead of the use of a mini-batch MC strategy for the problem we are studying. Eect of the penalty parameter λ In Figure 5, we assess the inuence of the penalty coecient λ involved in p λ on the optimal solutions obtained with the SMD. The experiments have been set up with 3 and 8 assets in the portfolio. The uctuations of the curves in Figure 5 are due to the early stopping strategy we adopted to limit the computational cost, but the following conclusions are highly trustable since the evolution of each curve is rather clear. We have indicated in the legend of each subgure in Figure 5 the ER of each individual asset that is estimated on-line with our algorithm, as well as the CV@R. We observe that for small values of λ, the optimal portfolio with 3 and 8 assets favor the most risky ones and do not weight the risk-free component. It is of course the opposite behaviour for large values of λ when essentially the CIR component gathers the essential weight of the portfolio. Finally, in the intermediary regime, the CV@R penalty generates a product diversication, which is recovered both on the left and on the right of Figure 5, which is also exactly the objective of the CV@R penalty.

Real dataset

We briey detail how we shall use our algorithm to estimate optimal portfolio. This short numerical paragraph should be understood as a proof of concept study since we believe that it illustrates the ability of our method. We use the daily dataset of the Treasury Yield 5 Years ( ˆFVX) to calibrate the null risk asset as well as several assets that are detailed in Figure 6, between 2014 and 2016 (to avoid the period of negative interest rates). We model the ˆFVX time series as a daily realization of a CIR stochastic model whereas risky assets are considered as correlated GBM. Parameters are estimated following the optimal martingale method as it is reported in [START_REF] Going-Jaeschke | Parameter estimation and bessel processes in nancial models and numerical analysis in hamiltonian dynamics[END_REF] as it outperforms the discretized log-likelihood maximization approach (for discretely observed CIR realizations). We refer to Equation 4.2 for the drift coecients and Equation given at the end of Section 4.3 of [START_REF] Going-Jaeschke | Parameter estimation and bessel processes in nancial models and numerical analysis in hamiltonian dynamics[END_REF].

Concerning the GBM, we estimate their coecients (µ, σ and ρ) with the log-return series, which is a standard procedure described among other in Section 3.2.3 of [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF] for example. Markowitz ecient frontier with CV@R Once the parameters are estimated, we assume that their values are kept xed for the next daily observation so that we are able to use our SMD Monte-Carlo approach to estimate the optimal portfolio associated to these nancial assets. We then use our SMD to estimate the Markowitz ecient frontier, adapted with the CV@R as a measure of risk of the portfolio.

We represent we obtain in Figure 7 for two dierent level of quantiles (α = 5% and α = 1%).

We use for the construction of our frontier a discretization of the SDE and compute the return of the portfolio for T=30 days. The red point represents the position of the "optimal" portfolio according to the Sharpe ratio. The computational cost for the construction of the Markowitz frontier with 10 assets and a time window h ∝ n -3 appears to be reasonnably fast since we compute these frontiers on standard PC with Python 3.0 in less than 10 minutes each. As indicated in Figure 7, SMD produces each time an optimal portfolio that really balances the CV@R and the ER. Lemma B.1 (Three points lemma). For any triple of points (x, y, z), one has:

D Φ (x, z) = D Φ (x, y) + D Φ (y, z) -⟨∇Φ(z) -∇Φ(y), x -y⟩.
Lemma B.2 (Gradient of the Bregman divergence). For any pair of points (x, y), one has:

∇ x D Φ (x, y) = ∇Φ(x) -∇Φ(y).
With the help of these lemmas, we derive the proof of the almost sure convergence result stated in Theorem 2. The main issue generated by our biased simulation setting is that ĝ involved in Equations ( 4), ( 5) and ( 6) acts non-linearly on 1 ⟨ Ẑk+1 ,U k ⟩≥θ k in our stochastic approximation term, which lead to a signicant amount of theoretical diculties.

Proof of Theorem 2. Let us recall that throughout the proof, the position of the algorithm is

X k = (U k , θ k
) and that we denote by x = (u, θ) a generic point of the set X = ∆ m × R.

Proof of i).

The proof is divided into three steps. We rst write a biased descent inequality, and then control the size of the bias induced by our model before using the Robbins Siegmund Theorem.

Step 1: Biased descent inequality. The main point here is to obtain a recursive inequality on

D Φ (x ⋆ λ , X k )
, where x ⋆ λ is a point that minimizes p λ . Our starting point is the denition of X k+1 :

X k+1 = arg min x∈X ⟨ĝ k+1 , x -X k ⟩ + D Φ (x, X k ) η k+1 ,
where we recall that ĝk+1 is the stochastic approximation of the sub-gradients:

     ĝk+1,1 = -Ẑk+1 + λ 1-α Ẑk+1 1 ⟨ Ẑk+1 ,U k ⟩≥θ k , ĝk+1,2 = 1 -1 1-α 1 ⟨ Ẑk+1 ,U k ⟩≥θ k .
The rst order condition entails that:

∀x ∈ X η k+1 ⟨ĝ k+1 , x -x k+1 ⟩ + ⟨∇ x D Φ (X k+1 , X k ), x -X k+1 ⟩ ≥ 0.
Using Lemma B.2 on the second term, we deduce that:

∀x ∈ X ⟨∇Φ(X k+1 ) -∇Φ(X k ), x -X k+1 ⟩ ≥ η k+1 ⟨ĝ k+1 , X k+1 -x⟩.
We then apply Lemma B.1 with y = X k+1 and z = X k and obtain that:

η k+1 ⟨ĝ k+1 , X k+1 -x⟩ ≤ ⟨∇Φ(X k+1 ) -∇Φ(X k ), x -X k+1 ⟩ ≤ D Φ (x, X k ) -D Φ (x, X k+1 ) -D Φ (X k+1 , X k )
Using the strong convexity of D Φ stated in Inequality (3), we deduce that:

η k+1 ⟨ĝ k+1 , X k+1 -x⟩ ≤ D Φ (x, X k ) -D Φ (x, X k+1 ) -(θ k+1 -θ k ) 2 - 1 2 ∥U k+1 -U k ∥ 2 ,
which can be written as:

D Φ (x, X k+1 ) ≤D Φ (x, X k ) -(θ k+1 -θ k ) 2 - 1 2 ∥U k+1 -U k ∥ 2 -η k+1 ⟨ĝ k+1 , X k+1 -x⟩ ≤D Φ (x, X k ) -(θ k+1 -θ k ) 2 - 1 2 ∥U k+1 -U k ∥ 2 -η k+1 ⟨ĝ k+1 , X k -x⟩ -η k+1 ⟨ĝ k+1 , X k+1 -X k ⟩ (19)
Let us rst handle the last term of the right hand side:

η k+1 ⟨ĝ k+1 , X k+1 -X k ⟩ = η k+1 ĝk+1,2 (θ k+1 -θ k ) + η k+1 ⟨ĝ k+1,1 , U k+1 -U k ⟩.
The Young inequality |ab| ≤ a 2 2c + cb 2 2 leads to:

η k+1 ĝk+1,2 (θ k+1 -θ k ) ≤ η 2 k+1 ĝ2 k+1,2 4 + (θ k+1 -θ k ) 2 ≤ C α η 2 k+1 + (θ k+1 -θ k ) 2 . ( 20
)
where in the last line we use that associated with |ĝ k+1,2 | ≤ max(1, α 1-α ) (see Equation ( 6)). In the same way, we also have:

η k+1 ⟨ĝ k+1,1 , U k+1 -U k ⟩ ≤ η 2 k+1 ∥ĝ k+1,1 ∥ 2 2 + 1 2 ∥U k+1 -U k ∥ 2 . ( 21
)
We use Equations ( 20) and ( 21) in Inequality [START_REF] Cox | A theory of the term structure of interest rates[END_REF] and obtain that:

D Φ (x, X k+1 ) ≤ D Φ (x, X k ) + η 2 k+1 C α + ∥ĝ k+1,1 ∥ 2 -η k+1 ⟨ĝ k+1 , X k -x⟩. ( 22
)
We need to handle the biased drift that comes from the sampled random variables Z k+1 involved in ĝk+1 . We write:

ĝk+1 = ∇p λ (X k ) + (E[ĝ k+1 | F k ] -∇p λ (X k )) + (ĝ k+1 -E[ĝ k+1 | F k ]) = ∇p λ (X k ) -b k+1 + ∆M k+1 , with ∆M k+1 := ĝk+1 -E[ĝ k+1 | F k ]
is a martingale increment and b k+1 stands for the bias:

b k+1 = ∇p λ (U k , θ k ) -E[ĝ k+1 | F k ] = E[ Ẑk+1 | F k ] -E[Z] + λ 1-α E[Z1 ⟨Z,U k ⟩≥θ k ] -E[ Ẑk+1 1 ⟨ Ẑk+1 ,U k ⟩≥θ k | F k ] P(⟨Z, U k ⟩ ≥ θ k ) -P(⟨ Ẑk+1 , U k ⟩ ≥ θ k |F k ) . ( 23
)
Now, choosing x = x ⋆ λ and using this decomposition into Equation [START_REF] Duchi | Ecient projections onto the l 1 -ball for learning in high dimensions[END_REF], we obtain the main descent inequality, which will be the cornerstone of our analysis:

D Φ (x ⋆ λ , X k+1 ) ≤D Φ (x ⋆ λ , X k ) -η k+1 ⟨∇p λ (X k ), X k -x ⋆ λ ⟩ -η k+1 ⟨b k+1 , X k -x ⋆ λ ⟩ -η k+1 ⟨∆M k+1 , X k -x ⋆ λ ⟩ + η 2 k+1 1 4 + ∥ĝ k+1,1 ∥ 2 . ( 24
)
The Cauchy-Schwarz inequality entails:

η k+1 |⟨b k+1 , X k -x ⋆ λ ⟩| ≤ η k+1 ∥b k+1 ∥∥X k -x ⋆ λ ∥. ( 25
)
Now, the Young inequality and Equation (3) implies that:

∥X k -x ⋆ λ ∥ ≤ 1 + ∥X k -x ⋆ λ ∥ 2 2 ≤ 1 2 + 2D Φ (x ⋆ λ , X k ). ( 26
)
We then plug Equations ( 26) and ( 25) into [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF] and obtain that:

D Φ (x ⋆ λ , X k+1 ) ≤ D Φ (x ⋆ λ , X k ) 1 + 2η k+1 ∥b k+1 ∥ -η k+1 ⟨∇p λ (X k ), X k -x ⋆ λ ⟩ -η k+1 ⟨∆M k+1 , X k -x ⋆ λ ⟩ + η 2 k+1 1 4 + ∥ĝ k+1,1 ∥ 2 + 1 2 η k+1 ∥b k+1 ∥. ( 27
)
Step 2: Bias upper bound. We are led to estimate the size of ∥b k+1 ∥. For this purpose, we introduce the Kolmogorov distance d Kol , dened on real distributions and dened by:

d Kol (µ, ν) = sup a∈R µ(] -∞, a]) -ν(] -∞, a]) .
The Kolmogorov and Wasserstein distance for real distributions (see e.g. Theorem 3.1 of [START_REF] Chen | Stein's method for normal approximation[END_REF]) are related as follows:

d Kol (µ, ν) ≤ 2 W 1 (µ, ν).
Using the expression of b k+1 = (b k+1,1 , b k+1,2 ) in ( 23), its second coordinates veries :

E[∥b k+1,2 ∥] ≤ d Kol L(⟨ Ẑk+1 , U k ⟩, L(⟨Z, U k ⟩) ≤ 2 W 1 L(⟨ Ẑk+1 , U k ⟩), L(⟨Z, U k ⟩) .
Using now the variational characterization of the W 1 distance on 1-Lipschitz function, and using that for any k, U k ∈ ∆ m , we have:

∥U k ∥ 2 ≤ ∥U k ∥ 2 1 = 1,
we then deduce that:

sup u∈∆m W 1 L(⟨ Ẑk+1 , U k ⟩), L(⟨Z, U k ⟩) ≤ W 1 (L( Ẑk+1 ), L(Z)).
We then obtain from Assumption 1 Equation ( 7) that:

E[∥b k+1,2 ∥] ≤ 2 δ k+1 .
We are led to study an upper bound of ∥b k+1,1 ∥. We observe that:

E[∥b k+1,1 ∥] ≤ ∥E[ Ẑk+1 -Z]∥ + λ 1 -α E[ Ẑk+1 1 ⟨ Ẑk+1 ,U k ⟩≥θ k -Z1 ⟨Z,U k ⟩≥θ k ] ≤ δ k+1 + λ 1 -α E[ Ẑk+1 1 ⟨ Ẑk+1 ,U k ⟩≥θ k -Z1 ⟨Z,U k ⟩≥θ k ] .
We simply observe that in this case, the second term of the right hand side of the previous inequality is upper bounded by υ k+1 thanks to Assumption 1 Equation ( 8). We deduce that:

E[∥b k+1 ∥] ≤ E[∥b k+1,1 ∥] + E[∥b k+1,2 ∥] ≤ 2 δ k+1 + δ k+1 + λ 1 -α υ k+1 . ( 28
)
Step 3: Conclusion of the proof. We apply the Robbins-Siegmund Lemma to Equation ( 27): we observe that the compactness of ∆ m and the existence of a second order moment for Z leads to the existence of C > 0 such that:

E[∥ĝ k+1 ∥ 2 ] ≤ C < +∞. (29)
Finally, from Equation ( 28) and our assumptions on (δ k+1 ) k≥0 and (υ k+1 ) k≥0 , one has:

k≥0 η k+1 E[∥b k+1 ∥] < +∞. ( 30
)
From Equation ( 30),we deduce the existence of a second order moment and since η k+1 ⟨∆M k+1 , X kx ⋆ λ ⟩ is a F k Martingale-increment, the Robbins-Siegmund Lemma (see e.g. [START_REF] Robbins | A convergence theorem for non negative almost supermartingales and some applications[END_REF]) implies that:

D Φ (x ⋆ λ , X k ) a.s -→ k→∞ D ∞ ∈ L 1 and k≥0 η k+1 ⟨∇p λ (X k ), X k -x ⋆ λ ⟩ < +∞ a.s.

This ends the proof of i). ⋄

Proof of ii) This result is almost straightforward starting from i). We verify from

D Φ (x ⋆ λ , X k ) a.s -→ k→∞ D ∞ ∈ L 1 that (X k )
k≥1 is an almost surely bounded sequence. Hence, almost surely the Cesaro average sequence ( Xη k ) k≥0 is a Cauchy sequence and therefore converges towards Xη ∞ . We then use the convexity of p λ and observe that

p λ ( Xη k ) -min(p λ ) ≤ k i=0 η i k j=0 η j p λ (X k ) -min(p λ ) = k i=0 η i (p λ (X k ) -min(p λ )) k j=0 η j -→ 0 a.s.
The continuity of p λ then implies that almost surely:

Xη k -→ Xη ∞ and p λ ( Xη ∞ ) = min(p λ ). ⋄ Proof of iii)
The proof of this last point is more intricate and needs the introduction of more sophisticated tools of dynamical systems to obtain the convergence of the overall sequence (X k ) k≥0 towards a point of arg min(p λ ).

Step 1: continuous time trajectory. We follow the roadmap of [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] that introduces the asymptotic pseudo-trajectories for stochastic algorithms, and [START_REF] Mertikopoulos | On the convergence of gradient-like ows with noisy gradient input[END_REF] that adapts these tools to the mirror descent.

We introduce the mirror map Q and the Fenchel conjugate Φ * dened on ∈ ∆ m × R by: Note that the Fenchel coupling veries that for any (x, y) ∈ {∆ m × R} 2 :

Q(y) =
Φ(x) + Φ * (y) -⟨y, x⟩ = D Φ (x, Q(y)). (31)
The Mirror Descent ordinary dierential equation is then dened as:

ẏ = -∇p λ (x) x = Q(y) . (32)
From Equation (3), we know that Φ is 1/2 strongly convex, which implies rst that Q is a 2-Lipschitz continuous function (see Proposition 2.2 of [START_REF] Mertikopoulos | On the convergence of gradient-like ows with noisy gradient input[END_REF]). Second, -∇p λ is a continuous bounded function, which implies with the Cauchy-Arzela-Peano theorem that any solution of the Cauchy problem always exists. 1 We consider any x ⋆ any minimizer of p λ and (x t , y t ) a solution of [START_REF] Fujino | Large-scale multiclass support vector machine training via euclidean projection onto the simplex[END_REF]. Following the arguments of Nemirovski and Yudin [START_REF] Nemirovskij | Problem complexity and method eciency in optimization[END_REF], we dene the Fenchel coupling as:

V x ⋆ (t) = Φ(x ⋆ ) + Φ * (y t ) -⟨y t , x ⋆ ⟩.
Note that using (31), V x ⋆ (t) = D ϕ (x ⋆ , y t ). The Fenchel inequality yields V x ⋆ (t) ≥ 0 for any time t. We observe that :

V ′ x ⋆ (t) = -⟨∇p λ (x t ), x t -x ⋆ ⟩ ≤ p λ (x ⋆ ) -p λ (x t ) ≤ 0 such that V x ⋆ is non-increasing with time.
In turn, it implies the compactness of any trajectory (and the convergence with a O(t -1 ) of the value of p λ (x t ) towards min(f ) where xt is the Cesaro average of the trajectory). From any trajectory (x t ) t≥0 , we consider x ∞ any adherence point. If x ∞ does not belong to arg min(p λ ), then we can nd a compact neighbourhood V 1 of x ∞ and an increasing sequence

(t k ) k≥0 such that x t k ∈ V 1 .
Using the continuity of p λ we can nd a second set V 2 and a radius r > 0 such that

∀x ∈ V 1 : B(x, r) ⊂ V 2 ,
and for which we have:

∀x ∈ V 2 : -⟨∇p λ (x), x -x ⋆ ⟩ ≤ (min(p λ ) -p λ (x)) ≤ -a < 0.
We use that Q is 2-Lipschitz to obtain that:

∀s > 0 : ∥x t k +s -x t k ∥ = ∥Q(y t k +s ) -Q(y t k )∥ ≤ 2∥y t k +s -y t k ∥ ≤ 2 t k +s t k ∥ ẏu ∥du ≤ 2s sup x∈∆m×R ∥∇p λ (x)∥ ≤ 2s(α + E[∥Z∥]) 1 + λ 1 -α ,
where the last line comes from Equations ( 4) and ( 5). We can deduce that:

∀s ≤ s r := r 2(α + E[∥Z∥])(1 + λ 1-α ) : x t k +s ∈ B(x t k , r) ⊂ V 2 .
Finally, we write that:

V x ⋆ (t k + s r ) -V x ⋆ (0) ≤ k j=0 t j +sr t j V ′ x ⋆ (u)du ≤ -k × (as r ) -→ -∞ as k -→ +∞.
This contradicts the compactness of the trajectories and leads to a contradiction. We deduce that x ∞ ∈ arg min(p λ ) and that any adherence point of (x t ) t≥0 belongs to arg min(p λ ).

Finally, we prove that the whole trajectory converges towards x ∞ . The construction of V x ⋆ holds for any x ⋆ in arg min{p λ }, in particular for x ⋆ = x ∞ . For this particular choice, we verify that:

∀t ≥ t k 0 ≤ V x∞ (t) ≤ V x∞ (t k ),
because V is non-increasing and larger than 0.

Since by construction

x t k converges to x ⋆ , then D Φ (x ∞ , x t k ) = V x ⋆ (t k ) -→ 0 as k -→ +∞, which in turn implies that V x∞ (t) -→ 0 as t -→ +∞.
We then conclude that the whole trajectory converges towards x ∞ .

Step 2: stochastic algorithm convergence. To derive the asymptotic behaviour of the SMD, we follow the argument of Proposition 4.1 in [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] . We consider

(X k , Y k ) that satises the stochastic recursion Y k+1 = Y k -η k+1 ĝk+1 X k+1 = Q(Y k+1 )
We then dene the continuous time ane interpolation of (X k , Y k ) k≥1 using the natural timescale τ 0 = 0 and τ k+1 = τ k + η k+1 and introduce ( Xt , Ŷt ) t≥0 the stochastic process dened

as ∀t ∈ [τ k , τ k+1 ] ( Xt , Ŷt ) = (X k , Y k ) + t -τ k τ k+1 -τ k (X k+1 , Y k+1 ) -(X k , Y k ) .
We still use the decomposition derived from ( 23):

ĝk+1 = ∇p λ (X k ) + ∆M k+1 + b k+1 ,
where the martingale increment ∆M k+1 satises (see ( 29)):

sup k≥0 E[∥∆M k+1 ∥ 2 ] < +∞,
and the bias term satises under Equation [START_REF] Kloeden | Numerical solution of SDE[END_REF] and assumption of iii) that:

k≥0 E[∥b k+1 ∥] ≤ k≥0 2 δ k+1 + δ k+1 + λ 1 -α υ k+1 ≤ +∞.
Then, the Markov inequality and the Borel-Cantelli lemma yields b k -→ 0 a.s. as k -→ +∞.

We then use Remark 4.5 and Proposition 4.2 of [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] and deduce that Assumption A1 of Proposition 4.1 [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] holds. In the meantime, Assumption A2 of Proposition 4.1 [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] is valid from the proof of i) and the convergence of

D Φ (x ⋆ λ , X k ).
Since Q is a 2-Lipschitz function, we then deduce that ( Xt , Ŷt ) t≥0 is an aymptotic pseudo-trajectory of the deterministic mirror descent o.d.e. stated in [START_REF] Fujino | Large-scale multiclass support vector machine training via euclidean projection onto the simplex[END_REF]. In particular, (X t ) t≥0 converges using Step 1 of iii), which concludes the proof of iii).

B.2 Proof of the non-asymptotic behaviour (Theorem 3).

Step

1: control of E[D Φ (x ⋆ λ , X k )].
For the sake of convenience, we denote by:

E[D Φ (x ⋆ λ , X k )] = D k Φ .
Taking the expectation in [START_REF] Kiwiel | Breakpoint searching algorithms for the continuous quadratic knapsack problem[END_REF] and using ⟨∇p λ (X k ), X k -x ⋆ λ ⟩ > 0, we get:

D k+1 Φ ≤ D k Φ + 2η k+1 E D Φ (x ⋆ λ , X k )∥b k+1 ∥ + η 2 k+1 1 4 + E[∥ĝ k+1,1 ∥ 2 ] + 1 2 η k+1 E[∥b k+1 ∥] ≤ D k Φ 1 + 2η k+1 2 δ k+1 + δ k+1 + λυ k+1 1 -α + η 2 k+1 1 4 + E[∥ĝ k+1,1 ∥ 2 ] + 1 2 η k+1 2 δ k+1 + δ k+1 + λυ k+1 1 -α ,
where we used (28) and a conditional expectation argument. The inequality can be written as:

D k+1 Φ ≤ D k Φ (1 + a k+1 ) + b k+1 , using that E∥ĝ k+1,1 ∥ 2 ≤ M , a C > 0 exists such that: a k+1 = 2η k+1 2 δ k+1 + δ k+1 + λυ k+1 1-α b k+1 = C η 2 k+1 + η k+1 2 δ k+1 + δ k+1 + λυ k+1 1-α (33)
Now, denote:

V k = D k Φ k i=1 (1 + a i ) - k j=0 b j j i=1 (1 + a i )
, we verify that:

V k+1 = D k+1 Φ k+1 i=1 (1 + a i ) - k+1 j=1 b j j i=1 (1 + a i ) ≤ D k Φ (1 + a k+1 ) + b k+1 k+1 i=1 (1 + a i ) - k+1 j=1 b j j i=1 (1 + a i ) ≤ D k Φ k i=1 (1 + a i ) + b k+1 k+1 i=1 (1 + a i ) - k+1 j=1 b j j i=1 (1 + a i ) ≤ D k Φ k i=1 (1 + a i ) - k j=1 b j j i=1 (1 + a i ) = V k .
From this inequality, we deduce that:

D k Φ k i=1 (1 + a i ) - k j=1 b j j i=1 (1 + a i ) ≤ V 0 = D 0 Φ .
Let us remark that:

D 0 Φ = D Φ (X 0 , x ⋆ λ ) = (θ 0 -θ ⋆ λ ) 2 2 + D φ (U 0 , u ⋆ λ ) = (θ 0 -V@R α (u ⋆ λ )) 2 2 + D φ (U 0 , u ⋆ λ ) ≤ (θ 0 -V@R α (u ⋆ λ )) 2 2 + log m := {∆ 0 Φ } 2 , (34)
where we used the complementary slackness condition (see Equation [START_REF] Chen | Stein's method for normal approximation[END_REF] in the Appendix A) and the upper bound on the entropic Bregman divergence used in [START_REF] Lan | Validation analysis of mirror descent stochastic approximation method[END_REF]. We nally deduce that:

D k Φ ≤ D 0 Φ k i=1 (1 + a i ) + k j=1 b j j i=1 (1 + a i ) k i=1 (1 + a i ). (35)
Step 2: non-asymptotic control of the position of the SMD. Going back again to [START_REF] Kiwiel | Breakpoint searching algorithms for the continuous quadratic knapsack problem[END_REF],

we can now use the convergence of D Φ (x ⋆ λ , X k ) → 0 to get:

η k+1 ⟨∇p(X k ), X k -x ⋆ ⟩ ≤ D Φ (x ⋆ , X k ) 1 + 2η k+1 ∥b k+1 ∥ -D Φ (x ⋆ , X k+1 ) + η k+1 ⟨E[ĝ k+1 |F k ] -ĝk+1 , X k -x ⋆ ⟩ + η 2 k+1 1 4 + ∥ĝ k+1,1 ∥ 2 + 1 2 η k+1 ∥b k+1 ∥.
From the convexity of p λ , we deduce that

⟨∇p λ (X k ), X k -x ⋆ λ ⟩ ≥ p λ (X k ) -p λ (x ⋆ λ )
, and taking the expectation in both sides, we deduce that:

η k+1 E[p λ (X k )] -p λ (x ⋆ λ ) ≤ D k Φ -D k+1 Φ + a k+1 D k Φ + b k+1 .
Summing these inequalities from k = 0 to n, we get:

n k=0 η k+1 E[p λ (X k )] -p λ (x ⋆ λ ) ≤ n k=0 D k Φ -D k+1 Φ + n k=0 a k+1 D k Φ + b k+1 ≤ D 0 Φ -D n+1 Φ + n k=0 a k+1 D k Φ + b k+1 ≤ D 0 Φ + n k=0 a k+1 D k Φ + b k+1 .
Following the convexity argument of ii) Theorem 2 and the denition of Xη n in [START_REF] Bardou | Recursive computation of value-at-risk and conditional value-at-risk using mc and qmc[END_REF], we deduce that:

E[p λ ( Xη n )] -p λ (x ⋆ λ ) ≤ n k=0 η k+1 n j=0 η j+1 E[p λ (X k )] -p λ (x ⋆ λ ) ≤ 1 n j=0 η j+1 D 0 Φ + n k=0 a k+1 D k Φ + b k+1 , (36)
which ends the proof. □

Proof of Corollary 2.1. We now consider a xed step sequence for the algorithm and the discretization stopped at a horizon n. More precisely, we assume η k+1 = η > 0 for 0 ≤ k ≤ n and that [START_REF] Kloeden | Numerical solution of SDE[END_REF] translates:

E[∥b k+1 ∥] ≤ ω, ∀ 0 ≤ k ≤ n.
We aim at choosing constant values for step sequences adapted to the nite horizon n.

Within this framework the step sequences (a k+1 ) and (b k+1 ) dened in [START_REF] Markowitz | Portfolio selection[END_REF] become

a k+1 = 2ηω and b k+1 = C(η 2 + ηω).
Let us now see how both η and ω should depend on n. With these choices, [START_REF] Nemirovskij | Problem complexity and method eciency in optimization[END_REF] becomes:

D k Φ ≤ D 0 Φ (1 + 2ηω) k + C(η 2 + ηω) (1 + 2ηω) k -1 2ηω .
Equation ( 36) then yields :

E[p λ ( Xη n )] -p λ (x ⋆ λ ) ≤ 1 (n + 1)η D 0 Φ + n k=0 2ηωD k Φ + Cη(η + ω) ≤ 1 (n + 1)η D 0 Φ 1 + 2ηω n k=0 (1 + 2ηω) k + Cη(η + ω) n + 1 + n k=0 (1 + 2ηω) k+1 -1 ≤ D 0 Φ (n + 1)η 1 + 2ηω (1 + 2ηω) n+1 -1 2ηω + Cη(η + ω) η(n + 1) (1 + 2ηω) n+1 -1 2ηω + n + 1 ≤ D 0 Φ e 2ηω(n+1) (n + 1)η + C(η + ω) e 2ηω(n+1) 2ηω(n + 1) + 1 (37)
where we used in the last inequality that 1 + x ≤ e x . We then choose (η, ω) such that:

η = D 0 Φ 2 √ n + 1 and ω = 1 √ n + 1 D 0 Φ
and obtain, using the fact that ηω(n + 1) ≤ 1/2, that a large enough C exists such that

E[p λ ( Xη n )] -p λ (x ⋆ λ ) ≤ C D 0 Φ √ n + 1 .
To conclude, we plug-in the upper bound [START_REF] Mertikopoulos | On the convergence of gradient-like ows with noisy gradient input[END_REF] for D 0 Φ in the previous equations. Starting from (37), we also observe that the choice η = ω = 1/ √ n + 1 yields:

E[p λ ( Xη n )] -p λ (x ⋆ λ ) ≤ C D 0 Φ √ n + 1 .
C Proofs associated with the simulation of the portfolio

The main goal of this section is to prove Proposition 3.1 and Proposition 3.2. To this aim, we rst need technical properties on the drift-implicit Euler Scheme coupled with the Riemann integral approximation and on the density of the portfolio Z.

C.1 Properties of the CIR process and its discretization Conditional moments of the CIR process We recall some useful properties satised by the CIR process and on some important objects that are related to the CIR introduced in Equation ( 10).

Proposition C.1. (r t ) t≥0 is a Markov process and satises:

i) If 2ab > σ 2 0 , then (r t ) t≥0 remains strictly positive almost surely. ii) The (conditional) expectation and variance are given by:

E[r t+s |r s ] = r s e -at + b(1 -e -at ), (38)
and

V(r t+s |r s ) = r s σ 2 0 a (e -at -e -2at ) + bσ 2 0 2a (1 -e -at ) 2 . ( 39
)
iii) The CIR possesses an exponential integrability:

∀λ < a 2 2σ 2 0 ∀t ≥ 0 : E exp λ t 0 r s ds ≤ C(λ, t) < +∞.
The two rst items of the previous proposition are classical, and may be found in [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF]. The exponential integrability of the CIR process is more recent and may be traced back to [START_REF] Cozma | Exponential integrability properties of Euler dicretization schemes for the cox-ingersoll-ross process[END_REF] (Proposition 3.2). We emphasize that within realistic situations, the values of a and σ for the CIR process are in general linked within a reasonable ratio of 5 ≤ a σ ≤ 10, so that the limiting support of λ is generally signicantly larger than 1. Note also that this last result holds at any horizon time but it may be explicited for t = 1 with a larger size of λ.

Proposition C.1 will be useful for deriving some important properties related to the bias of our approximation strategy.

Weak error associated to the implicit Euler scheme We recall the following statement of [START_REF] Alfonsi | Strong order one convergence of a drift implicit euler scheme: Application to the cir process[END_REF] proved in [START_REF] Dereich | An Euler-type method for the strong approximation of the cox ingersoll ross process[END_REF], that provides an upper bound of the weak error associated to the implicit Euler scheme [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF]. In the previous reference, a strong error rate of order 1 is also provided with stronger assumptions on the parameters. Proposition C.2 (Theorem 1 in [START_REF] Alfonsi | Strong order one convergence of a drift implicit euler scheme: Application to the cir process[END_REF]). Assume that 2ab > σ 2 0 , then for any p ∈ [1, 2ab σ 2 ), a constant K p exists such that:

E max k∈{0,...,N } rkh -r kh p 1/p ≤ K p √ h.
Control of the error in the approximation of the integral of the CIR Recall the denition of ∆ h in (47)

∆ h = 1 0 r s ds - 1 N N -1 k=0 rkh .
Let us now prove two technical lemmas rst on the estimation of E(∆ 2 h ) then on exponential moments.

Lemma C.1. Assume that the CIR parameters satisfy ab > σ 2 , then there exists a constant C such that, for any h = N -1 > 0,

E[|∆ h | 2 ] ≤ Ch.
Proof of Lemma C.1. We bound ∆ 2 h using rst the Jensen inequality and then the Cauchy-Schwarz inequality, we obtain that:

∆ 2 h = 1 N N N -1 k=0 (k+1)h kh (r s -rkh )ds 2 ≤ 1 N N k=1 N (k+1)h kh (r s -rkh )ds 2 ≤ N k=1 (k+1)h kh (r s -rkh ) 2 ds.
Now taking the expectation and decomposing between the approximation error from the numerical scheme at time kh and the regularity of the CIR process we obtain that:

E[∆ 2 h ] ≤ 2 N k=1 kh (k-1)h E[(r s -r kh ) 2 ] + E[(r kh -rkh ) 2 ] ds ≤ 2 N k=1 kh (k-1)h E[(r s -r kh ) 2 ]ds + 2h N -1 k=0 E[(r kh -rkh ) 2 ]. (40)
Let us start with the rst term in the right hand side of the equation above. We use the explicit expression of conditional expectation and variances of the CIR r t stated in ii) of Proposition C.1. Thus, for any s = kh + u with 0 ≤ u ≤ h, the bias-variance decomposition associated with Equations ( 38) and ( 39) leads to:

E (r s -r kh ) 2 |r kh = E[r s |r kh ] -r kh 2 + V[r s |r kh ] = (b -r kh ) 2 (1 -e -au ) 2 + σ 2 a r kh e -as (1 -e -au ) + b σ 2 2a (1 -e -au ) 2 .
We then use again Equation [START_REF] Ch | Some remarks on the value-at-risk and the conditional value-at-risk[END_REF] to obtain a recursive expression of E(r kh ):

E[r (k+1)h ] = E E[r (k+1)h |r kh ] = E[r kh )e -ah + b(1 -e -ah ] = e -a(k+1)h (E[r 0 ] -b) + b. Therefore N k=1 E[r kh ] = bN + (E[r 0 ] -b) N k=1 e -akh = bN + (E[r 0 ] -b) 1 -e -a 1 -e -ah .
From Equation ( 39) and a conditional expectation argument we deduce that:

E[(b -r (k+1)h ) 2 ] = E[r kh ] σ 2 a (e -ah -e -2ah ) + bσ 2 2a (1 -e -ah ) 2 .
We then obtain that:

N k=1 kh (k-1)h E[(r s -r kh ) 2 ]ds = N k=1 E[(b -r kh ) 2 ] h 0 (1 -e -au ) 2 du + σ 2 a N k=1 E[r kh ] h 0 e -au (1 -e -au )du + b σ 2 2a N h 0 (1 -e -au ) 2 du.
We then compute that:

h 0 (1 -e -au ) 2 du = a 2 h 3 3 + O a (h 4 ) and h 0 e -au (1 -e -au )du = ah 2 2 - a 2 h 3 2 + O a (h 4 ).
Finally we deduce that for some C > 0

N k=1 kh (k-1)h E[(r s -r kh ) 2 ] = Ch + O a,b (h 2 ). (41)
We now turn to the second term above, that we estimate using the strong error rate of Proposition C.2 for p = 2. Assuming ab > σ 2 0 , we get:

E[(r kh -rkh ) 2 ] ≤ E[ max 1≤k≤N (r kh -rkh ) 2 ] ≤ Ch.
Since this bound is uniform in k, we deduce that :

N k=1 E[(r kh -rkh ) 2 ] ≤ C. (42)
Combining estimates ( 41) and [START_REF] Rockafellar | Random variables, monotone relations, and convex analysis[END_REF] in equation [START_REF] Robbins | A stochastic approximation method[END_REF] gives the conclusion.

Lemma C.2. For any choice of q ≤ 2, if a > 2 √ 2σ 0 then:

sup N E[e q|∆ h | ] < +∞. (43)
Proof of Lemma C.2. We use the Cauchy-Schwarz inequality to upper-bound the expectation:

E[e q|∆ h | ] = E e q 1 0 rsds-1 N N k=1 rkh
≤ E e q 1 0 rsds e q 1 N N k=1 rkh

≤ E e 2q 1 0 rsds

1 2 × E e 2q N N k=1 rkh 1 2 .
We apply iii) of Proposition C.1 and observe that since a > 2 √ 2σ 0 , then 2q ≤ 4 ≤ a 2 2σ 2 0 and:

E e 2q 1 0 rsds 1 2 < +∞.
For the second part, we rely on the recursive expression for rkh stated in (C.2), that we bound using:

∀ϵ > 0 (u + v) 2 = u 2 + v 2 + 2uv ≤ u 2 (1 + 1/ϵ) + v 2 (1 + ϵ). Thus r(k+1)h = √ rkh + σ 0 2 ∆B (k) 0 2(1 + ah 2 ) + √ rkh + σ 0 2 ∆B (k) 0 2 4(1 + ah 2 ) 2 + (4ab -σ 2 0 )h 8(1 + ah 2 ) 2 ≤ (1 + ϵ 2 ) √ rkh + σ 0 2 ∆B (k) 0 2(1 + ah 2 ) 2 + (1 + 1 2ϵ ) √ rkh + σ 0 2 ∆B (k) 0 2 4(1 + ah 2 ) 2 + (4ab -σ 2 0 )h 8(1 + ah 2 ) ≤ (2 + ϵ 2 + 1 2ϵ ) √ rkh + σ 0 2 ∆B (k) 0 2 4(1 + ah 2 ) 2 + (1 + 1 2ϵ ) (4ab -σ 2 0 )h 8(1 + ah 2 ) ≤ (2 + ϵ 2 + 1 2ϵ )(1 + 2τ ) 4(1 + ah 2 ) 2 rkh + (2 + ϵ 2 + 1 2ϵ )(1 + 1 2τ ) 4(1 + ah 2 ) 2 σ 2 0 4 {∆B (k) 0 } 2 + (1 + 1 2ϵ ) (4ab -σ 2 0 )h 8(1 + ah 2 )
.

Choosing ϵ = 1 and τ = 1 7 < 1 6 , we deduce that:

r(k+1)h ≤ 1 (1 + ah 2 ) 2 rkh + 27 32(1 + ah 2 ) 2 σ 2 0 4 {∆B (k) 0 } 2 + 3(4ab -σ 2 0 )h 16(1 + ah
2 ) Using the independence and stationarity of the Brownian motion, we set ξ k+1 a standard Gaussian random variable such that ξ k+1 √ h = ∆B .

A straightforward recursion yields:

r(k+1)h ≤ r0 + α k j=0 ξ 2 j+1 + βk.
In particular, we can use this to upper bound the sum with:

1 N N k=1 rkh ≤ r0 + α N j=1 ξ 2 j+1 + βN.
We can now bound the Laplace transform, using the independence of the ξ j and the Laplace transform of a χ

2 distribution E(e λξ 2 ) = (1 -2λ) -1/2 , ∀λ < 1/2. E e 2q N N k=1 rkh ≤ e 2qr 0 +2qN β N j=1 E e 2qαξ 2 j ≤ e 2qr 0 +2qN β 1 1 -4qαN 1/2 ≤ e 2qr 0 +2qN β e 2qαN ,
where the last line derives from 1/(1 -t) ≤ e t . From our choice of α and β, we observe that αN < +∞ and βN < +∞. We then conclude that E e 2q Proposition C.3 (Density of the portfolio). Assume that Σ is positive denite, then the distribution of S = (Y 1 , S 1 , . . . , S m ′ ) at time 1 is absolutely continuous with respect to the Lebesgue measure. Moreover, the density p may be written as:

p : (z 1 , . . . , z m ) -→ M (z 1 , . . . , z m ) (σ 1 . . . σ m ′ ) det(L)z 1 . . . z m ,
where M is a bounded continuous function, whose bound is independent from m and Σ.

Proof. Let us argue rst in the case of uncorrelated Brownian motions in equations [START_REF] Bardou | Computing var and cvar using stochastic approximation and adaptive unconstrained importance sampling[END_REF], i.e.

when Σ is the identity matrix. The correlated case will again be dealt with using Cholesky decomposition. Note that

Y 1 (t) = Y 1 (0) exp( t 0 r s ).
The density of the integral of the CIR process can be derived from Pitman and Yor [START_REF] Pitman | A decomposition of bessel bridges[END_REF], wherein the authors give an expression of the Laplace transform of the integral of a Bessel process. This representation can in turn be used to derive the Laplace transform of the Y 1 (t), which upon inversion, gives its density, see e.g. Gulisashvili and Stein [START_REF] Gulisashvili | Asymptotic behavior of the stock price distribution density and implied volatility in stochastic volatility models[END_REF]. Next, since the assets S 1 (t), . . . , S m ′ (t) are geometric brownian motions, their density can be computed explicitly, and we obtain the density of the vector (Y 1 , S 1 , . . . , S m ′ ) by multiplying these densities together, thanks to the independence of the driving Brownian motions. Furthermore, to justify the existence of a density for the pair ( 1 0 r s ds, W 0 (1)), we can once again cite [START_REF] Gulisashvili | Asymptotic behavior of the stock price distribution density and implied volatility in stochastic volatility models[END_REF]. This density is denoted by q and veries for any bounded function Ψ:

E Ψ( 1 0 r s ds, W 0 (1)) = Ψ(I, w 0 )q(I, w 0 )dIdw 0 .
Below, we denote by q 1 the marginal distribution of 1 0 r s ds, whose density is bounded on R (see among others e.g. [START_REF] Pitman | A decomposition of bessel bridges[END_REF] and [START_REF] Gulisashvili | Asymptotic behavior of the stock price distribution density and implied volatility in stochastic volatility models[END_REF]). To describe the other components, we introduce:

c k = µ k - {σ k } 2 2 ,
and using a conditional expectation argument, we have:

E[f (Y 1 , S 1 . . . , S m ′ )] = E f Y 1 (0)e I , S 1 (0)e c 1 +σ 1 (( L W ) 1 +ℓ 0,1 W 0 ) , . . . , S m ′ (0)e c m ′ +σ m ′ (( L W ) m ′ +ℓ 0,m ′ W 0 ) = R + ×R R m-1 q(I, w 0 )γ m ′ ( w) f Y 1 (0)e I , S 1 (0)e c 1 +σ 1 ( L w) 1 +ℓ 0,1 w 0 , S m ′ (0)e c m ′ +σ m ′ ( L w) m ′ +ℓ 0,m ′ w 0 d wdIdw 0 ,
where γ m-1 refers to the m -1 dimensional standard centered Gaussian distribution. For a xed value of w 0 , we consider the change of variables:

z 1 = Y 1 (0)e I and z i+1 = S i (0)e c i +σ i ( L w) i +σ i ℓ 0,i w 0 , which is equivalent to I(z 1 ) = log( z 1 Y 1 (0) ) and w(z 2 , . . . , z m , w 0 ) = L-1      {σ 1 } -1 log z 2 S 1 (0) -c 1 -ℓ 0,1 w 0 . . . {σ m ′ } -1 log zm S m ′ (0) -c m ′ -ℓ 0,m ′ w 0      .
Since L is invertible, the Jacobian of the change of variable is then given by:

dId w = 1 (σ 1 . . . σ m ′ ) det L dz 1 dz 2 . . . dz m z 1 . . . z m .
We observe that det L = det L, which leads to:

E[f (Y 1 , S 1 . . . , S m ′ )] = {R + } m f (z 1 , . . . , z m ) R q log z 1 Y 1 (0) , w 0 γ m-1 ( w(z 2 , . . . , z m , w 0 ))dw 0 (σ 1 . . . σ m ′ ) det L × (z 1 . . . z m ) dz 1 . . . dz m .
The previous expression then guarantees that (Y 1 , S 1 , . . . , S m ′ ) has a distribution uniformly continuous with respect to the Lebesgue measure on {R + } m , whose density is of the form:

(z 1 , . . . , z m ) -→ M (z 1 , . . . , z m ) (σ 1 . . . σ m ′ ) det(L)z 1 . . . z m ,
where

M (z 1 , . . . , z m ) = R q log z 1 Y 1 (0)
, w 0 γ m-1 ( w(z 2 , . . . , z m , w 0 ))dw 0 .

The Gaussian density γ m-1 being continuous and bounded by (2π) -m ′ /2 , we then have:

M (z 1 , . . . , z m ) ≤ R q log z 1 Y 1 (0) , w 0 dw 0 = q 1 (I(z 1 )) ≤ ∥q 1 ∥ ∞ < +∞.
Therefore, M is a bounded continuous function, whose bound is independent from m and Σ.

Below, we will need to upper bound the probability of sliced events that are dened as follows. Consider u = (u 1 , v) a vector of the simplex S m and ϵ > 0, we introduce:

Ω ϵ (θ, u) = θ ≤ ⟨Z, u⟩ ≤ θ + ϵ .
From Proposition C.3, we deduce the following property.

Proposition C.4 (Sliced events). Consider a portfolio Z dened by Equation [START_REF] Bardou | Computing var and cvar using stochastic approximation and adaptive unconstrained importance sampling[END_REF]. Assume σ + = sup 1≤i≤m σ i < +∞, and that the correlation matrix Σ is invertible. Assume that the initialization is given by Z i (0) = 1. A constant K exists such that for any ϵ > 0, for any ρ ∈ (0, 1):

P[Ω ε (θ, u)] ≤ Kme {σ + } 2 m 2 (1-ρ) 4ρ 2 ϵ 1-ρ 2 , ∀θ ∈ R, ∀u ∈ S m .
Proof. We consider η > 0 small enough (whose dependency with ϵ and m will be precised later on) and we observe that:

P[Ω ϵ (θ, u)] = E 1 Ωϵ(θ,u) = E 1 Ωϵ(θ,u) 1 min Z i ≥η + 1 Ωϵ(θ,u) 1 ∃i :Z i ≤η ≤ E 1 Ωϵ(θ,u) 1 min Z i ≥η + m i=1 P[Z i ≤ η]. (44)
Step 1: Let us rst consider the second term in Equation [START_REF] Rockafellar | Conditional value-at-risk for general loss distributions[END_REF]. We then consider two separate cases.

If i = 1 and considering η ≤ Z 1 (0), we then observe that Z 1 ≥ Z 1 (0) since r is a nonnegative process, so that P[Z 1 ≤ η] = 0.

Assume now that i ′ = i -1 ∈ {1, . . . , m ′ }, using the notations of Proposition C.3, we observe that:

Z i ′ ≤ η ⇐⇒ Z i ′ (0)e c i ′ +σ i ′ B i ′ ≤ η ⇐⇒ B i ′ ≤ -{σ i ′ } -1 log 1 η + c i ′ , so that P[Z i ′ ≤ η] = P B i ′ ≤ -{σ i ′ } -1 log 1 η + c i ′ = P N (0, 1) ≥ {σ i ′ } -1 log 1 η + c i ′ ≤ σ i ′ e - log 1 η +c i ′ 2 2{σ i ′ } 2 √ 2π log 1 η + c i ′
, where the last line comes from standard estimation of the Gaussian tail. Choosing η small enough, we then deduce that a constant K exists such that:

∀i ∈ {2, . . . , m}

P[Z i ≤ η] ≤ Ke - (log 1 η ) 2 4{σ i } 2 . (45)
Dening σ + = sup 1≤i≤m σ i , a union bound then leads to:

m i=1 P[Z i ≤ η] ≤ Kme - (log 1 η ) 2 4{σ + } 2 = Km(η log 1 η 4{σ + } 2 ).
Step 2: We study the rst term of (44) and for any vector q = (q 1 , . . . , q m ) and for any integer ℓ, we denote by q -ℓ the vector q = (q 1 , . . . , q ℓ-1 , q ℓ+1 , . . . , q m ). For any u ∈ S m , we can nd ℓ such that u ℓ ≥ 1 m , so that: We use the Cauchy-Schwarz inequality in order to bound the inner integral: where we use the fact that we integrate over [η, +∞[ and that the function M is bounded. We nally obtain that for a positive constant K:

E[1 Ωϵ(θ,u) 1 min Z i ≥η ] = E[1 θ-⟨Z -ℓ ,u -ℓ ⟩≤Z ℓ u ℓ ≤θ+ϵ-⟨Z -ℓ ,u -ℓ ⟩ 1 min Z i ≥η ] ≤ E[1 (u ℓ ) -1 [θ-⟨Z -ℓ ,u -ℓ ⟩]≤Z ℓ ≤(u ℓ ) -1 [θ-⟨Z -ℓ ,u -ℓ ⟩]+(u ℓ ) -1 ϵ 1 min Z i ≥η ].
E[1 Ωϵ(θ,u) 1 min Z i ≥η ] ≤ √ mϵ K η m/2 . ( 46
)
Step 3: Final bound We now choose η to balance the size of Equation [START_REF] Sharpe | Mutual fund performance[END_REF] and Equation ( 46) in the decomposition [START_REF] Rockafellar | Conditional value-at-risk for general loss distributions[END_REF]. Fix ρ ∈ (0, 1), we choose η such that η = ϵ ρ/m , so that:

(46) ≤ K √ m η -m ϵ = K √ mϵ 1-ρ 2 .
We also verify with this choice of η that: .

To balance both contributions, we aim at nding the smallest constant C such that:

- m 2 4{σ + } 2 ρ 2 (log(ϵ)) 2 ≤ 1 -ρ 2 log(ϵ) + C, which is equivalent to m 2 4{σ + } 2 ρ 2 (log(ϵ)) 2 + 1 -ρ 2 log(ϵ) ≥ -C.
Now we consider the left hand side as a polynomial of degree 2 in log(ϵ) and obtain that the inequality is satised for any ϵ when the constant C is chosen as:

C = m 2 {σ + } 2 (1 -ρ) 4ρ 2 .
Finally we obtain that:

≤ Kme {σ + } 2 m 2 (1-ρ) 4ρ 2 ϵ 1-ρ 2 ,
which concludes the proof.

Using |1 -e x | ≤ |x|e |x| and a three way Holder's inequality, we obtain that for p, q and r such that 1 p + 1 q + 1 r = 1 then:

E(| Ŷ h 1 -Y 1 |) ≤ E |Y 1 | r 1/r E |∆ h | p 1/p E e q∆ h 1/q .
Using Lemma C.1 and C.2, we deduce that for any choice of r, q such that 1 q + 1 r = 1/2, we have:

E[| Ŷ h 1 -Y 1 |] ≤ C q E[|Y 1 | r ] 1/r √ h, (48)
for C q large enough (independent from h), which concludes the proof. Taking the expectation and applying the Cauchy-Schwarz inequality in the second term we obtain that B can be bounded by:

B ≤ √ m E[(|Y | + S 2 ) 2 ] 1/2 E |1 ⟨Z,w⟩≥θ -1 ⟨ Ẑ,w⟩≥θ | 1/2 + √ mE ∥Y -Ŷ ∥ . (49)
We apply Proposition 3.1 and deduce that the second term is of the order √ h. We then focus on the dierence of indicator functions. Let us remark that: These two events are handled similarly, we write:

Y -Ŷ = Y 1 (1 -e ∆ h ),
where ∆ h is given by [START_REF] Zhou | Stochastic mirror descent in variationally coherent optimization problems[END_REF]. Let us x ϵ > 0 and consider the event Ω ε = {u|Y -Ŷ | ≤ ε}, then

P(A 1 ) = P(A 1 ∩ Ω ε ) + P(A 1 ∩ Ω c ε ),
and now, observe that:

A 1 ∩ Ω ε ⊂ {θ ≤ ⟨Z, w⟩ ≤ θ + ε},
which lead to control the probability the the portfolio Z belongs to a slice. Consider w = (u, v) a vector of the simplex ∆ m and ϵ > 0, we introduce:

Ω ϵ (θ, w) = θ ≤ ⟨Z, w⟩ ≤ θ + ϵ .

An upper bound for this probability was obtained in Proposition C.4. A constant K exists such that for any ϵ > 0, for any ρ ∈ (0, 1):

P[Ω ε (θ, u)] ≤ Kme {σ + } 2 m 2 (1-ρ) 4ρ 2 ϵ 1-ρ 2 , ∀θ ∈ R, ∀ρ ∈ (0, 1) ∀u ∈ ∆ m ′ .
We nally study A 1 ∩ Ω c ε and write that for any ζ > 0, and using |1 -e x | ≤ |x|e |x| :

P(A 1 ∩ Ω c ε ) ≤ P[uY (1 -e ĥ) ≥ ε] ≤ P[Y (1 -e ∆ h ) ≥ ε] ≤ P[Y ≥ ζ] + P[|1 -e ∆ h | ≥ ζ -1 ε] ≤ P[Y ≥ ζ] + P[|∆ h |e |∆ h | ≥ ζ -1 ε] ≤ E[Y ] ζ + ζ E[|∆ h |e |∆ h | ] ε ≤ E[Y r ] ζ r + ζ E[∆ 2 h ] 1/2 E[e 2|∆ h | ] 1/2 ε ,
where we applied rst a union bound and then the Markov and the Cauchy-Schwarz inequalities.

We nally obtain from Lemma C.2 and C.1 that:

P(A 1 ∩ Ω c ε ) ≤ E[Y r ] ζ r + Cζ √ h ϵ . (50)
Combining Proposition C.4 and (50), we deduce that:

E 1 ⟨Z,w⟩≥θ -1 ⟨ Ẑ,w⟩≥θ ≤ C Kme {σ + } 2 m 2 (1-ρ) 4ρ 2 ϵ 1-ρ 2 + E[Y r ] ζ r + Cζ √ h ϵ .
Now using this bound as well as (48) in (49) we obtain that:

B ≤ K √ me {σ + } 2 m 2 (1-ρ) 4ρ 2 ϵ 1-ρ 2 + 1 ζ r + ζ √ h ϵ + √ h , (51)
where the parameters r > 1, ε > 0, ζ > 0 and ρ ∈ (0, 1) have to be chosen. Note that the constant depends heavily on the number of assets m and on the choice of ρ. We optimize only on the two parameters ε > 0 and ζ > 0 the quantity:

ϵ 1-ρ 2 + 1 ζ r + ζ √ h ϵ + √ h.
It is immediate to verify that ζ should be chosen such that

ζ -r = ζ √ h ε =⇒ ζ = ε 1 1+r √ h 1 1+r
.

In the meantime, the optimal value of ε satises:

ε 1-ρ 2 = ζ √ h ε =⇒ ε 3-ρ 2 = ε 1 1+r √ h 1-1 1+r .
Some straightforward computations show that ε = h 2r 1-ρ+r (3-ρ) , which in turn implies that a constant K large enough exists such that:

B ≤ K r,ρ √ me {σ + } 2 m 2 (1-ρ) 4ρ 2 h 1 2( 3-ρ 1-ρ + 1 r ) .
It is then easy to see that the rate may be arbitrarily close to h 1/6 by picking r large enough and ρ close to 0. More precisely, for any e > 0, we can nd ρ and r such that:

B ≤ K e √ me {σ + } 2 m 2 4e 2 h 1 6 -e .

  Data: Parameters of the CIR: a, b, σ 0 and of the G.B.M. (µ i , σ i ); Correlation matrix Σ Result: A sample Ẑ1 = ( Ŷ (h) 1 , S 1 1 , . . . , S m ′ 1 ) Compute the Cholesky factorization of Σ (see Equation (12)): L and L . Choose r0 .

  approximation given in Equation (15); Simulate the G.B.M. S 1 1 , . . . , S m ′ 1 with Equation (13).

Figure 1 :

 1 Figure 1: Time evolution of the return of the discretized trajectories associated to the assets of the synthetic portfolio (CIR + GBM).

Figure 2 :

 2 Figure 2: Evolution of the weights of the assets with the number of iterations of the SMD (left) and of the PSGD (right) on p λ with λ = 0.9 and α = 0.05. The composition of optimal portfolio may be dierent as represented above with SMD and PSGD.

Figure 3 :

 3 Figure3: Evolution of the CV@R and ER of the optimal portfolio with the number of iterations of the SMD (orange) and of the PSGD (blue) on p λ with λ = 0.9 and α = 0.05.

n = 10 4 n = 10

 10

Figure 4 :

 4 Figure 4: Evolution of the V@R estimation (left) of the portfolio with the computational time of the SMD and of the MC-MD p λ with λ = 0.95 and α = 0.05. The composition of optimal portfolio (right) converges to similar values as represented on the right with SMD and MC-MD but with a dierent cost in computational time.

Figure 5 :

 5 Figure 5: Composition of the optimal portfolio when λ increases. Left: 3 assets, Right: 8 assets. ER and CV@R are indicated in the legend box of each subgure.

Figure 6 :

 6 Figure 6: Sample of the ˆFVX series and of the assets that composed our portfolio, from the Yahoo! database, bewtween the 01-01-2014 and 01-01-2018 (2016 for theˆFVX time series).

Figure 7 :

 7 Figure 7: Markowitz Ecient frontier. Left: α = 5% and Right: α = 1%. Red point: optimal Sharpe ratio location. Blue points: location of the assets involved in the portfolio. Pink: risk free asset.

  arg max x∈∆m×R ⟨y, x⟩ -Φ(x) and Φ * (y) = max x∈∆m×R ⟨y, x⟩ -Φ(x).

E 1 1 ∞ η 1 (

 111 Now, we recall the joint density p(z 1 , . . . z m ) dened in Proposition C.3 and rewrite the previous inequality as:Ωϵ(θ,u) 1 min Z i ≥η ≤ [η,+∞[ m 1 (u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]≤z ℓ ≤(u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]+(u ℓ ) -1 ϵ p(z 1 , . . . , z m )dz 1 . . . dz m ≤ [η,+∞[ m-u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]≤z ℓ ≤(u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]+(u ℓ ) -1 ϵ p(z 1 , . . . , z m )dz ℓ .

∞ η 1 (E[ 1 E[ 1 M

 111 u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]≤z ℓ ≤(u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]+(u ℓ ) -1 ϵ p(z 1 , . . . , z m )dz ℓ ≤ ∞ η 1 (u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]≤z ℓ ≤(u ℓ ) -1 [θ-⟨z -ℓ ,u -ℓ ⟩]+(u ℓ ) -1 ϵ dz ℓ ∞ η p(z 1 , . . . , z m ) 2 dz ℓ ≤ √ mϵ ∞ η p(z 1 , . . . , z m ) 2 dz ℓ ,where the last line derives from u -1 ℓ ϵ ≤ mϵ. The previous inequality combined with the Jensen inequality |g|dµ ≤ g 2 dµ and the Tonelli theorem yield:Ωϵ(θ,u) 1 min Z i ≥η ] ≤ √ mϵ [η,+∞[ m-1 ∞ η p(z 1 , . . . , z m ) 2 dz ℓ i̸ =ℓ dz i ≤ √ mϵ [η,+∞[ m p(z 1 , . . . , z m ) 2 dz 1 . . . dz m .The expression of p obtained in Proposition C.3 yields: Ωϵ(θ,u)1 min Z i ≥η ] ≤ √ mϵ [η,+∞[ m M (z 1 , . . . , z m ) 2 (σ 1 . . . σ m ′ ) 2 det(L) 2 (z 1 . . . z m ) 2 ,dz 1 . . . dz m (z 1 , . . . , z m ) (σ 1 . . . σ m ′ ) det(L)(z 1 . . . z m ) , dz 1 . . . dz m ,

2 4{σ + } 2 ρ 2

 22 (log(ϵ))2

C. 3 ∈Lip 1 {E(f ( Ŷ h 1 ) 1 0

 3111 Proof of Proposition 3.1Proof. Let us rst rewrite the Wasserstein distance as:W 1 (L( Ŷ h 1 ), L(Y 1 )) = sup f -E(f (Y 1 ))}.Since the G.B.M. are exactly simulated and the test functionf is chosen among 1-Lipschitz functions it is sucient to bound E[| Ŷ h 1 -Y 1 |].Recall the denition of the approximation Ŷ h in Equation (15), we obtainE[|Y 1 -Ŷ h 1 |] = E Y 0 e

C. 4 Z = (Y 0 e 1 0 2 ≤ Y 1 2 ≤

 41212 Upper bound of the bias termProof of Proposition 3.2. Below, we alleviate the notations and denote by S i = S i (1), which is the return of asset i at time 1. From the exact simulation of the G.B.M. S = (S 1 , . . . , S m ′ ), we have that Z = (Y, S) whereas Ẑ = Ẑh = ( Ŷ h , S) where Ŷ h is a biased approximation of Y 1 dened in[START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent ows, with applications[END_REF]. We denote by w = (u, v) ∈ R + × (R + ) m the repartition of the portfolio and recall that w 1 = u + v 1 = 1. The goal of this section is to control the bias in terms of a function of the discretization parameter h:B = E[⟨Z, w⟩1 ⟨Z,w⟩≥θ -⟨ Ẑ, w⟩1 ⟨ Ẑ,w⟩≥θ ] 2 = E[⟨Z1 ⟨Z,w⟩≥θ -Ẑ1 ⟨ Ẑ,w⟩≥θ , w⟩] 2 .From the re-investment condition, we know that w 1 = 1 and the Cauchy-Schwarz inequality yields:B ≤ √ mE[ Z1 ⟨Z,w⟩≥θ -Ẑ1 ⟨ Ẑ,w⟩≥θ 2 ].Since we discretize only the rst component, we have that: rsds , S 1 , . . . , S m ′ ), and Ẑ = (Y 0 e 1 N N k=1 rkh , S 1 , . . . , S m ′ ),We can rewrite the dierence as:Z1 ⟨Z,w⟩≥θ -Ẑ1 ⟨ Ẑ,w⟩≥θ ⟨Z,w⟩≥θ -Ŷ 1 ⟨ Ẑ,w⟩≥θ 2 + S(1 ⟨Z,w⟩≥θ -1 ⟨ Ẑ,w⟩≥θ ) Y (1 ⟨Z,w⟩≥θ -1 ⟨ Ẑ,w⟩≥θ ) 2 + (Y -Ŷ )1 ⟨ Ẑ,w⟩≥θ 2 + S(1 ⟨Z,w⟩≥θ -1 ⟨ Ẑ,w⟩≥θ ) 2 ≤ |Y -Ŷ | + (|Y | + S 2 )|1 ⟨Z,w⟩≥θ -1 ⟨ Ẑ,w⟩≥θ |.

  ⟨Z, w⟩ = Y u + ⟨S, v⟩ and ⟨ Ẑ, w⟩ = Ŷ u + ⟨S, v⟩,thus E 1 ⟨Z,w⟩≥θ -1 ⟨ Ẑ,w⟩≥θ = E 1 ⟨S,v⟩≥θ-Y u -1 ⟨S,v⟩≥θ-Ŷ u = P(A 1 ) + P(A 2 ),whereA 1 = {θ ≤ ⟨Z, w⟩ ≤ θ + u(Y -Ŷ )} and A 2 = {θ + u(Y -Ŷ ) ≤ ⟨Z, w⟩ ≤ θ}.

  Step-size sequence (η n ) n∈N and U 0 ∈ R, θ 0 ∈ R; α ∈ (0, 1) Result: Two sequences: X k = (U k , θ k ) k≥0 for k = 0, . . . , do Simulate the random Z k+1 satisfying (7) and (8); Compute a stochastic approximation ĝk+1 of ∇p λ (U k , θ k ) with:

	Data:
	Under a suitable assumption on the biased sequences (δ k+1 ) k≥0 and (υ k+1 ) k≥0 ,
	we derive the next result.
	Theorem 2 (Almost sure convergence of the biased SMD). Assume that
	η k+1 = +∞, and	η 2 k+1 < +∞,
	k≥0	k≥0
	and that the bias sequences (δ k+1 ) k≥0 and (υ k+1 ) k≥0 introduced in Assumption 1 satisfy
	η k+1 ( δ k+1 + υ k+1 ) < +∞,
	k≥0
	then:

  and 3.2.

SMD at xed time horizon n According to our previous results, we can now consider the choice of parameters induced by Corollary 2.1. Recall that we assume a constant step-size sequence η and discretizaton step-size h. Corollary 2.1 combined with Propositions 3.1 and 3.2 induce that h should be chosen as:

It is not clear that ∇p λ is a Lipschitz continuous function, which prevents from the direct use of Theorem

2.4 of[START_REF] Mertikopoulos | On the convergence of gradient-like ows with noisy gradient input[END_REF]. Indeed a such Lipschitz property depends on the distribution of the random variable Z. Nevertheless, the nature of the gradient dynamics prevents the explosion of the solutions, which compacties the trajectories and permits to by-pass the uncertainty on this Lipschitz property.

Acknowledgments

The authors gratefully thank Gersende Fort for her numerous insights and stimulating discussions about some earlier form of this work. This work is partially supported by the Fondation Simone et Cino Del Duca through the project OpSi-MorE, and by the French Agence Nationale de la Recherche (ANR), project under reference ANR-PRC-

CE23 MASDOL.

A Lagrangian formulation

This appendix is dedicated to the proof of the Lagrangian unconstrained formulation used to deal with the CV@R α constraint.

Proof of Proposition 2.1. We observe that J(u) = -m i=1 u i E[Z i ] is a continuous convex function, dened on the convex simplex ∆ m and P M is a convex optimization problem dened with a constraint g M dened through a CV@R α inequality. It is well known that CV@R α induces a coherent measure of risk, which implies the two important properties of positive homogeneity and sub-additivity. We refer to [START_REF] Ch | Some remarks on the value-at-risk and the conditional value-at-risk[END_REF] for this important remark, and to [START_REF] Artzner | Thinking coherently[END_REF][START_REF] Artzner | Coherent measures of risk[END_REF] for seminal contributions on risk measurements in mathematical nance. The important consequence of a such coherence is the convexity of the function that denes the constraint in P M . We verify that ∀λ ∈ (0, 1),

We are led to dene the convex function g M that induces the constraint in P M by:

We shall observe that the collection of functions (g M ) M ∈R denes a family of convex constraints and ∃M 0 ∈ R such that for all M > M 0 , the constraint g M is strictly feasible: there exists u ∈ ∆ m that is an interior point such that g M (u) < 0. For a such M > M 0 , we then dene the Lagrangian function:

From the Slater condition, we deduce that the strong duality involved in P M holds: there is no duality gap and primal and dual optimization problems coincide. In particular, the Karush-Kuhn-Tucker conditions are satised: the solution u ⋆ M of the primal problem P M veries:

and the complementary slackness is veried for the pair (u ⋆ M , λ ⋆ M ):

We now establish the monotonous property. We consider a pair (M 1 , M 2 ) and consider (u ⋆ 1 , λ ⋆ 1 ) and (u ⋆ 2 , λ ⋆ 2 ) derived from the previous Lagrangian equation. We know that the complementary slackness condition holds for both (u ⋆ 1 , λ ⋆ 1 ) and (u ⋆ 2 , λ ⋆ 2 ). Equations ( 16) and ( 17) entail:

where we used in the last line Equation ( 17)

1 ≥ 0 and λ ⋆ 2 ≥ 0, we deduce from Equation ( 18) that:

It nally implies that

We now consider the penalized criterion with λ 1 > λ 2 , we verify that:

It implies that:

Obviously, we also observe that v λ solves P M for M = CV@R α (v λ ), e.g.:

B Proof of the theoretical results on the SMD B.1 Almost sure convergence of the algorithm Below, we will use some standard results that are valid for any Bregman divergence D Φ , whose statements are given below. We refer to [START_REF] Nemirovskij | Problem complexity and method eciency in optimization[END_REF] for further details.