<figure><img src="image_1.png" /></figure>

HAL
open science

CV@R penalized portfolio optimization with biased
stochastic mirror descent

Manon Costa, Sébastien Gadat, Lorick Huang

► To cite this version:

Manon Costa, Sébastien Gadat, Lorick Huang. CV@R penalized portfolio optimization with biased stochastic mirror descent. Finance and Stochastics, 2025, 29 (3), pp.609-664. 10.1007/s00780-025-00568-2 . hal-05147991

HAL Id: hal-05147991

https://hal.science/hal-05147991v1

Submitted on 7 Jul 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
“CV@R penalized portfolio optimization with biased
stochastic mirror descent”

Sébastien Gadat, Manon Costa and Lorick Huang
# CV@R penalized portfolio optimization with biased stochastic mirror descent

M. Costa,\* S. Gadat,† L. Huang‡

July 28, 2023

## Abstract

This article studies and solves the problem of optimal portfolio allocation with CV@R penalty when dealing with imperfectly simulated financial assets. We use a Stochastic biased Mirror Descent to find optimal resource allocation for a portfolio whose underlying assets cannot be generated exactly and may only be approximated with a numerical scheme that satisfies suitable error bounds, under a risk management constraint. We establish almost sure asymptotic properties as well as the rate of convergence for the averaged algorithm. We then focus on the optimal tuning of the overall procedure to obtain an optimized numerical cost. Our results are then illustrated numerically on simulated as well as real data sets.

**Keywords:** Stochastic Mirror Descent, biased observations, risk management constraint, portfolio selection, discretization.

**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 OpSiMorE, and by the French Agence Nationale de la Recherche (ANR), project under reference ANR-PRC-CE23 MASDOL.

## 1 Introduction

Given a portfolio of financial assets, what is the best way to allocate resources under a risk management constraint? Without any fixed 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 profits 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 financial crises. Risks can be quantified with various metrics, and policies developed using one or another metric could lead to different 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

<footnote>\*Institut de Mathématiques de Toulouse. CNRS UMR 5219. Université Paul Sabatier, 118 route de Narbonne, F-31062 Toulouse cedex 09. Email: manon.costa@math.univ-toulouse.fr</footnote>

<footnote>†Toulouse School of Economics, CNRS UMR 5314, Université Toulouse 1 Capitole, Esplanade de l'Université, Toulouse, France. Institut Universitaire de France. Email: sebastien.gadat@tse-fr.eu</footnote>

<footnote>‡INSA de Toulouse, IMT UMR CNRS 5219, Université de Toulouse, 135 avenue de Rangueil 31077 Toulouse Cedex 4 France. Email: lorick.huang@insa-toulouse.fr</footnote>
literature) of the portfolio, which stands for a coherent measure of risk of an investment [6, 44]. We refer to Equation (2) in the next section for a formal definition. 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 fixed threshold.

Portfolio optimization under CV@R constraint has already been investigated in previous works, we can notably cite Rockafellar and Uryasev [43] who develop a Monte-Carlo strategy to approximate the CV@R of the portfolio. These results have been refined in subsequent works of the authors, notably in [29], [44] and more recently in [42], 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 [13] and [38], 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 [10], [11], and the book [9], as well as [23] all dealing with Robbins-Monro stochastic approximation. However, as in [23], 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 [7] or [31] (see also the references therein) for an overview of the techniques used for biased stochastic approximation algorithms.

In mathematical finance, assets are commonly modeled using stochastic differential 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 effect of this bias should be addressed to finely tune the optimization algorithm. Discretization for SDEs has been studied since the pioneering work of Talay and Tubaro [46], where strong regularity is assumed on the coefficients. As results developed, regularity assumptions have been dropped, and sharp bounds have been derived for both weak and strong error (see e.g. [28]) on specific stochastic models for specific discretization schemes (see e.g. [1, 2] 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) [19] process and we use this rate for the growth ratio of one component of our portfolio (see Subsection 3.1 below the definition). 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 [1], that is known to preserve positivity (see Section 3.2 and Proposition 12 below for more details). Thankfully, rate of convergence for the approximation of the CIR process has been widely studied, see e.g. [21, 3] and these quantitative results may be exploited to control the bias and its propagation in the specific 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, \dots, v_n)$ such that $v_1 + \cdots + v_n = V$. Equivalently, we consider $u = (u_1, \dots, u_n)$ in the simplex: $u_i \in [0, 1]$ and $u_1 + \cdots + 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 financial institutions.

In order to determine a strategy $u = (u_1, \dots, u_n)$ in the simplex, one natural idea would be to derive an algorithm without the constraint $u_1 + \cdots + 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 [35]. 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 [36] that initiates the use of a noisy stochastic gradient in the mirror descent, [30] for specific derivations on convex optimization problems on the simplex and more recently Zhou et al [47] 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. [15]).

In summary, to find an investment strategy while managing risk, we will use a stochastic mirror descent with an approximation of the portfolio, combined with a fine 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 (O.D.E. 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 influence 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 fine tuning of the numerical implicit Euler scheme.

In what follows, we carefully study the influence 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 [33] that crosses the risk (built with the CV@R) and the expected return to identify the optimal investment strategy. With this approach, our main results are

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

## 2.1 Financial model and constrained optimization

We consider the relative return of a portfolio, described by a random vector $Z = (Z_1, \dots, 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 fixed time horizon $t_0 + T$ where $T = 1$, i.e. each random variable $Z_i$ corresponds to $Z_i = \frac{A_i(T)}{A_i(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, \dots, u_m)$ that belongs to the $m-1$ dimensional simplex denoted by $\Delta_m$ and defined as:

$$ \Delta_m := \left\{ u \in \mathbb{R}_+^m : \sum_{i=1}^m u_i = 1 \right\}. $$
We are interested in a constrained optimization of the mean of the random variable $\langle Z, u \rangle$ defined by:

$$ \langle Z, u \rangle = \sum_{i=1}^{m} u_i Z_i. \qquad (1) $$

The constraint we are using involves a risk measure, named the conditional value at risk, which quantifies the mean value of a loss given that a loss has occurred. More precisely, if $\alpha$ refers to a chosen risk level, $V@R_{\alpha}(u)$ corresponds to the classical statistical $\alpha$-quantile defined as:

$$ V@R_{\alpha}(u) := \sup\{q \in \mathbb{R} : \mathbb{P}(\langle Z, u \rangle \le q) \le \alpha\}, $$

while $CV@R_{\alpha}$, that defines the active constraint we will use, is the mean value of the loss when $\langle Z, u \rangle$ is below $V@R_{\alpha}$, namely:

$$ CV@R_{\alpha}(u) = \mathbb{E}[\langle Z, u \rangle | \langle Z, u \rangle \le V@R_{\alpha}(u)]. \qquad (2) $$

Note that here, we implicitly assume that the risk level $\alpha$ is chosen such that $V@R_{\alpha}(u)$ is negative, which implies that $CV@R_{\alpha}$ is a positive quantity.

Then, for any fixed positive level $M$, the optimization problem of $u \mapsto \mathbb{E}[\langle Z, u \rangle]$ with $CV@R_{\alpha}$ constraints we are interested in, is the problem defined by:

$$ \begin{aligned} \mathcal{P}_M &:= \arg\max_{u \in \Delta_m} \left\{ \sum_{i=1}^{m} u_i \mathbb{E}[Z_i] : \text{CV@R}_{\alpha}(u) \le M \right\} \\ &= \arg\min_{u \in \Delta_m} \left\{ -\sum_{i=1}^{m} u_i \mathbb{E}[Z_i] : \text{CV@R}_{\alpha}(u) \le M \right\}. \end{aligned} \qquad (3) $$

For our purpose, we also introduce the unconstrained penalized optimization problem: we consider for any $\lambda > 0$ the solution $v_{\lambda}$ of the following (convex) optimization problem:

$$ Q_{\lambda} := \arg \min_{u \in \Delta_m} \left\{ -\sum_{i=1}^{m} u_i \mathbb{E}[Z_i] + \lambda \mathbb{C}V@R_{\alpha}(u) \right\}. $$

The (standard) following result is key for our forthcoming analysis.

**Proposition 1.** The collection of convex problems $(\mathcal{P}_M)_{M>0}$ and $(Q_{\lambda})_{\lambda>0}$ are equivalent. More precisely, for any $M > 0$ that defines a feasible constraint, a solution $u_M^*$ exists such that:

$$ \exists \lambda_M^* > 0 \quad u_M^* = \arg \min_{u \in \Delta_m} \left\{ -\sum_{i=1}^{m} u_i \mathbb{E}[Z_i] + \lambda_M^* \mathbb{C}V@R_{\alpha}(u) \right\}. $$

Moreover, $\lambda_M^*$ is a decreasing function of $M$. Oppositely, any solution $v_{\lambda}$ of $Q_{\lambda}$ solves $\mathcal{P}_M$ with $M = CV@R_{\alpha}(v_{\lambda})$.

We emphasize that this latter result is the keystone of our forthcoming analysis: instead of solving directly $\mathcal{P}_M$, we will instead define a collection of penalized problems $Q_{\lambda}$ for several values of $\lambda$. Proposition 1 then establishes that this approach leads to optimal portfolio allocation under $CV@R$ constraints. We also highlight that the relationship between $M$ and $\lambda_M^*$ remains unknown and challenging to obtain. The proof of this Lagrangian formulation is deferred to Appendix A. Using the result of [43, 29] and in particular the convex representation of the $CV@R_{\alpha}$, we observe that:

$$ \mathbb{C}V@R_{\alpha}(u) = \min_{\theta \in \mathbb{R}} \psi_{\alpha}(u, \theta), $$
where $\psi$ is the convex coercive Lipschitz continuous and differentiable function defined by:

$$ \psi_{\alpha}(u, \theta) = \theta + \frac{1}{1-\alpha} \mathbb{E}[\langle Z, u \rangle - \theta]_{+}, \quad (4) $$

where $\lfloor x \rfloor_{+} = \max(0, x)$. We emphasize that despite the definition with $\lfloor \cdot \rfloor_{+}$, when $Z$ has an absolutely continuous distribution with respect to the Lebesgue measure, $\psi_{\alpha}$ is differentiable (we refer to [10] for further details on this remark).

Then, we deduce from Proposition 1 and the above representation that the collection of the optimization problems $(\mathcal{P}_M)_{M>0}$ are then described equivalently by the convex unconstrained problem:

$$ Q_{\lambda} = \arg \min_{(u, \theta) \in \Delta_{m} \times \mathbb{R}} \{p_{\lambda}(u, \theta)\}, \quad (5) $$

where the key function $p_{\lambda}$ is defined by:

$$ p_{\lambda}(u, \theta) = - \sum_{i=1}^{m} u_i \mathbb{E}[Z_i] + \lambda \psi_{\alpha}(\theta, u). \quad (6) $$

We emphasize that $(p_{\lambda})_{\lambda>0}$ forms a collection of convex functions defined on $\Delta_m \times \mathbb{R}$, coercive w.r.t. $\theta$ and that $Q_{\lambda}$ enables to avoid the CV@R$_{\alpha}$ constraint with a suitable reparametrization of the problem.

## 2.2 Biased stochastic Mirror Descent

To solve the minimization problem (5) we are led to use stochastic approximation theory, that originates in the pioneering works of [40] and [41] since the (convex) function $p_{\lambda}$ is written as an expectation.

However, we encounter here specific difficulties. First we need to handle the minimization over the simplex $\Delta_m$. Second, the random variables $Z$ involved in $p_{\lambda}$ 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_{\lambda}$ is differentiable, it is the expectation of a non-differentiable function of $(u, \theta)$ which will require some specific attention in the following algorithms.

### 2.2.1 Deterministic case

Mirror Descent (MD below) originates from the pioneering work of [35] 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 $\Delta_m$ (see e.g. [30],[16]). The MD approach has the nice feature to define a smooth evolution that lives inside the constrained set without adding some supplementary projection step and “pushes” the frontiers of the simplex at an infinite distance from any point strictly inside $\Delta_m$. Though we will need to encompass the stochastic framework, since $p_{\lambda}$ is defined through an expectation which is not explicit, for the sake of clarity we describe first the deterministic version.

For this purpose, we first introduce the strongly convex negative entropy function over the simplex $\Delta_m$ of probability distributions, defined by:

$$ \forall u \in \Delta_m, \quad \varphi(u) = \sum_{i=1}^{m} u_i \log u_i. $$

If $\langle \cdot, \cdot \rangle$ refers to the standard Euclidean inner product, the Bregman divergence $D_{\varphi}$ is defined by:

$$ \forall (u, v) \in \Delta_m^2, \quad D_{\varphi}(u, v) = \varphi(u) - \varphi(v) - \langle \nabla \varphi(v), u - v \rangle. $$

This Bregman divergence will induce the natural metric associated to the component $u$ involved in the problem $Q_{\lambda}$. In the same way, we also introduce the standard Euclidean square distance over $\mathbb{R}^2$,
which may be viewed as a Bregman divergence. It then leads to a Bregman divergence over pairs $(\theta, u)$ and $(\theta', v)$ defined as:

$$ \forall (\theta, \theta') \in \mathbb{R}^2, \quad \forall (u, v) \in \Delta_m^2, \quad \mathcal{D}_{\Phi}((u, \theta), (v, \theta')) = \frac{(\theta - \theta')^2}{2} + \mathcal{D}_{\varphi}(u, v), $$

which is the Bregman divergence associated to the strongly convex function:

$$ \Phi(u, \theta) = \varphi(u) + \frac{\theta^2}{2}. $$

We emphasize that the strong convexity of $\Phi$ induces the following lower bound:

$$ \forall (\theta, \theta') \in \mathbb{R}^2, \quad \forall (u, v) \in \Delta_m^2, \quad \mathcal{D}_{\Phi}((\theta, u), (\theta', v)) \geq (\theta - \theta')^2 + \frac{1}{2} \|u - v\|^2. \quad (7) $$

The (deterministic) MD with a step-size sequence $(\eta_k)_{k \ge 1}$ consists in minimizing from $k$ to $k + 1$ the first order Taylor approximation of $p_\lambda$ penalized by the Bregman divergence. To alleviate the notations, we will denote by $X_k = (U_k, \theta_k)$ the position of the algorithm and generally $x = (u, \theta)$ an element of $\Delta_m \times \mathbb{R}$. We observe that the MD iterative step corresponds to:

$$ X_{k+1} = \arg \min_{x \in \Delta_m \times \mathbb{R}} \left\{ \langle \nabla p_{\lambda}(X_k), x - X_k \rangle + \frac{1}{\eta_{k+1}} \mathcal{D}_{\Phi}(x, X_k) \right\}. $$

A remarkable feature is that this latter minimization can be made explicit :

$$ X_{k+1} = \begin{pmatrix} U^{k+1} \\ \theta^{k+1} \end{pmatrix} \quad \text{with} \quad \begin{cases} U^{k+1} = \frac{U^k e^{-\eta_{k+1} \partial_u p_\lambda(U^k, \theta^k)}}{\|U^k e^{-\eta_{k+1} \partial_u p_\lambda(U^k, \theta^k)}\|_1}, \\ \theta^{k+1} = \theta^k - \eta_{k+1} \partial_\theta p_\lambda(U^k, \theta^k) \end{cases}, \quad (8) $$

where the first equation has to be understood within a $m$ dimensional vector structure.

For a fixed $\lambda > 0$, we define $x_\lambda^*$ the minimizer of $p_\lambda$. Following the work of [12, 35], it can be shown that an averaged version of the sequence $(X_k)_{k \ge 1}$ defined by

$$ \tilde{\chi}_n = \left( \sum_{k=0}^{n} \eta_k \right)^{-1} \sum_{k=0}^{n} \eta_k X_k, $$

satisfies the next error bound:

**Theorem 2** (Proposition 1 in [30]). For any $\lambda > 0$ and $L = \arg \max_{x \in \Delta_m \times \mathbb{R}} \{\|\nabla p_\lambda(x)\|\}$, then:

$$ p_{\lambda}(\tilde{\chi}_n) - p_{\lambda}(x_{\lambda}^*) \le \frac{\{\Delta_{\Phi}^0\}^2 + L^2 \sum_{k=1}^{n} \eta_k^2}{2 \sum_{k=1}^{n} \eta_k}, $$

$$ \text{where } \{\Delta_{\Phi}^0\}^2 = \frac{(\theta_0 - V @ R_{\alpha}(u_{\lambda}^*))^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 (9) and (10) stated below, which implies:

$$ L^2 \le \frac{m(1 - \alpha + \lambda)^2 \mathbb{E}[\|Z\|^2] + \lambda^2 \alpha^2}{(1 - \alpha)^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 finely tune the step-size sequence to obtain a $O(n^{-1/2})$ upper bound. We refer to [30] for further details.
### 2.2.2 Stochastic Mirror Descent (SMD) using biased simulations

In our problem, $p_λ$ and $\nabla 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 $\partial_u$ and $\partial_\theta$ the partial derivatives with respect to $u$ and $\theta$, and using that $p_λ$ is expectation of a convex function, we then verify that:

$$ \partial_u p_\lambda(u, \theta) = -\mathbb{E}[Z] + \frac{\lambda \mathbb{E}[Z \mathbf{1}_{\langle Z, u \rangle \ge \theta}]}{1 - \alpha}, \quad (9) $$

and

$$ \partial_{\theta} p_{\lambda}(u, \theta) = \lambda \left[ 1 - \frac{1}{1-\alpha} \mathbb{E}[\mathbf{1}_{\langle Z, u \rangle \ge \theta}] \right]. \quad (10) $$

**Remark 3.** We formulated the above convergence and rate estimation results for $\nabla 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 $(\hat{Z}^k)_{k \ge 0}$ that are also sampled independently from the previous positions of the algorithm. The expressions (9) and (10) lead to a natural (possibly biased) stochastic approximation of $\nabla p_λ$ with the help of the sequence $(\hat{Z}^k)_{k \ge 0}$. Assuming that the algorithm is at step $k$ at position $(U_k, \theta_k)$, we introduce the stochastic approximation of the sub-gradients:

$$ \left\{ \begin{aligned} \hat{g}_{k+1,1} &= -\hat{Z}^{k+1} + \frac{\lambda}{1-\alpha} \hat{Z}^{k+1} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k} \\ \hat{g}_{k+1,2} &= \lambda \left[ 1 - \frac{1}{1-\alpha} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k} \right] \end{aligned} \right. . \quad (11) $$

We now describe the necessary assumptions on the sequence $(\hat{Z}^k)_{k \ge 0}$ to build a consistent SMD algorithm.

#### Assumption (H₁) Assumptions on the biased simulations

Let $\mathcal{L}(\hat{Z}^{k+1})$ and $\mathcal{L}(Z)$ denote the distributions of $\hat{Z}^{k+1}$ and $Z$ respectively. We assume that the sequence $(\hat{Z}^k)_{k \ge 0}$ satisfies both:

$$ \mathcal{W}_1(\mathcal{L}(\hat{Z}^{k+1}), \mathcal{L}(Z)) \le \delta_{k+1}, \quad (12) $$

and

$$ \forall u \in \mathbf{\Delta}_m, \quad \forall \theta \in \mathbb{R}, \quad \left\| \mathbb{E} \left[ \langle Z, u \rangle \mathbf{1}_{\langle Z, u \rangle \ge \theta} - \langle \hat{Z}^{k+1}, u \rangle \mathbf{1}_{\langle \hat{Z}^{k+1}, u \rangle \ge \theta} \mid \mathcal{F}_k \right] \right\| \le v_{k+1}. \quad (13) $$

where $\mathcal{W}_1$ stands for the Wasserstein-1 distance.

The sequences $(\delta_{k+1})_{k \ge 0}$ and $(v_{k+1})_{k \ge 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 $\delta_{k+1} = v_{k+1} = 0$ for all integer $k$. The Wasserstein distance used in $\delta_{k+1}$ is indeed an easiest way to upper bound the bias naturally involved in our approximation procedure. Indeed, Equation (29) (see below) will involve the Kolmogorov distance instead, denoted by $d_{Kol}$ below, which quantifies the difference 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 straightforward inequality $d_{Kol} \le 2\sqrt{W_1}$, we prefer to introduce Assumption (H₁) in terms of $W_1$.

We emphasize that both $(\delta_{k+1})_{k \ge 0}$ and $(v_{k+1})_{k \ge 0}$ might heavily depend on $m$ the dimension of the vector $Z$. A specific example will be detailed below in our work.

The method we propose is then defined with Algorithm 1. The next result states the asymptotic almost sure convergence of the sequence $(X_k)_{k \ge 0}$ constructed in Algorithm 1 towards the target point $x_\lambda^* = (u_\lambda^*, V@R_\alpha(u_\lambda^*))$. Under a suitable assumption on the biased sequences $(\delta_{k+1})_{k \ge 0}$ and $(v_{k+1})_{k \ge 0}$, we derive the next result.
**Data:** Step-size sequence {$\eta_n$, $n \in \mathbb{N}$} and $U_0 \in \mathbb{R}$, $\theta_0 \in \mathbb{R}$; $\alpha \in (0, 1)$

**Result:** Two sequences: $X_k = (U_k, \theta_k)_{k \ge 0}$

1. for $k = 0, \dots, \mathbf{do}$

2. Simulate the random $Z^{k+1}$ satisfying (12) and (13);

3. Compute a stochastic approximation $\hat{g}_{k+1}$ of $\nabla p_{\lambda}(U_k, \theta_k)$ with:

$$ 
\begin{cases}
\hat{g}_{k+1,1} = -\hat{Z}^{k+1} + \frac{\lambda}{1-\alpha} \hat{Z}^{k+1} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k} \\
\hat{g}_{k+1,2} = \lambda \left[ 1 - \frac{1}{1-\alpha} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k} \right].
\end{cases}
$$

Update the algorithm $X_{k+1} = \arg \min_{x \in \Delta_m \times \mathbb{R}} \{\langle \hat{g}_{k+1}, x - X_k \rangle + \frac{1}{\eta_{k+1}} \mathcal{D}_{\Phi}(x, X_k)\}$ using:

$$ X_{k+1} = (U_{k+1}, \theta_{k+1}), \quad \begin{cases} U^{k+1} = \frac{U^k e^{-\eta_{k+1} \hat{g}_{k+1,1}}}{\|U^k e^{-\eta_{k+1} \hat{g}_{k+1,1}}\|_1}, \\ \theta^{k+1} = \theta^k - \eta_{k+1} \hat{g}_{k+1,2} \end{cases}. $$

<caption>**Algorithm 1:** Biased SMD</caption>

**Theorem 4** (Almost sure convergence of the biased SMD). Assume that

$$ \sum_{k \ge 0} \eta_{k+1} = +\infty, \quad \text{and} \quad \sum_{k \ge 0} \eta_{k+1}^2 < +\infty, $$

and that the bias sequences $(\delta_{k+1})_{k \ge 0}$ and $(v_{k+1})_{k \ge 0}$ introduced in (H$_1$) satisfy

$$ \sum_{k \ge 0} \eta_{k+1} (\sqrt{\delta_{k+1}} + v_{k+1}) < +\infty, $$

then:

$$ i) \quad \sum_{k \ge 0} \eta_{k+1} [p_{\lambda}(X_k) - \min(p_{\lambda})] \le \sum_{k \ge 0} \eta_{k+1} \langle \nabla p_{\lambda}(X_k), X_k - x_{\lambda}^* \rangle < +\infty \quad \text{a.s.} $$

ii) The Cesaro average $\bar{X}_k^\eta$ defined by

$$ \bar{X}_k^\eta := \left( \sum_{i=0}^k \eta_i \right)^{-1} \left( \sum_{i=0}^k \eta_i X_i \right) \quad (14) $$

is almost surely convergent and

$$ p_{\lambda}(\bar{X}_k^\eta) \to \min(p_{\lambda}) \quad \text{a.s.} $$

iii) Assume that:

$$ \sum_{k \ge 0} \sqrt{\delta_{k+1}} + v_{k+1} < +\infty, $$

then the sequence $(X_k)_{k \ge 1}$ almost surely converges and verifies that:

$$ \lim_{k \to +\infty} p_{\lambda}(X_k) = \min(p_{\lambda}) \quad \text{a.s.} $$
The results obtained in i) and ii) are rather classical in the stochastic optimization community. The difficulty here is to derive a “Lyapunov” function and then to quantify the effect of the bias induced by $\hat{Z}^{k+1}$ satisfying (12) and (13). To do so, we will follow the arguments developed in [30] 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 Sigmund theorem. The last item is much more difficult to obtain, roughly speaking we use the O.D.E. method introduced in [15]. 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 [15]. In particular, we will need to use and adapt some fine results of [34] to our framework. The proof is given in Appendix B.

Asymptotic results such as the one stated in Theorem 4 are commonly criticized because algorithms are always ran within a finite number of iterations. In particular, Theorem 4 does not provide any insight on the rate of convergence of the method, which may crucially influence the practical usefulness of the algorithm. Nevertheless, it shows that Algorithm 1 converges and minimizes the function $p_\lambda$ up to a Cesaro averaging procedure (see i) and ii) of Theorem 4). 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 quantified by $\delta_{k+1}$ and $v_{k+1}$.

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

**Theorem 5** (Finite-time guarantees). Let us consider $(X_k)_{k \ge 0}$ defined in Algorithm 1 and its Cesaro averaging counterpart introduced in (14), then for any $n > 1$,

$$ \mathbb{E}[p_{\lambda}(\bar{X}_{n}^{\eta})] - p_{\lambda}(x_{\lambda}^{*}) \le \left( \sum_{j=0}^{n-1} \eta_{j+1} \right)^{-1} \left( \mathbf{D}_{\Phi}^{0} + \sum_{k=0}^{n-1} [a_{k+1} \mathbf{D}_{\Phi}^{k} + b_{k+1}] \right) $$

where

$$ a_{k+1} = 2\eta_{k+1} \left( 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda v_{k+1}}{1-\alpha} \right), \quad b_{k+1} = C(\eta_{k+1}^2 + a_{k+1}) \quad \text{and} \quad \mathbf{D}_{\Phi}^{k} = \mathbb{E}\mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k}). $$

Finally, we aim at applying these guarantees in the case of a finite horizon of computation $n$ and in this view, $\mathbf{D}_{\Phi}^{k}$ can be controlled using a recursion inequality (41). We consider the specific case of the SMD with a constant step-size sequence stopped at iteration $n$:

$$ \eta_{k+1} = \eta > 0, \quad \forall 0 \le k \le n. $$

We will also assume a constant upper bound of the bias in the simulation of the random variables $\hat{Z}^k$: we denote by $\omega$ 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 fixed values of $\delta_{k+1}$ and $v_{k+1}$ such that:

$$ 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda v_{k+1}}{1-\alpha} = \omega > 0, \quad \forall 1 \le k \le n. $$

**Corollary 6.** For a given $n \in \mathbb{N}$, if $(\eta, \omega)$ are chosen such that $\eta = \frac{\Delta_{\Phi}^{0}}{2\sqrt{n+1}}$ and $\omega = \frac{1}{\sqrt{n+1}\Delta_{\Phi}^{0}}$ with $\{\Delta_{\Phi}^{0}\}^{2} = \frac{(\theta_{0} - V@R_{\alpha}(u_{\lambda}^{*}))^{2}}{2} + \log m$, then there exists $C > 0$ large enough such that:

$$ \mathbb{E}[p_{\lambda}(\hat{X}_{n}^{\eta})] - p_{\lambda}(x_{\lambda}^{*}) \le C \frac{|\theta_{0} - V@R_{\alpha}(u_{\lambda}^{*})| + \sqrt{\log m}}{\sqrt{n+1}}. $$

This previous corollary is built using the optimal tuning of the parameters $\eta$ and $\omega$ 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 $\Delta_{\Phi}^0$. Nevertheless, any upper bound of $\Delta_{\Phi}^0$ may be used to derive a strategy and an upper bound of the excess risk (see Equation (43) in Appendix (B.2)). In the meantime, the choice:

$$ \eta = \frac{1}{2\sqrt{n+1}} \quad \text{and} \quad \omega = \frac{1}{\sqrt{n+1}}, $$

yields:

$$ \mathbb{E} [p_{\lambda}(\hat{X}_{n}^{\eta})] - p_{\lambda}(x_{\lambda}^{*}) \leq C \frac{|\theta_{0} - \text{V@R}_{\alpha}(u_{\lambda}^{*})|^{2} + \log m}{\sqrt{n+1}}, $$

which is slightly worse than the upper bound stated in Corollary (6).

# 3 Portfolio hedging under CV@R$_{α}$ constraint

We aim to apply our optimization strategy to the specific situation of portfolio hedging. As a consequence, we precise the structure of Z, a portfolio of m financial 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 field of application described in Section 2.

## 3.1 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 \ge 0}$ with no drift, and a family $S = (S^1, \dots, 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 [19] as a diffusion 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 \mapsto r_t$, whose stochastic differential equation depends on a triple $(a, b, \sigma_0)$ and is given by:

$$ dr_t = a(b - r_t)dt + \sigma_0 \sqrt{r_t} dB_0(t), \quad (15) $$

where $(B_0(t))_{t \ge 0}$ stands for a standard real Brownian motion. The parameter b stands for the long-time mean of the short rate while a quantifies the strength of the mean-reversion effect. The volatility $\sigma_0$ is multiplied by $\sqrt{r_t}$, and the condition $2ab > \sigma_0^2$ guarantees that the interest rate remains positive with probability 1. We refer to [24] 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_t^1, \dots, S_t^{m'})$ are then described by the following system of stochastic differential equations:

$$ \forall t \ge 0 \qquad \begin{cases} dY_t = r_t Y_t dt, \\ dS_t^i = \mu_i S_t^i dt + \sigma_i S_t^i dB_i(t), \quad \forall i \in \{1, \dots, m'\}, \end{cases} \quad (16) $$

where $\mathbf{B} = (B_0, B_1, \dots, B_{m'})$ refers to a multivariate Brownian motion with correlated components. For the sake of simplicity, these components are assumed to satisfy:

$$ \mathbb{E}[B_i(t)B_j(t)] = \rho_{i,j}t, $$

but more general correlation structures could be handled in our framework with further efforts. The correlation matrix is the symmetric positive matrix, denoted by $\Sigma = (\rho_{i,j})_{1 \le (i,j) \le m}$.

The first process $(Y_t)_{t \le 1}$ corresponds to the neutral risk process whereas each geometric Brownian motion $(S_t^i)_{t \le 1}$ corresponds to some specific risky assets parametrized by known parameters $\mu_i$ and $\sigma_i$. We refer to Pitman and Yor [39] and to Gulisashvili and Stein [26] 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 (15) and (16)) are known and we are simply interested in the optimal hedging strategy, i.e. we are looking for the optimal $u_{\lambda}^*$ associated to the optimal mean return penalized by $\lambda CV@R_{\alpha}$ defined in Equation (5) 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 (H$_2$) Assumptions on the portfolio parameters

We assume that:

i) the CIR parameters satisfy $ab > \sigma_0^2$ and $a > 2\sqrt{2}\sigma_0$.

ii) the correlation matrix of the Brownian motions $\Sigma$ in invertible.

These assumptions on the coefficients defining 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).

## 3.2 Simulation of the portfolio

We are interested in an efficient simulation method that satisfies assumptions (H$_1$): both a good approximation of the law of $Z_1$ (12) and a bias upper bound of the form (13).

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

$$ \forall t \ge 0 \qquad S_t^i = S_0^i \exp(\mu_i t + \sigma_i B_i(t)). \qquad (17) $$

In the meantime, we can also observe that some exact simulation method exists for the CIR process sampled at a given fixed time. Nevertheless, the CIR process induces some numerical difficulties.

* First, the SDE (15) is not explicitly solvable at any time $t \in [0, 1]$ so that the integral of the CIR between 0 and 1 needs to be approximated.

* Second the correlations between the different 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 financial modeling, is hardly tractable with the use of an exact simulation of the CIR (with the help of $\chi^2$ square distributions, see, e.g. [4]).

* 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 fixed 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, \xi_{k+1}), $$

where $F$ corresponds to an iterative update that uses the position $(r_{kh}, Y_{kh})$ and $\xi_{k+1}$ is a two dimensional Gaussian innovation with a specific covariance. Nevertheless, such an attractive point of view induces some theoretical difficulties because of the form of the function $Z1_{(Z,u)\ge\theta}$ we need to approximate. This function is non-smooth and unbounded, which generates some technical difficulties 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), \dots, 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), \dots, B_{m'}(t))$. More precisely, the Cholesky decomposition used on $\Sigma$ yields a matrix $L$ such that:

$$LL^T = \Sigma. \qquad (18)$$

We now define $B(t) = LW(t)$. There are two important consequences of constructing $(B_i(t))_{i \ge 0}$ for $i = 0, \dots, m'$ in this way. The first one is that the Cholesky method yields a lower triangular matrix $L$ whose first line is $1, 0, \dots, 0$:

$$L = \begin{pmatrix} 1 & 0 & \dots & 0 \\ \ell_{01} & \ell_{11} & \dots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ \ell_{0m'} & \dots & \dots & \ell_{m'm'} \end{pmatrix} = \begin{pmatrix} 1 & 0 & \dots & 0 \\ \ell_{01} & \ell_{11} & \dots & \tilde{L} \\ \vdots & \vdots & \ddots & \vdots \\ \ell_{0m'} & \dots & \dots & \ell_{m'm'} \end{pmatrix}.$$

In particular, we see that we can take $B_0(t) = W_0(t)$ and $\tilde{W}(t) = (W_1(t), \dots, W_{m'}(t))$ is independent of $W_0(t)$. The $m-1$ components may be simply written as:

$$\left\{ \begin{aligned} B_1(t) &= \ell_{01}W_0(t) + \ell_{11}W_1(t) \\ B_2(t) &= \ell_{02}W_0(t) + \ell_{12}W_1(t) + \ell_{22}W_2(t) \\ \vdots \\ B_{m'}(t) &= \ell_{0m'}W_0(t) + \ell_{1m'}W_1(t) + \dots + \ell_{m'm'}W_{m'}(t). \end{aligned} \right.$$

The second important and standard consequence is that these stochastic processes $(B_i(t))_{i \ge 0}$ for $i = 0, \dots, 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_1^i = S_i$, $B_i(1) = B_i$ and $W_i(1) = W_i$. The asset $i$ at time $1$ for all $i \in \{1, \dots, m'\}$ is a Geometric Brownian motion driven by $(B_i(t))_{i \ge 0}$, and we have:

$$ \begin{aligned} S_i &= S_0^i \exp \left( \left( \mu_i - \frac{(\sigma_i)^2}{2} \right) + \sigma_i B_i \right) \\ &= S_0^i \exp \left( \left( \mu_i - \frac{(\sigma_i)^2}{2} \right) + \sigma_i \left[ \ell_{0i} W_0 + \ell_{1i} W_1 + \dots + \ell_{ii} W_i \right] \right) \\ &= S_0^i \exp \left( \left( \mu_i - \frac{(\sigma_i)^2}{2} \right) + \sigma_i (\tilde{L} \tilde{W})_i \right) \times e^{\sigma_0 \ell_{0i} W_0}, \end{aligned} \qquad (19) $$

since $[\ell_{1i}W_1 + \dots + \ell_{ii}W_i] = (\tilde{L}\tilde{W})_i$ when $\tilde{W} = (W_1, \dots, 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 coefficients in the SDE are not uniformly elliptic and bounded, as assumed in the seminal works of Bally and Talay [8]. Besides, a classical explicit Euler scheme generates positivity issues (because of the square root). However, many authors, notably Alfonsi [1, 2] 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 [1] and studied in numerous articles [21, 3]. The drift-implicit Euler scheme
on a discrete time grid $(kh)_{0 \le k \le N}$ can be written by considering the SDE satisfied by $y_t = \sqrt{r_t}$ which leads to

$$ \hat{y}_{(k+1)h} = \hat{y}_{kh} + \left( \frac{4ab - \sigma_0^2}{8\hat{y}_{(k+1)h}} - \frac{a}{2}\hat{y}_{(k+1)h} \right) h + \frac{\sigma_0}{2}\Delta B_0^{(k)}, $$

where $\Delta B_0^{(k)} = B_0((k+1)h) - B_0(kh)$. This implicit scheme can be solved explicitly on $\mathbb{R}_+$ from iteration $k$ to iteration $k+1$, which ensures the positivity of the scheme $(\hat{r}_{kh})_{0 \le k \le N}$. In particular, the update from $kh$ to $(k+1)h$ is given by:

$$ \hat{r}_{(k+1)h} = \left( \frac{\sqrt{\hat{r}_{kh}} + \frac{\sigma_0}{2} \Delta B_0^{(k)}}{2\left(1 + \frac{ah}{2}\right)} + \sqrt{\frac{\left(\sqrt{\hat{r}_{kh}} + \frac{\sigma_0}{2} \Delta B_0^{(k)}\right)^2}{4\left(1 + \frac{ah}{2}\right)^2} + \frac{(4ab - \sigma_0^2)h}{8\left(1 + \frac{ah}{2}\right)}} \right)^2. \quad (20) $$

Convergence results for this discretization scheme will be recalled in Appendix C.

We then use the integral representation $Y_t = Y_0 \exp\left(\int_0^t r_s ds\right)$ and propose to approximate $Y_1$ with a Riemann integral approximation between 0 and 1:

$$ \hat{I}_h := \frac{1}{N} \sum_{k=1}^{N} \hat{r}_{kh}. $$

This leads to the definition of our approximation:

$$ \hat{Y}_1^{(h)} := Y_0 \exp(\hat{I}_h) = Y_0 \exp\left(\frac{1}{N} \sum_{k=1}^{N} \hat{r}_{kh}\right). \quad (21) $$

**Overall simulation algorithm** Combining all these steps, we are led to define the following algorithm to compute a biased simulation of $Z$.

**Data:** Parameters of the CIR: $a, b, \sigma_0$ and of the G.B.M. $(\mu_i, \sigma_i)$; Correlation matrix $\Sigma$

**Result:** A sample $\hat{Z}_1 = (\hat{Y}_1^{(h)}, S_1^1, \dots, S_1^{m'})$

1. Compute the Cholesky factorization of $\Sigma$ (see Equation (18)): $L$ and $\tilde{L}$.

2. Choose $\hat{r}_0$.

3. **for** $k = 1, \dots, N$ **do**

4. **if** Compute the recursive approximation $\hat{r}_{kh}$ using Equation (20).

5. Set $\hat{Y}_1^{(h)}$ as the Riemann approximation given in Equation (21);

6. Simulate the G.B.M. $S_1^1, \dots, S_1^{m'}$ with Equation (19).

**Algorithm 2:** Approximation of the portfolio $\hat{Z}_1$

In order to propose a convergent algorithm derived from Theorem 4, 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 (21). 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 (21). It is necessary to assess the accuracy of our method in terms of the value of $h$ to make (12) and (13) explicit. The simulation given in Algorithm 2 satisfies the next statements.

**Proposition 7.** Assume that the portfolio parameters satisfy Assumptions (H$_2$), then using the above discretization scheme presented in Algorithm 2, a constant $C$ exists (dependent on the CIR parameters) such that:

$$ W_1(\mathcal{L}(\hat{Z}_1), \mathcal{L}(Z_1)) = W_1(\mathcal{L}(\hat{Y}_1^{(h)}), \mathcal{L}(Y_1)) \le C\sqrt{h}. $$
**Proposition 8.** Assume that the portfolio parameters satisfy Assumptions (H₂), then for any $\epsilon > 0$, there exists a constant $K_\epsilon$ independent of $h$ and $m$ such that:

$$ \left\| \mathbb{E} \left[ \langle Z_1, w \rangle \mathbf{1}_{\langle Z_1, w \rangle \ge \theta} - \langle \hat{Z}_1, w \rangle \mathbf{1}_{\langle \hat{Z}_1, w \rangle \ge \theta} \right] \right\|_2 \le K_\epsilon \sqrt{me^{\frac{(\sigma^+)^2 m^2}{4\epsilon^2}}} h^{\frac{1}{6}-\epsilon}. $$

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 finance. They are not so easy to obtain in particular regarding the influence 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 efficient bias reduction strategy.

### 3.3 Stochastic Mirror Descent with sampling approximation

We finally 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 4, 5 and Propositions 7 and 8.

**SMD at fixed time horizon** $n$ According to our previous results, we can now consider the choice of parameters induced by Corollary 6. Recall that we assume a constant step-size sequence $\eta$ and discretizaton step-size $h$. Corollary 6 combined with Propositions 7 and 8 induce that $h$ should be chosen as:

$$ h^{1/4} + h^{\frac{1}{6}-\epsilon} \sim 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 4). Recall that the sequence $(\eta_k)_{k \ge 1}$ scales the step-size of stochastic gradient descent and assume that

$$ \eta_k = k^{-\alpha}, \quad \alpha \in (\frac{1}{2}, 1]. $$

We now chose the step sequence for approximating the CIR process as:

$$ h_k = h_0^{(m)} k^{-\beta}, \quad \beta > 0. $$

According to Propositions 7 and 8, using the notations of Theorem 4, we deduce that:

$$ \delta_k = k^{-\frac{\beta}{2}} \quad \text{and} \quad v_k = K_\epsilon \sqrt{me^{\frac{(\sigma^+)^2 m^2}{4\epsilon^2}}} h^{-\beta(\frac{1}{6}-\epsilon)}. $$

Assumption in Theorem 4 $\sum_{k \ge 0} \eta_{k+1}(\sqrt{\delta_{k+1}} + v_{k+1}) < +\infty$ now reads:

$$ \sum_{k \ge 1} k^{-\alpha} \left( k^{-\frac{\beta}{4}} + k^{-\beta(\frac{1}{6}-\epsilon)} \right) < \infty, $$

which is equivalent to:

$$ \alpha + \frac{\beta}{6} > 1, \quad \text{with} \quad \alpha \in (\frac{1}{2}, 1]. $$

With these conditions on $\alpha$ and $\beta$, we can now derive the non-asymptotic bound Theorem 5. We set $\epsilon = \alpha + \frac{\beta}{6} - 1 > 0$, using the notations of Theorem 5, we deduce that:

$$ a_{k+1} \sim k^{-1-\epsilon} \quad \text{and} \quad b_{k+1} \sim k^{-(2\alpha\wedge(1+\epsilon))}. $$
Using comparisons between series and integral, we find that:

$$ \mathbf{D}_{\phi}^{k} \leq D_{\phi}^{0} \exp(k^{-\epsilon}) + k^{-(2\alpha - \epsilon + 1)}. $$

If we use the notation $\lesssim$ 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:

$$ \begin{aligned} \mathbb{E}[p_{\lambda}(\bar{X}_{n}^{\eta})] - p_{\lambda}(x_{\lambda}^{*}) &\leq \frac{1}{\sum_{j=0}^{n-1} \eta_{j+1}} \left( \mathbf{D}_{\phi}^{0} + \sum_{k=0}^{n-1} [a_{k+1} \mathbf{D}_{\phi}^{k} + b_{k+1}] \right) \\ &\lesssim \frac{1}{n^{1-\alpha}} \left( \mathbf{D}_{\phi}^{0} + \sum_{k=0}^{n-1} \left[ k^{-1-\epsilon} [\mathbf{D}_{\phi}^{0} \exp(k^{-\epsilon})] + k^{-(2\alpha+\epsilon-1)} \right] + k^{-(2\alpha \wedge (1+\epsilon))} \right) \\ &\lesssim n^{\alpha-1} \mathbf{D}_{\phi}^{0} (1 + \exp(n^{-\epsilon})) + n^{-\alpha-2} + n^{-\alpha \wedge (1+\epsilon-\alpha)} \\ &\lesssim n^{\alpha-1} \mathbf{D}_{\phi}^{0} + n^{-\alpha \wedge \frac{\beta}{6}}. \end{aligned} $$

We point out that choosing $\alpha = 1/2$ and $\beta > 3$ allows for a convergence rate $O(n^{-1/2})$, which is similar to the previous case.

# 4 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.

## 4.1 Description of the simulation study

**Bregman divergence vs Euclidean projection** In order to assess the specificity and efficiency 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 definition is straightforward (see [22, 32, 18]). In what follows, we have used the recent algorithm introduced in [22] to project any points into the simplex of probability measure for which a fast implementation is available in [32] and for which a complexity bound of $O(m^2)$ is established in [18]. That being said, it is also shown in [18] 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 [27]. As reported in Table 1 of [18], the observed complexity of several projection methods, including the one of [18], is $O(m)$. We have chosen to compare SMD with PSGD with the projection method of [18] as indicated in Section 4 of [18], 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 (9) and (10). 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.
**Diversification and influence of $\lambda$** To undertake the influence of the penalty parameter $\lambda$, and in particular its effect on the composition of the optimal portfolio, we vary its value in a regular grid $(\lambda_{\min} + s \frac{\lambda_{\max} - \lambda_{\min}}{\Lambda})_{0 \le s \le \Lambda}$ defined with the help of a lower and upper value $(\lambda_{\min}, \lambda_{\max})$ and a number of points $\Lambda + 1$ in the grid. According to Proposition 1, it is expected that large values of $\lambda$ essentially favor low risk portfolio composition whereas small values of $\lambda$ allow for risky portfolio with higher expected returns.

**Efficient frontier construction** Since the value of $\lambda$ has to be determined from a practical point of view, we decided to use the Markowitz efficient frontier [33]. The efficient 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 [29] to build our Markowitz efficient frontier. The frontier is constructed as follow, for an integer $s \in [0, \Lambda]$, we compute with our SMD algorithm an approximation of the optimal portfolio associated to this penalty parameter $\lambda_s = \lambda_{\min} + s \frac{\lambda_{\max} - \lambda_{\min}}{\Lambda}$, i.e. the set of weights $(\hat{u}_{\lambda_s}, \hat{\theta}_{\lambda_s})$ that approximately minimizes $p_{\lambda_s}$:

$$ (\hat{u}_{\lambda_s}, \hat{\theta}_{\lambda_s}) = \arg \min_{(u, \theta) \in \Delta_m \times \mathbb{R}} p_{\lambda_s}. $$

For convenience, we denote in $p_{\lambda}(u, \theta) = -J(u) + \lambda \psi_{\alpha}(\theta, u)$ where $J(u) = - \sum_{i=1}^{m} u_i \mathbb{E}[Z_i]$. We then draw the frontier obtained ($CV@R_{\alpha}(\hat{u}_{\lambda_s}), J(\hat{u}_{\lambda_s})$) that may be computed on-line all along the iterations of the SMD algorithm with a Cesaro averaging strategy: if $(X_k)_{k \ge 0} = (u_k, \theta_k)_{k \ge 0}$ corresponds to a random sequence built with Algorithm 1 stopped at iteration $n$ and if $(\hat{Z}^k)_{k \ge 0}$ is the biased portfolio simulation sequence used all along Algorithm 1, then we may approximate the expected return $J(\hat{u}_{\lambda_s})$ and the risk $CV@R_{\alpha}(\hat{u}_{\lambda_s})$ by:

$$ \frac{1}{n} \sum_{k=1}^{n} \langle u_k, \hat{Z}^k \rangle \quad \text{and} \quad \frac{1}{n} \sum_{k=1}^{n} \left( \theta_k + \frac{1}{1-\alpha} \lfloor \langle \hat{Z}^k, u_k \rangle - \theta_k \rfloor_+ \right). $$

We emphasize that since this construction requires the optimization of the portfolio for each value of $\lambda$, it is mandatory to design an efficient 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. [45]) 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, \dots, 0)$, the Sharpe ratio criterion selects the portfolio that maximizes the ratio between the difference $J(\hat{u}_{\lambda_s}) - J(1_{CIR})$ and the risk measured by $CV@R_{\alpha}(\hat{u}_{\lambda_s})$.

## 4.2 Synthetic dataset

In this first set of simulations, we implement our SMD with purely artificial 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 financial application purposes. In a second stage in Section 4.3, we will consider correlated assets.

**Comparison between SMD and PSGD** Our first experiment consists in comparing the behaviour of the SMD and the PSGD. The value of $\lambda$ used here is $\lambda = 0.9$.

Figure 2 shows that there is no real evidence that SMD and PSGD lead to very different results in terms of the rate of convergence of $p_{\lambda}(X_k) \rightarrow \min(p_{\lambda})$, 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_{\lambda}$ may have an infinite set of minima, with a continuum of critical points. Even in this theoretically
<figure><img src="image_3.png" /><figcaption>Figure 1: Time evolution of the return of the discretized trajectories associated to the assets of the synthetic portfolio (CIR + GBM).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>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 different as represented above with SMD and PSGD.</figcaption></figure>

difficult situation, the trajectory of SMD (and the one of PSGD) almost surely converges as indicated in Theorem 4 and as shown in Figure 2.

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 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 amplified when $m$ increases since the projection algorithm requires more and more operations in larger dimensions.

<table><caption>Table 1: Ratio of computation with $n$ iterations on $m$ assets between PSGD and SMD (rounded at the closest integer).</caption><thead><tr><th></th><th>n = 10<sup>4</sup></th><th>n = 10<sup>5</sup></th></tr></thead><tbody><tr><td>m = 5</td><td>2</td><td>3</td></tr><tr><td>m = 20</td><td>3</td><td>5</td></tr><tr><td>m = 40</td><td>4</td><td>6</td></tr><tr><td>m = 80</td><td>5</td><td>7</td></tr></tbody></table>

All the more, as indicated in Figure 3 (the same conclusion can be drawn from Figure 2), the
<figure><img src="image_5.png" /><figcaption>Figure 3: 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_\lambda$ with $\lambda = 0.9$ and $\alpha = 0.05$.</figcaption></figure>

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) amplifies the benefits 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. [47]).

**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 sufficient to overcome the additional supplementary cost (in terms of seconds) generated by the Monte-Carlo step. 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.

**Effect of the penalty parameter $\lambda$** In Figure 5, we assess the influence of the penalty coefficient $\lambda$ involved in $p_\lambda$ on the optimal solutions obtained with the SMD. The experiments have been set up with 3 and 8 assets in the portfolio. The fluctuations 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 subfigure 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 $\lambda$, 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 $\lambda$ when essentially the CIR component gathers the essential weight of the portfolio. Finally, in the intermediary regime, the CV@R penalty generates a product diversification, which is recovered both on the left and on the right of Figure 5, which is also exactly the objective of
<figure><img src="image_6.png" /><figcaption>Figure 4: Evolution of the variance estimation (left) of the assets with the number of iterations of the SMD and of the MC-MD $p_λ$ with $λ = 0.95$ and $α = 0.05$. The composition of optimal portfolio (right) converge to a similar value as represented on the right with SMD and MC-MD but with a cost in computational time.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>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 subfigure.</figcaption></figure>

the CV@R penalty.

## 4.3 Real dataset

We briefly 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.

**Yahoo! dataset and parameters estimation** In this final study, we consider some stock prices that are daily recorded on the Yahoo! website https://finance.yahoo.com/lookup and that are freely available with the yahoo-finance API of Python [37].

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 [25] as it outperforms the discretized log-likelihood
maximization approach (for discretely observed CIR realizations). We refer to Equation 4.2 for the drift coefficients and Equation given at the end of Section 4.3 of [25].

Concerning the GBM, we estimate their coefficients $(\mu, \sigma$ and $\rho$) with the log-return series, which is a standard procedure described among other in Section 3.2.3 of [24] for example.

<figure><img src="image_8.png" /><figcaption>Figure 6: Sample of the ^FVX series and of the assets that composed our portfolio, from the Yahoo! database, between the 01-01-2014 and 01-01-2018 (2016 for the ^FVX time series).</figcaption></figure>

**Markowitz efficient frontier with CV@R** Once the parameters are estimated, we assume that their values are kept fixed 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 financial assets. We then use our SMD to estimate the Markowitz efficient frontier, adapted with the CV@R as a measure of risk of the portfolio. We represent the optimal portfolio in Figure 7 for two different levels of quantiles ($\alpha = 5\%$ and $\alpha = 1\%$).

<figure><img src="image_9.png" /><figcaption>Figure 7: Markowitz Efficient frontier. Left: $\alpha = 5\%$ and Right: $\alpha = 1\%$. Red point: optimal Sharpe ratio location. Blue points: location of the assets involved in the portfolio. Pink: risk free asset.</figcaption></figure>

We use the construction of our frontier to discretize 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 \propto n^{-3}$ appears to be reasonably 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 an optimal portfolio that really balances the CV@R and the ER.
# 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 1.** We observe that $J(u) = -\sum_{i=1}^{m} u_i E[Z_i]$ is a continuous convex function, defined on the convex simplex $\Delta_m$ and $P_M$ is a convex optimization problem defined with a constraint $g_M$ defined 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 [38] for this important remark, and to [5, 6] for seminal contributions on risk measurements in mathematical finance. The important consequence of such coherence is the convexity of the function that defines the constraint in $P_M$. We verify that $\forall \lambda \in (0, 1), \forall (u, v) \in \Delta_m$:

$$ 
\begin{aligned}
\mathrm{CV@R}_{\alpha}(\lambda u + (1 - \lambda)v) &\leq \mathrm{CV@R}_{\alpha}(\lambda u) + \mathrm{CV@R}_{\alpha}((1 - \lambda)v) \\
&= \lambda \mathrm{CV@R}_{\alpha}(u) + (1 - \lambda)\mathrm{CV@R}_{\alpha}(v).
\end{aligned}
$$

We are led to define the convex function $g_M$ that induces the constraint in $P_M$ by:

$$ g_M : u \mapsto \mathrm{CV@R}_{\alpha}(u) - M. $$

We shall observe that the collection of functions $(g_M)_{M \in \mathbb{R}}$ defines a family of convex constraints and $\exists M_0 \in \mathbb{R}$ such that for all $M > M_0$, the constraint $g_M$ is strictly feasible: there exists $u \in \Delta_m$ that is an interior point such that $g_M(u) < 0$. For a such $M > M_0$, we then define the Lagrangian function:

$$ \mathcal{L}_M(u, \lambda) := J(u) + \lambda g_M(u). $$

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 satisfied: the solution $u_M^*$ of the primal problem $P_M$ verifies:

$$ \exists \lambda_M^* \geq 0 \quad u_M^* = \arg\min_{u \in \Delta_m} \mathcal{L}_M(u, \lambda_M^*), \qquad (22) $$

and the complementary slackness is verified for the pair $(u_M^*, \lambda_M^*)$:

$$ \lambda_M^* g_M(u_M^*) = 0. \qquad (23) $$

We now establish the monotonous property. We consider a pair $(M_1, M_2)$ and consider $(u_1^*, \lambda_1^*)$ and $(u_2^*, \lambda_2^*)$ derived from the previous Lagrangian equation. We know that the complementary slackness condition holds for both $(u_1^*, \lambda_1^*)$ and $(u_2^*, \lambda_2^*)$. Equations (22) and (23) entail:

$$ 
\begin{aligned}
J(u_1^*) &= J(u_1^*) + \lambda_1^* g_{M_1}(u_1^*) = \mathcal{L}_{M_1}(u_1^*, \lambda_1^*) \\
&\leq \mathcal{L}_{M_1}(u_2^*, \lambda_1^*) \\
&= J(u_2^*) + \lambda_1^* g_{M_1}(u_2^*) \\
&= J(u_2^*) + \lambda_2^* g_{M_1}(u_2^*) + (\lambda_1^* - \lambda_2^*) g_{M_1}(u_2^*) \\
&= J(u_2^*) + \lambda_2^* [g_{M_2}(u_2^*) + M_2 - M_1] + (\lambda_1^* - \lambda_2^*) g_{M_1}(u_2^*) \\
&= \mathcal{L}_{M_2}(u_2^*, \lambda_2^*) + \lambda_2^* [M_2 - M_1] + (\lambda_1^* - \lambda_2^*) g_{M_1}(u_2^*) \\
&\leq \mathcal{L}_{M_2}(u_1^*, \lambda_2^*) + \lambda_2^* [M_2 - M_1] + (\lambda_1^* - \lambda_2^*) g_{M_1}(u_2^*) \\
&= J(u_1^*) + \lambda_2^* g_{M_2}(u_1^*) + \lambda_2^* (M_2 - M_1) + (\lambda_1^* - \lambda_2^*) [\mathrm{CV@R}_{\alpha}(u_2^*) - M_1] \\
&= J(u_1^*) + \lambda_2^* g_{M_1}(u_1^*) + (\lambda_1^* - \lambda_2^*) [g_{M_2}(u_2^*) + M_2 - M_1] \\
&\leq J(u_1^*) + \lambda_2^* g_{M_1}(u_1^*) + (\lambda_1^* - \lambda_2^*) (M_2 - M_1) + \lambda_1^* g_{M_2}(u_2^*), \quad (24)
\end{aligned}
$$

where we used in the last line Equation (23) at point $u_2^*$. Since $g_{M_1}(u_1^*) \leq 0$ and $g_{M_2}(u_2^*) \leq 0$ while $\lambda_1^* \geq 0$ and $\lambda_2^* \geq 0$, we deduce from Equation (24) that:

$$ J(u_1^*) \leq J(u_1^*) + (\lambda_1^* - \lambda_2^*)(M_2 - M_1). $$
It finally implies that

$$M_1 < M_2 \Rightarrow \lambda_1^* \ge \lambda_2^*.$$

We now consider the penalized criterion with $\lambda_1 > \lambda_2$, we verify that:

$$ 
\begin{aligned}
J(v_{\lambda_1}) + \lambda_1 \text{CV@R}_\alpha(v_{\lambda_1}) &\le J(v_{\lambda_2}) + \lambda_1 \text{CV@R}_\alpha(v_{\lambda_2}) \\
&= J(v_{\lambda_2}) + \lambda_2 \text{CV@R}_\alpha(v_{\lambda_2}) + (\lambda_1 - \lambda_2) \text{CV@R}_\alpha(v_{\lambda_2}) \\
&\le J(v_{\lambda_1}) + \lambda_2 \text{CV@R}_\alpha(v_{\lambda_1}) + (\lambda_1 - \lambda_2) \text{CV@R}_\alpha(v_{\lambda_2}) \\
&= J(v_{\lambda_1}) + \lambda_1 \text{CV@R}_\alpha(v_{\lambda_1}) + (\lambda_1 - \lambda_2) [\text{CV@R}_\alpha(v_{\lambda_2}) - \text{CV@R}_\alpha(v_{\lambda_1})].
\end{aligned}
$$

It implies that:

$$\lambda_1 > \lambda_2 \Rightarrow \text{CV@R}_\alpha(v_{\lambda_1}) \le \text{CV@R}_\alpha(v_{\lambda_2}).$$

Obviously, we also observe that $v_\lambda$ solves $\mathcal{P}_M$ for $M = \text{CV@R}_\alpha(v_\lambda)$, e.g.:

$$v_\lambda = \arg \min_{v \in S_m} \{J(v) : \text{CV@R}_\alpha(v) \le \text{CV@R}_\alpha(v_\lambda)\}.$$

□

# 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_\Phi$, whose statements are given below. We refer to [35] for further details.

**Lemma 9** (Three points lemma). For any triple of points $(x, y, z)$, one has:

$$D_\Phi(x, z) = D_\Phi(x, y) + D_\Phi(y, z) - \langle \nabla \Phi(z) - \nabla \Phi(y), x - y \rangle.$$

**Lemma 10** (Gradient of the Bregman divergence). For any pair of points $(x, y)$, one has:

$$\nabla_x D_\Phi(x, y) = \nabla \Phi(x) - \nabla \Phi(y).$$

With the help of these lemmas, we derive the proof of the almost sure convergence result stated in Theorem 4. The main issue generated by our biased simulation setting is that $\hat{g}$ involved in Equations (9), (10) and (11) acts non-linearly on $1_{\{\hat{g}^{k+1}, U_k\} \ge \theta_k}$ in our stochastic approximation term, which lead to a significant amount of theoretical difficulties.

*Proof of Theorem 4.* Let us recall that throughout the proof, the position of the algorithm is $X_k = (U_k, \theta_k)$ and that we denote by $x = (u, \theta)$ a generic point of the set $\mathcal{X} = \Delta_m \times \mathbb{R}$.

**Proof of i).** The proof is divided into three steps. We first 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_\Phi(x_\lambda^*, X_k)$, where $x_\lambda^*$ is a point that minimizes $p_\lambda$. Our starting point is the definition of $X_{k+1}$:

$$X_{k+1} = \arg \min_{x \in \mathcal{X}} \left\{ \langle \hat{g}_{k+1}, x - X_k \rangle + \frac{D_\Phi(x, X_k)}{\eta_{k+1}} \right\},$$
where we recall that $\hat{g}_{k+1}$ is the stochastic approximation of the sub-gradients:

$$ \begin{cases} \hat{g}_{k+1,1} = -\hat{Z}^{k+1} + \frac{\lambda}{1-\alpha}\hat{Z}^{k+1} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k}, \\ \hat{g}_{k+1,2} = 1 - \frac{1}{1-\alpha}\mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k}. \end{cases} $$

The first order condition entails that:

$$ \forall x \in \mathcal{X} \quad \eta_{k+1} \langle \hat{g}_{k+1}, x - x_{k+1} \rangle + \langle \nabla_x \mathcal{D}_\Phi(X_{k+1}, X_k), x - X_{k+1} \rangle \ge 0. $$

Using Lemma 10 on the second term, we deduce that:

$$ \forall x \in \mathcal{X} \quad \langle \nabla \Phi(X_{k+1}) - \nabla \Phi(X_k), x - X_{k+1} \rangle \ge \eta_{k+1} \langle \hat{g}_{k+1}, X_{k+1} - x \rangle. $$

We then apply Lemma 9 with $y = X_{k+1}$ and $z = X_k$ and obtain that:

$$ \eta_{k+1} \langle \hat{g}_{k+1}, X_{k+1} - x \rangle \le \langle \nabla \Phi(X_{k+1}) - \nabla \Phi(X_k), x - X_{k+1} \rangle \\ \le \mathcal{D}_{\Phi}(x, X_k) - \mathcal{D}_{\Phi}(x, X_{k+1}) - \mathcal{D}_{\Phi}(X_{k+1}, X_k) $$

Using the strong convexity of $\mathcal{D}_\Phi$ stated in Inequality (7), we deduce that:

$$ \eta_{k+1} \langle \hat{g}_{k+1}, X_{k+1} - x \rangle \le \mathcal{D}_{\Phi}(x, X_k) - \mathcal{D}_{\Phi}(x, X_{k+1}) - (\theta_{k+1} - \theta_k)^2 - \frac{1}{2} \|U_{k+1} - U_k\|^2, $$

which can be written as:

$$ \begin{aligned} \mathcal{D}_{\Phi}(x, X_{k+1}) &\le \mathcal{D}_{\Phi}(x, X_k) - (\theta_{k+1} - \theta_k)^2 - \frac{1}{2} \|U_{k+1} - U_k\|^2 - \eta_{k+1} \langle \hat{g}_{k+1}, X_{k+1} - x \rangle \\ &\le \mathcal{D}_{\Phi}(x, X_k) - (\theta_{k+1} - \theta_k)^2 - \frac{1}{2} \|U_{k+1} - U_k\|^2 \\ &\quad - \eta_{k+1} \langle \hat{g}_{k+1}, X_k - x \rangle - \eta_{k+1} \langle \hat{g}_{k+1}, X_{k+1} - X_k \rangle \end{aligned} \quad (25) $$

Let us first handle the last term of the right hand side:

$$ \eta_{k+1} \langle \hat{g}_{k+1}, X_{k+1} - X_k \rangle = \eta_{k+1} \hat{g}_{k+1,2} (\theta_{k+1} - \theta_k) + \eta_{k+1} \langle \hat{g}_{k+1,1}, U_{k+1} - U_k \rangle. $$

The Young inequality $|ab| \le \frac{a^2}{2c} + \frac{cb^2}{2}$ leads to:

$$ \begin{aligned} |\eta_{k+1} \hat{g}_{k+1,2} (\theta_{k+1} - \theta_k)| &\le \frac{\eta_{k+1}^2 \hat{g}_{k+1,2}^2}{4} + (\theta_{k+1} - \theta_k)^2 \\ &\le C_\alpha \eta_{k+1}^2 + (\theta_{k+1} - \theta_k)^2. \end{aligned} \quad (26) $$

where in the last line we use that associated with $|\hat{g}_{k+1,2}| \le \max(1, \frac{\alpha}{1-\alpha})$ (see Equation (11)). In the same way, we also have:

$$ |\eta_{k+1} \langle \hat{g}_{k+1,1}, U_{k+1} - U_k \rangle| \le \frac{\eta_{k+1}^2 \|\hat{g}_{k+1,1}\|^2}{2} + \frac{1}{2} \|U_{k+1} - U_k\|^2. \quad (27) $$

We use Equations (26) and (27) in Inequality (25) and obtain that:

$$ \mathcal{D}_{\Phi}(x, X_{k+1}) \le \mathcal{D}_{\Phi}(x, X_k) + \eta_{k+1}^2 [C_{\alpha} + \|\hat{g}_{k+1,1}\|^2] - \eta_{k+1} \langle \hat{g}_{k+1}, X_k - x \rangle. \quad (28) $$

We need to handle the biased drift that comes from the sampled random variables $Z_{k+1}$ involved in $\hat{g}_{k+1}$. We write:

$$ \begin{aligned} \hat{g}_{k+1} &= \nabla p_{\lambda}(X_k) + (\mathbb{E}[\hat{g}_{k+1} | \mathcal{F}_k] - \nabla p_{\lambda}(X_k)) + (\hat{g}_{k+1} - \mathbb{E}[\hat{g}_{k+1} | \mathcal{F}_k]) \\ &= \nabla p_{\lambda}(X_k) - \mathfrak{b}_{k+1} + \Delta M_{k+1}, \end{aligned} $$
with $\Delta M_{k+1} := \hat{g}_{k+1} - \mathbb{E}[\hat{g}_{k+1} | \mathcal{F}_k]$ is a martingale increment and $\mathfrak{b}_{k+1}$ stands for the bias:

$$ \mathfrak{b}_{k+1} := \nabla p_{\lambda}(U_k, \theta_k) - \mathbb{E}[\hat{g}_{k+1} | \mathcal{F}_k] \\ = \left( \mathbb{E}[\hat{Z}^{k+1} | \mathcal{F}_k] - \mathbb{E}[Z] + \frac{\lambda}{1-\alpha} \left( \mathbb{E}[Z \mathbf{1}_{\langle Z, U_k \rangle \ge \theta_k}] - \mathbb{E}[\hat{Z}^{k+1} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k} | \mathcal{F}_k] \right) \mathbb{P}(\langle Z, U_k \rangle \ge \theta_k) - \mathbb{P}(\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k | \mathcal{F}_k) \right). \quad (29) $$

Now, choosing $x = x_{\lambda}^{*}$ and using this decomposition into Equation (28), we obtain the main descent inequality, which will be the cornerstone of our analysis:

$$ \begin{aligned} \mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k+1}) \le \mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k}) & - \eta_{k+1} \langle \nabla p_{\lambda}(X_{k}), X_{k} - x_{\lambda}^{*} \rangle - \eta_{k+1} \langle \mathfrak{b}_{k+1}, X_{k} - x_{\lambda}^{*} \rangle \\ & - \eta_{k+1} \langle \Delta M_{k+1}, X_{k} - x_{\lambda}^{*} \rangle + \eta_{k+1}^{2} \left[ \frac{1}{4} + \|\hat{g}_{k+1,1}\|^{2} \right]. \end{aligned} \quad (30) $$

The Cauchy-Schwarz inequality entails:

$$ \eta_{k+1} |\langle \mathfrak{b}_{k+1}, X_k - x_{\lambda}^{*} \rangle| \le \eta_{k+1} \|\mathfrak{b}_{k+1}\| \|X_k - x_{\lambda}^{*}\|. \quad (31) $$

Now, the Young inequality and Equation (7) implies that:

$$ \|X_k - x_{\lambda}^{*}\| \le \frac{1 + \|X_k - x_{\lambda}^{*}\|^2}{2} \le \frac{1}{2} + 2\mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_k). \quad (32) $$

We then plug Equations (32) and (31) into (30) and obtain that:

$$ \begin{aligned} \mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k+1}) \le \mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k}) (1 + 2\eta_{k+1} \|\mathfrak{b}_{k+1}\|) & - \eta_{k+1} \langle \nabla p_{\lambda}(X_{k}), X_{k} - x_{\lambda}^{*} \rangle \\ & - \eta_{k+1} \langle \Delta M_{k+1}, X_{k} - x_{\lambda}^{*} \rangle + \eta_{k+1}^{2} \left[ \frac{1}{4} + \|\hat{g}_{k+1,1}\|^{2} \right] + \frac{1}{2} \eta_{k+1} \|\mathfrak{b}_{k+1}\|. \end{aligned} \quad (33) $$

**Step 2: Bias upper bound.** We are led to estimate the size of $\|\mathfrak{b}_{k+1}\|$. For this purpose, we introduce the Kolmogorov distance $d_{Kol}$, defined on real distributions and defined by:

$$ d_{Kol}(\mu, \nu) = \sup_{a \in \mathbb{R}} |\mu([ - \infty, a ]) - \nu([ - \infty, a ])|. $$

The Kolmogorov and Wasserstein distance for real distributions (see e.g. Theorem 3.1 of [17]) are related as follows:

$$ d_{Kol}(\mu, \nu) \le 2\sqrt{W_1(\mu, \nu)}. $$

Using the expression of $\mathfrak{b}_{k+1} = (\mathfrak{b}_{k+1,1}, \mathfrak{b}_{k+1,2})$ in (29), its second coordinates verifies:

$$ \mathbb{E}\|\mathfrak{b}_{k+1,2}\| \le d_{Kol} \left( \mathcal{L}(\langle \hat{Z}^{k+1}, U_k \rangle, \mathcal{L}(\langle Z, U_k \rangle)) \right) \le 2\sqrt{W_1 \left( \mathcal{L}(\langle \hat{Z}^{k+1}, U_k \rangle), \mathcal{L}(\langle Z, U_k \rangle) \right)}. $$

Using now the variational characterization of the $W_1$ distance on 1-Lipschitz function, and using that for any $k, U_k \in \Delta_m$, we have: $\|U_k\|^2 \le \|U_k\|_1^2 = 1$, we then deduce that:

$$ \sup_{u \in \Delta_m} W_1 \left( \mathcal{L}(\langle \hat{Z}^{k+1}, U_k \rangle), \mathcal{L}(\langle Z, U_k \rangle) \right) \le W_1(\mathcal{L}(\hat{Z}^{k+1}), \mathcal{L}(Z)). $$

We then obtain from Assumption (H$_1$) Equation (12) that:

$$ \mathbb{E}\|\mathfrak{b}_{k+1,2}\| \le 2\sqrt{\delta_{k+1}}. $$

We are led to study an upper bound of $\|\mathfrak{b}_{k+1,1}\|$. We observe that:

$$ \begin{aligned} \mathbb{E}\|\mathfrak{b}_{k+1,1}\| & \le \|\mathbb{E}[\hat{Z}^{k+1} - Z]\| + \frac{\lambda}{1-\alpha} \left\| \mathbb{E}[\hat{Z}^{k+1} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k} - Z \mathbf{1}_{\langle Z, U_k \rangle} \ge \theta_k] \right\| \\ & \le \delta_{k+1} + \frac{\lambda}{1-\alpha} \left\| \mathbb{E}[\hat{Z}^{k+1} \mathbf{1}_{\langle \hat{Z}^{k+1}, U_k \rangle \ge \theta_k} - Z \mathbf{1}_{\langle Z, U_k \rangle} \ge \theta_k] \right\|. \end{aligned} $$
We simply observe that in this case, the second term of the right hand side of the previous inequality is upper bounded by $v_{k+1}$ thanks to $(\mathbf{H}_1)$ Equation (13). We deduce that:

$$ \mathbb{E}[\|\mathbf{b}_{k+1}\|] \le \mathbb{E}[\|\mathbf{b}_{k+1,1}\|] + \mathbb{E}[\|\mathbf{b}_{k+1,2}\|] \le 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda}{1-\alpha}v_{k+1}. \quad (34) $$

*Step 3: Conclusion of the proof.* We apply the Robbins-Siegmund Lemma to Equation (33): we observe that the compactness of $\Delta_m$ and the existence of a second order moment for $Z$ leads to the existence of $C > 0$ such that:

$$ \mathbb{E}[\|\hat{g}_{k+1}\|^2] \le C < +\infty. \quad (35) $$

Finally, from Equation (34) and our assumptions on $(\delta_{k+1})_{k \ge 0}$ and $(v_{k+1})_{k \ge 0}$, one has:

$$ \sum_{k \ge 0} \eta_{k+1} \mathbb{E} \| \mathbf{b}_{k+1} \| < +\infty. \quad (36) $$

From Equation (36), the existence of a second order moment and since $\eta_{k+1} \langle \Delta M_{k+1}, X_k - x_\lambda^* \rangle$ is a $\mathcal{F}_k$ Martingale-increment, the Robbins-Siegmund Lemma (see e.g. [41]) implies that:

$$ \mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k}) \xrightarrow{k \to \infty} D_{\infty} \in L^{1} \quad \text{and} \quad \sum_{k \ge 0} \eta_{k+1} \langle \nabla p_{\lambda}(X_{k}), X_{k} - x_{\lambda}^{*} \rangle < +\infty \quad \text{a.s.} $$

This ends the proof of i).

◇

**Proof of ii)** This result is almost straightforward starting from i). We verify from $\mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k}) \xrightarrow{k \to \infty} D_{\infty} \in L^{1}$ that $(X_{k})_{k \ge 1}$ is an almost surely bounded sequence. Hence, almost surely the Cesaro average sequence $(\bar{X}_{k}^{\eta})_{k \ge 0}$ is a Cauchy sequence and therefore converges towards $\bar{X}_{\infty}^{\eta}$. We then use the convexity of $p_{\lambda}$ and observe that

$$ p_{\lambda}(\bar{X}_{k}^{\eta}) - \min(p_{\lambda}) \le \sum_{i=0}^{k} \left[ \frac{\eta_{i}}{\sum_{j=0}^{k} \eta_{j}} \right] p_{\lambda}(X_{k}) - \min(p_{\lambda}) = \frac{\sum_{i=0}^{k} \eta_{i} (p_{\lambda}(X_{k}) - \min(p_{\lambda}))}{\sum_{j=0}^{k} \eta_{j}} \to 0 \quad \text{a.s.} $$

The continuity of $p_{\lambda}$ then implies that almost surely:

$$ \bar{X}_{k}^{\eta} \to \bar{X}_{\infty}^{\eta} \quad \text{and} \quad p_{\lambda}(\bar{X}_{\infty}^{\eta}) = \min(p_{\lambda}). $$

◇

**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 \ge 0}$ towards a point of arg min($p_\lambda$).

*Step 1: continuous time trajectory.* We follow the roadmap of [14] that introduces the asymptotic pseudo-trajectories for stochastic algorithms, and [34] that adapts these tools to the mirror descent. We introduce the mirror map $Q$ and the Fenchel conjugate $\Phi^*$ defined on $\Delta_m \times \mathbb{R}$ by:

$$ Q(y) = \arg \max_{x \in \Delta_m \times \mathbb{R}} \langle y, x \rangle - \Phi(x) \quad \text{and} \quad \Phi^*(y) = \max_{x \in \Delta_m \times \mathbb{R}} \langle y, x \rangle - \Phi(x). $$

Note that the Fenchel coupling verifies that for any $(x, y) \in \{\Delta_m \times \mathbb{R}\}^2$:

$$ \Phi(x) + \Phi^*(y) - \langle y, x \rangle = D_{\Phi}(x, Q(y)). \quad (37) $$
The Mirror Descent ordinary differential equation is then defined as:

$$ \begin{cases} \dot{y} = -\nabla p_{\lambda}(x) \\ x = Q(y) \end{cases}. \qquad (38) $$

From Equation (7), we know that $\Phi$ is 1/2 strongly convex, which implies first that $Q$ is a 2-Lipschitz continuous function (see Proposition 2.2 of [34]). Second, $-\nabla p_{\lambda}$ is a continuous bounded function, which implies with the Cauchy-Arzela-Peano theorem that any solution of the Cauchy problem always exists.<fnref n="1" /> We consider any $x^*$ any minimizer of $p_{\lambda}$ and $(x_t, y_t)$ a solution of (38). Following the arguments of Nemirovski and Yudin [35], we define the Fenchel coupling as:

$$ V_{x^*}(t) = \Phi(x^*) + \Phi^*(y_t) - \langle y_t, x^* \rangle. $$

Note that using (37), $V_{x^*}(t) = D_{\phi}(x^*, y_t)$. The Fenchel inequality yields $V_{x^*}(t) \ge 0$ for any time $t$. We observe that:

$$ V'_{x^*}(t) = -\langle \nabla p_{\lambda}(x_t), x_t - x^* \rangle \le p_{\lambda}(x^*) - p_{\lambda}(x_t) \le 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_{\lambda}(\bar{x}_t)$ towards $\min(f)$ where $\bar{x}_t$ is the Cesaro average of the trajectory). From any trajectory $(x_t)_{t \ge 0}$, we consider $x_{\infty}$ any adherence point. If $x_{\infty}$ does not belong to $\min(p_{\lambda})$, then we can find a compact neighbourhood $\mathcal{V}_1$ of $x_{\infty}$ and an increasing sequence $(t_k)_{k \ge 0}$ such that $x_{t_k} \in \mathcal{V}_1$. Using the continuity of $p_{\lambda}$ we can find a second set $\mathcal{V}_2$ and a radius $r > 0$ such that

$$ \forall x \in \mathcal{V}_1: \quad B(x, r) \subset \mathcal{V}_2, $$

and for which we have:

$$ \forall x \in \mathcal{V}_2: \quad -\langle \nabla p_{\lambda}(x), x - x^* \rangle \le (\min(p_{\lambda}) - p_{\lambda}(x)) \le -a < 0. $$

We use that $Q$ is 2-Lipschitz to obtain that:

$$ \begin{aligned} \forall s > 0: \quad & \|x_{t_k+s} - x_{t_k}\| = \|Q(y_{t_k+s}) - Q(y_{t_k})\| \\ & \le 2\|y_{t_k+s} - y_{t_k}\| \\ & \le 2 \int_{t_k}^{t_k+s} \|\dot{y}_u\| du \\ & \le 2s \sup_{x \in \Delta_m \times \mathbb{R}} \|\nabla p_{\lambda}(x)\| \\ & \le 2s(\alpha + \mathbb{E}[\|Z\|]) \left(1 + \frac{\lambda}{1-\alpha}\right), \end{aligned} $$

where the last line comes from Equations (9) and (10). We can deduce that:

$$ \forall s \le s_r := \frac{r}{2(\alpha + \mathbb{E}[\|Z\|])(1 + \frac{\lambda}{1-\alpha})} : \quad x_{t_k+s} \in B(x_{t_k}, r) \subset \mathcal{V}_2. $$

Finally, we write that:

$$ V_{x^*}(t_k + s_r) - V_{x^*}(0) \le \sum_{j=0}^{k} \int_{t_j}^{t_j + s_r} V'_{x^*}(u) du \le -k \times (as_r) \to -\infty \quad \text{as} \quad k \to +\infty. $$

<footnote n="1">¹It is not clear that $\nabla p_{\lambda}$ is a Lipschitz continuous function, which prevents from the direct use of Theorem 2.4 of [34]. 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 compactifies the trajectories and permits to by-pass the uncertainty on this Lipschitz property.</footnote>
This contradicts the compactness of the trajectories and leads to a contradiction. We deduce that $x_{\infty} \in \arg\min(p_{\lambda})$ and that any adherence point of $(x_t)_{t \ge 0}$ belongs to $\arg\min(p_{\lambda})$.

Finally, we prove that the whole trajectory converges towards $x_{\infty}$. The construction of $V_{x^*}$ holds for any $x^*$ in $\arg\min\{p_{\lambda}\}$, in particular for $x^* = x_{\infty}$. For this particular choice, we verify that:

$$ \forall t \ge t_k \quad 0 \le V_{x_{\infty}}(t) \le V_{x_{\infty}}(t_k), $$

since $V$ is non-increasing and larger than 0.

Since by construction $x_{t_k}$ converges to $x^*$, then $D_{\Phi}(x_{\infty}, x_{t_k}) = V_{x^*}(t_k) \to 0$ as $k \to +\infty$, which in turn implies that $V_{x_{\infty}}(t) \to 0$ as $t \to +\infty$. We then conclude that the whole trajectory converges towards $x_{\infty}$.

**Step 2:** stochastic algorithm convergence. To derive the asymptotic behaviour of the SMD, we follow the argument of Proposition 4.1 in [14]. We consider $(X_k, Y_k)$ that satisfies the stochastic recursion

$$ \begin{cases} Y_{k+1} = Y_k - \eta_{k+1} \hat{g}_{k+1} \\ X_{k+1} = Q(Y_{k+1}) \end{cases} $$

We then define the continuous time affine interpolation of $(X_k, Y_k)_{k \ge 1}$ using the natural time-scale $\tau_0 = 0$ and $\tau_{k+1} = \tau_k + \eta_{k+1}$ and introduce $(\hat{X}_t, \hat{Y}_t)_{t \ge 0}$ the stochastic process defined as

$$ \forall t \in [\tau_k, \tau_{k+1}] \quad (\hat{X}_t, \hat{Y}_t) = (X_k, Y_k) + \frac{t - \tau_k}{\tau_{k+1} - \tau_k} ((X_{k+1}, Y_{k+1}) - (X_k, Y_k)). $$

We still use the decomposition derived from (29):

$$ \hat{g}_{k+1} = \nabla p_{\lambda}(X_k) + \Delta M_{k+1} + \mathfrak{b}_{k+1}, $$

where the martingale increment $\Delta M_{k+1}$ satisfies (see (35)):

$$ \sup_{k \ge 0} \mathbb{E}[\|\Delta M_{k+1}\|^2] < +\infty, $$

and the bias term satisfies under Equation (34) and assumption of iii) that:

$$ \sum_{k \ge 0} \mathbb{E}[\|\mathfrak{b}_{k+1}\|] \le \sum_{k \ge 0} \left( 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda}{1-\alpha} v_{k+1} \right) \le +\infty. $$

Then, the Markov inequality and the Borel-Cantelli lemma yields

$$ \mathfrak{b}_k \to 0 \quad \text{a.s. as } k \to +\infty. $$

We then use Remark 4.5 and Proposition 4.2 of [14] and deduce that Assumption A1 of Proposition 4.1 [14] holds. In the meantime, Assumption A2 of Proposition 4.1 [14] is valid from the proof of i) and the convergence of $D_{\Phi}(x_{\lambda}^*, X_k)$. Since $Q$ is a 2-Lipschitz function, we then deduce that $(\hat{X}_t, \hat{Y}_t)_{t \ge 0}$ is an asymptotic pseudo-trajectory of the deterministic mirror descent o.d.e. stated in (38). In particular, $(X_t)_{t \ge 0}$ converges using Step 1 of iii), which concludes the proof of iii). $\square$

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

**Step 1:** control of $\mathbb{E}[D_{\Phi}(x_{\lambda}^*, X_k)]$. For the sake of convenience, we denote by:

$$ \mathbb{E}D_{\Phi}(x_{\lambda}^*, X_k) = \mathbf{D}_{\Phi}^k. $$
Taking the expectation in (33) and using $\langle \nabla p_{\lambda}(X_k), X_k - x_{\lambda}^{*} \rangle > 0$, we get:

$$ 
\begin{aligned}
\mathbf{D}_{\Phi}^{k+1} &\le \mathbf{D}_{\Phi}^{k} + 2\eta_{k+1} \mathbb{E} (\mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_k) \|\mathfrak{b}_{k+1}\|) + \eta_{k+1}^{2} \left[ \frac{1}{4} + \mathbb{E} \|\hat{g}_{k+1,1}\|^{2} \right] + \frac{1}{2} \eta_{k+1} \mathbb{E} \|\mathfrak{b}_{k+1}\| \\
&\le \mathbf{D}_{\Phi}^{k} \left( 1 + 2\eta_{k+1} \left[ 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda\upsilon_{k+1}}{1-\alpha} \right] \right) \\
&\quad + \eta_{k+1}^{2} \left[ \frac{1}{4} + \mathbb{E} \|\hat{g}_{k+1,1}\|^{2} \right] + \frac{1}{2} \eta_{k+1} \left( 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda\upsilon_{k+1}}{1-\alpha} \right),
\end{aligned}
$$

where we used (34) and a conditional expectation argument. The inequality can be written as:

$$ \mathbf{D}_{\Phi}^{k+1} \le \mathbf{D}_{\Phi}^{k}(1 + a_{k+1}) + b_{k+1}, $$

using that $\mathbb{E}\|\hat{g}_{k+1,1}\|^2 \le M$, a $C > 0$ exists such that:

$$ \left\{ \begin{aligned} a_{k+1} &= 2\eta_{k+1} \left( 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda\upsilon_{k+1}}{1-\alpha} \right) \\ b_{k+1} &= C \left( \eta_{k+1}^2 + \eta_{k+1} \left( 2\sqrt{\delta_{k+1}} + \delta_{k+1} + \frac{\lambda\upsilon_{k+1}}{1-\alpha} \right) \right). \end{aligned} \right. \quad (39) $$

Now, denote:

$$ \mathcal{V}_k = \frac{\mathbf{D}_{\Phi}^k}{\prod_{i=1}^k (1 + a_i)} - \sum_{j=0}^k \frac{b_j}{\prod_{i=1}^j (1 + a_i)}, $$

we verify that:

$$ 
\begin{aligned}
\mathcal{V}_{k+1} &= \frac{\mathbf{D}_{\Phi}^{k+1}}{\prod_{i=1}^{k+1} (1 + a_i)} - \sum_{j=1}^{k+1} \frac{b_j}{\prod_{i=1}^{j} (1 + a_i)} \\
&\le \frac{\mathbf{D}_{\Phi}^{k}(1 + a_{k+1}) + b_{k+1}}{\prod_{i=1}^{k+1} (1 + a_i)} - \sum_{j=1}^{k+1} \frac{b_j}{\prod_{i=1}^{j} (1 + a_i)} \\
&\le \frac{\mathbf{D}_{\Phi}^{k}}{\prod_{i=1}^{k}(1 + a_i)} + \frac{b_{k+1}}{\prod_{i=1}^{k+1}(1 + a_i)} - \sum_{j=1}^{k+1} \frac{b_j}{\prod_{i=1}^{j}(1 + a_i)} \\
&\le \frac{\mathbf{D}_{\Phi}^{k}}{\prod_{i=1}^{k}(1 + a_i)} - \sum_{j=1}^{k} \frac{b_j}{\prod_{i=1}^{j}(1 + a_i)} = \mathcal{V}_k.
\end{aligned}
$$

From this inequality, we deduce that:

$$ \frac{\mathbf{D}_{\Phi}^{k}}{\prod_{i=1}^{k}(1 + a_i)} - \sum_{j=1}^{k} \frac{b_j}{\prod_{i=1}^{j}(1 + a_i)} \le \mathcal{V}_0 = \mathbf{D}_{\Phi}^{0}. $$

Let us remark that:

$$ 
\begin{aligned}
\mathbf{D}_{\Phi}^{0} &= D_{\Phi}(X_0, x_{\lambda}^{*}) = \frac{(\theta_0 - \theta_{\lambda}^{*})^2}{2} + D_{\varphi}(U_0, u_{\lambda}^{*}) \\
&= \frac{(\theta_0 - \mathrm{V}\mathrm@{R}_{\alpha}(u_{\lambda}^{*}))^2}{2} + D_{\varphi}(U_0, u_{\lambda}^{*}) \\
&\le \frac{(\theta_0 - \mathrm{V}\mathrm@{R}_{\alpha}(u_{\lambda}^{*}))^2}{2} + \log m := \{\Delta_{\Phi}^{0}\}^2, \quad (40)
\end{aligned}
$$

where we used the complementary slackness condition (see Equation (23) in the Appendix A) and the upper bound on the entropic Bregman divergence used in [30]. We finally deduce that:

$$ \mathbf{D}_{\Phi}^{k} \le \mathbf{D}_{\Phi}^{0} \prod_{i=1}^{k} (1 + a_i) + \left( \sum_{j=1}^{k} \frac{b_j}{\prod_{i=1}^{j} (1 + a_i)} \right) \prod_{i=1}^{k} (1 + a_i). \quad (41) $$
**Step 2: non-asymptotic control of the position of the SMD.** Going back again to (33), we can now use the convergence of $\mathcal{D}_{\Phi}(x_{\lambda}^{*}, X_{k}) \to 0$ to get:

$$ \eta_{k+1} \langle \nabla p(X_k), X_k - x_* \rangle \le \mathcal{D}_{\Phi}(x_*, X_k) (1 + 2\eta_{k+1} \|\mathfrak{b}_{k+1}\|) - \mathcal{D}_{\Phi}(x_*, X_{k+1}) \\ + \eta_{k+1} \langle \mathbb{E}[\hat{g}_{k+1} | \mathcal{F}_k] - \hat{g}_{k+1}, X_k - x_* \rangle + \eta_{k+1}^2 \left[ \frac{1}{4} + \|\hat{g}_{k+1,1}\|^2 \right] + \frac{1}{2} \eta_{k+1} \|\mathfrak{b}_{k+1}\|. $$

From the convexity of $p_{\lambda}$, we deduce that $\langle \nabla p_{\lambda}(X_k), X_k - x_{\lambda}^{*} \rangle \ge p_{\lambda}(X_k) - p_{\lambda}(x_{\lambda}^{*})$, and taking the expectation in both sides, we deduce that:

$$ \eta_{k+1} [\mathbb{E}p_{\lambda}(X_k) - p_{\lambda}(x_{\lambda}^{*})] \le (\mathbf{D}_{\Phi}^{k} - \mathbf{D}_{\Phi}^{k+1}) + a_{k+1}\mathbf{D}_{\Phi}^{k} + b_{k+1}. $$

Summing these inequalities from $k=0$ to $n$, we get:

$$ \begin{aligned} \sum_{k=0}^{n} \eta_{k+1} (\mathbb{E}(p_{\lambda}(X_k)) - p_{\lambda}(x_{\lambda}^{*})) &\le \sum_{k=0}^{n} (\mathbf{D}_{\Phi}^{k} - \mathbf{D}_{\Phi}^{k+1}) + \sum_{k=0}^{n} a_{k+1} \mathbf{D}_{\Phi}^{k} + b_{k+1} \\ &\le (\mathbf{D}_{\Phi}^{0} - \mathbf{D}_{\Phi}^{n+1}) + \sum_{k=0}^{n} a_{k+1} \mathbf{D}_{\Phi}^{k} + b_{k+1} \\ &\le \mathbf{D}_{\Phi}^{0} + \sum_{k=0}^{n} a_{k+1} \mathbf{D}_{\Phi}^{k} + b_{k+1}. \end{aligned} $$

Following the convexity argument of ii) Theorem 4 and the definition of $\bar{X}_{n}^{\eta}$ in(14), we deduce that:

$$ \begin{aligned} \mathbb{E}[p_{\lambda}(\bar{X}_{n}^{\eta})] - p_{\lambda}(x_{\lambda}^{*}) &\le \sum_{k=0}^{n} \frac{\eta_{k+1}}{\sum_{j=0}^{n} \eta_{j+1}} [\mathbb{E}(p_{\lambda}(X_k)) - p_{\lambda}(x_{\lambda}^{*})] \\ &\le \frac{1}{\sum_{j=0}^{n} \eta_{j+1}} \left( \mathbf{D}_{\Phi}^{0} + \sum_{k=0}^{n} [a_{k+1}\mathbf{D}_{\Phi}^{k} + b_{k+1}] \right), \end{aligned} \quad (42) $$

which ends the proof.

*Proof of Corollary 6.* We now consider fixed step sequences for the algorithm and the discretization stopped at a horizon $n$. More precisely, we assume $\eta_{k+1} = \eta > 0$ for $0 \le k \le n$ and that (34) translates:

$$ \mathbb{E}[\|\mathfrak{b}_{k+1}\|] \le \omega, \quad \forall 0 \le k \le n. $$

We aim at choosing constant values for step sequences adapted to the finite horizon $n$.

Within this framework the step sequences $(a_{k+1})$ and $(b_{k+1})$ defined in (39) become

$$ a_{k+1} = 2\eta\omega \quad \text{and} \quad b_{k+1} = C(\eta^2 + \eta\omega). $$

Let us now see how both $\eta$ and $\omega$ should depend on $n$. We these choices, equation (41) becomes:

$$ \mathbf{D}_{\Phi}^{k} \le \mathbf{D}_{\Phi}^{0}(1 + 2\eta\omega)^{k} + C(\eta^2 + \eta\omega) \frac{(1 + 2\eta\omega)^{k} - 1}{2\eta\omega}. $$
Equation (42) then yields :

$$
\begin{align*}
& \mathbb{E}p_{\lambda}(\bar{X}_{n}^{\eta}) - p_{\lambda}(x_{\lambda}^{*}) \\
\le {} & \frac{1}{(n+1)\eta} \left( \mathbf{D}_{\Phi}^{0} + \sum_{k=0}^{n} [2\eta\omega \mathbf{D}_{\Phi}^{k} + C\eta(\eta+\omega)] \right) \\
\le {} & \frac{1}{(n+1)\eta} \left( \mathbf{D}_{\Phi}^{0} \left( 1 + 2\eta\omega \sum_{k=0}^{n} (1+2\eta\omega)^k \right) + C\eta(\eta+\omega) \left( n+1 + \sum_{k=0}^{n} [(1+2\eta\omega)^{k+1} - 1] \right) \right) \\
\le {} & \frac{\mathbf{D}_{\Phi}^{0}}{(n+1)\eta} \left( 1 + 2\eta\omega \frac{(1+2\eta\omega)^{n+1}-1}{2\eta\omega} \right) + \frac{C\eta(\eta+\omega)}{\eta(n+1)} \left( \frac{(1+2\eta\omega)^{n+1}-1}{2\eta\omega} + n+1 \right) \\
\le {} & \frac{\mathbf{D}_{\Phi}^{0} e^{2\eta\omega(n+1)}}{(n+1)\eta} + C(\eta+\omega) \left[ \frac{e^{2\eta\omega(n+1)}}{2\eta\omega(n+1)} + 1 \right], \tag{43}
\end{align*}
$$

where we used in the last inequality that $1+x \le e^x$. We then choose $(\eta, \omega)$ such that:

$$
\eta = \frac{\sqrt{\mathbf{D}_{\Phi}^{0}}}{2\sqrt{n+1}} \quad \text{and} \quad \omega = \frac{1}{\sqrt{n+1}\sqrt{\mathbf{D}_{\Phi}^{0}}}
$$

and obtain, using the fact that $\eta\omega(n+1) \le 1/2$, that a large enough $C$ exists such that

$$
\mathbb{E}p_{\lambda}(\hat{X}_{n}^{\eta}) - p_{\lambda}(x_{\lambda}^{*}) \le C \frac{\sqrt{\mathbf{D}_{\Phi}^{0}}}{\sqrt{n+1}}.
$$

To conclude, we plug-in the upper bound (40) for $\sqrt{\mathbf{D}_{\Phi}^{0}}$ in the previous equations.
Starting from (43), we also observe that the choice $\eta = \omega = 1/\sqrt{n+1}$ yields:

$$
\mathbb{E}p_{\lambda}(\hat{X}_{n}^{\eta}) - p_{\lambda}(x_{\lambda}^{*}) \le C \frac{\mathbf{D}_{\Phi}^{0}}{\sqrt{n+1}}.
$$

□

C Proofs associated with the simulation of the portfolio

The main goal of this section is to prove Proposition 7 and Proposition 8. To this aim, we first need
technical properties on the drift-implicit Euler Scheme coupled with the Riemann integral approxima-
tion 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 satisfied by the CIR process and on some important objects that are related to the CIR introduced in Equation (15).

**Proposition 11.** $(r_t)_{t \ge 0}$ is a Markov process and satisfies:

i) If $2ab > \sigma_0^2$, then $(r_t)_{t \ge 0}$ remains strictly positive almost surely.

ii) The (conditional) expectation and variance are given by:

$$
\mathbb{E}[r_{t+s}|r_s] = r_s e^{-at} + b(1 - e^{-at}), \quad (44)
$$

and

$$
\mathbb{V}(r_{t+s}|r_s) = r_s \frac{\sigma_0^2}{a} (e^{-at} - e^{-2at}) + \frac{b\sigma_0^2}{2a} (1 - e^{-at})^2. \quad (45)
$$
iii) The CIR possesses an exponential integrability:

$$ \forall \lambda < \frac{a^2}{2\sigma_0^2}, \forall t \ge 0: \quad \mathbb{E}\left[\exp\left(\lambda \int_0^t r_s ds\right)\right] \le C(\lambda, t) < +\infty. $$

The two first items of the previous proposition are classical, and may be found in [24]. The exponential integrability of the CIR process is more recent and may be traced back to [20] (Proposition 3.2). We emphasize that within realistic situations, the values of $a$ and $\sigma$ for the CIR process are in general linked within a reasonable ratio of $5 \le \frac{a}{\sigma} \le 10$, so that the limiting support of $\lambda$ is generally significantly larger than 1. Note also that this last result holds at any horizon time but it may be explicit for $t = 1$ with a larger size of $\lambda$.

Proposition 11 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 [3] proved in [21], that provides an upper bound of the weak error associated to the implicit Euler scheme(20). In the previous reference, a strong error rate of order 1 is also provided with stronger assumptions on the parameters.

**Proposition 12** (Theorem 1 in [3]). Assume that $2ab > \sigma_0^2$, then for any $p \in [1, \frac{2ab}{\sigma_0^2})$, a constant $K_p$ exists such that:

$$ \left( \mathbb{E} \left[ \max_{k \in \{0, \dots, N\}} |\hat{r}_{kh} - r_{kh}|^p \right] \right)^{1/p} \le K_p \sqrt{h}. $$

**Control of the error in the approximation of the integral of the CIR** Recall the definition of $\Delta_h$ in (53)

$$ \Delta_h = \int_0^1 r_s ds - \frac{1}{N} \sum_{k=0}^{N-1} \hat{r}_{kh}. $$

Let us now prove two technical lemmas first on the estimation of $\mathbb{E}(\Delta_h^2)$ then on exponential moments.

**Lemma 13.** Assume that the CIR parameters satisfy $ab > \sigma^2$, then there exists a constant $C$ such that, for any $h = N^{-1} > 0$,

$$ \mathbb{E}(|\Delta_h|^2) \le Ch. $$

*Proof of Lemma 13.* We bound $\Delta_h^2$ using first the Jensen inequality and then the Cauchy-Schwarz inequality, we obtain that:

$$ \Delta_h^2 = \left[ \frac{1}{N} N \sum_{k=0}^{N-1} \int_{kh}^{(k+1)h} (r_s - \hat{r}_{kh}) ds \right]^2 \\ \le \frac{1}{N} \sum_{k=1}^{N} \left[ N \int_{kh}^{(k+1)h} (r_s - \hat{r}_{kh}) ds \right]^2 \\ \le \sum_{k=1}^{N} \int_{kh}^{(k+1)h} (r_s - \hat{r}_{kh})^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:

$$ \begin{aligned} \mathbb{E}(\Delta_h^2) & \le 2 \sum_{k=1}^{N} \left[ \int_{(k-1)h}^{kh} \mathbb{E}(r_s - r_{kh})^2 + \mathbb{E}(r_{kh} - \hat{r}_{kh})^2 \right] ds \\ & \le 2 \sum_{k=1}^{N} \int_{(k-1)h}^{kh} \mathbb{E}(r_s - r_{kh})^2 ds + 2h \sum_{k=0}^{N-1} \mathbb{E}(r_{kh} - \hat{r}_{kh})^2. \end{aligned} \quad (46) $$
Let us start with the first 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 11. Thus, for any $s = kh + u$ with $0 \le u \le h$, the bias-variance decomposition associated with Equations (44) and (45) leads to:

$$ 
\begin{aligned}
\mathbb{E}[(r_s - r_{kh})^2 | r_{kh}] &= (\mathbb{E}(r_s | r_{kh}) - r_{kh})^2 + \mathbb{V}(r_s | r_{kh}) \\
&= (b - r_{kh})^2(1 - e^{-au})^2 + \frac{\sigma^2}{a}r_{kh}e^{-as}(1 - e^{-au}) + b\frac{\sigma^2}{2a}(1 - e^{-au})^2.
\end{aligned}
$$

We then use again Equation (44) to obtain a recursive expression of $\mathbb{E}(r_{kh})$:

$$ 
\begin{aligned}
\mathbb{E}(r_{(k+1)h}) &= \mathbb{E}(\mathbb{E}(r_{(k+1)h}|r_{kh})) = \mathbb{E}(r_{kh})e^{-ah} + b(1 - e^{-ah}) \\
&= e^{-a(k+1)h}(\mathbb{E}(r_0) - b) + b.
\end{aligned}
$$

Therefore

$$ \sum_{k=1}^{N} \mathbb{E}(r_{kh}) = bN + (\mathbb{E}(r_0) - b) \sum_{k=1}^{N} e^{-akh} = bN + (\mathbb{E}(r_0) - b) \frac{1 - e^{-a}}{1 - e^{-ah}}. $$

From Equation (45) and a conditional expectation argument we deduce that:

$$ (\mathbb{E}(b - r_{(k+1)h})^2) = \mathbb{E}(r_{kh}) \frac{\sigma^2}{a} (e^{-ah} - e^{-2ah}) + \frac{b\sigma^2}{2a} (1 - e^{-ah})^2. $$

We then obtain that:

$$ \sum_{k=1}^{N} \int_{(k-1)h}^{kh} \mathbb{E}(r_s - r_{kh})^2 = \sum_{k=1}^{N} (\mathbb{E}(b - r_{kh})^2) \int_{0}^{h} (1 - e^{-au})^2 du + \frac{\sigma^2}{a} \sum_{k=1}^{N} (\mathbb{E}r_{kh}) \int_{0}^{h} e^{-au}(1 - e^{-au}) du \\ + b \frac{\sigma^2}{2a} N \int_{0}^{h} (1 - e^{-au})^2 du. $$

We then compute that:

$$ \int_{0}^{h} (1 - e^{-au})^2 du = \frac{a^2 h^3}{3} + O_a(h^4) \quad \text{and} \quad \int_{0}^{h} e^{-au}(1 - e^{-au}) du = \frac{ah^2}{2} - \frac{a^2 h^3}{2} + O_a(h^4). $$

Finally we deduce that for some $C > 0$

$$ \sum_{k=1}^{N} \int_{(k-1)h}^{kh} \mathbb{E}(r_s - r_{kh})^2 = Ch + O_{a,b}(h^2). \quad (47) $$

We now turn to the second term above, that we estimate using the strong error rate of Proposition 12 for $p = 2$. Assuming $ab > \sigma_0^2$, we get:

$$ \mathbb{E}(r_{kh} - \hat{r}_{kh})^2 \le \mathbb{E} \left[ \max_{1 \le k \le N} (r_{kh} - \hat{r}_{kh})^2 \right] \le Ch. $$

Since this bound is uniform in $k$, we deduce that:

$$ \sum_{k=1}^{N} \mathbb{E}(r_{kh} - \hat{r}_{kh})^2 \le C. \quad (48) $$

Combining estimates (47) and (48) in equation (46) gives the conclusion. □
**Lemma 14.** For any choice of $q \le 2$, if $a > 2\sqrt{2}\sigma_0$ then:

$$ \sup_N \mathbb{E}[e^{q|\Delta_h|}] < +\infty. \quad (49) $$

*Proof of Lemma 14.* We use the Cauchy-Schwarz inequality to upper-bound the expectation:

$$ \begin{aligned} \mathbb{E}[e^{q|\Delta_h|}] &= \mathbb{E}\left[e^{q\int_0^1 r_s ds - \frac{1}{N}\sum_{k=1}^N \hat{r}_{kh}}\right] \\ &\le \mathbb{E}\left[e^{q\int_0^1 r_s ds} e^{\frac{1}{N}\sum_{k=1}^N \hat{r}_{kh}}\right] \\ &\le \mathbb{E}\left(e^{2q\int_0^1 r_s ds}\right)^{\frac{1}{2}} \times \mathbb{E}\left(e^{\frac{2q}{N}\sum_{k=1}^N \hat{r}_{kh}}\right)^{\frac{1}{2}}. \end{aligned} $$

We apply iii) of Proposition 11 and observe that since $a > 2\sqrt{2}\sigma_0$, then $2q \le 4 \le \frac{a^2}{2\sigma_0^2}$ and:

$$ \mathbb{E}\left(e^{2q \int_0^1 r_s ds}\right)^{\frac{1}{2}} < +\infty. $$

For the second part, we rely on the recursive expression for $\hat{r}_{kh}$ stated in (12), that we bound using:

$$ \forall \epsilon > 0 \quad (u + v)^2 = u^2 + v^2 + 2uv \le u^2(1 + 1/\epsilon) + v^2(1 + \epsilon). $$

Thus

$$ \begin{aligned} \hat{r}_{(k+1)h} &= \left( \frac{\sqrt{\hat{r}_{kh}} + \frac{\sigma_0}{2} \Delta B_0^{(k)}}{2(1 + \frac{ah}{2})} + \sqrt{\frac{\left(\sqrt{\hat{r}_{kh}} + \frac{\sigma_0}{2} \Delta B_0^{(k)}\right)^2}{4(1 + \frac{ah}{2})^2} + \frac{(4ab - \sigma_0^2)h}{8(1 + \frac{ah}{2})}} \right)^2 \\ &\le (1 + \frac{\epsilon}{2}) \left( \frac{\sqrt{\hat{r}_{kh}} + \frac{\sigma_0}{2} \Delta B_0^{(k)}}{2(1 + \frac{ah}{2})} \right)^2 + (1 + \frac{1}{2\epsilon}) \left[ \frac{\left(\sqrt{\hat{r}_{kh}} + \frac{\sigma_0}{2} \Delta B_0^{(k)}\right)^2}{4(1 + \frac{ah}{2})^2} + \frac{(4ab - \sigma_0^2)h}{8(1 + \frac{ah}{2})} \right] \\ &\le (2 + \frac{\epsilon}{2} + \frac{1}{2\epsilon}) \frac{\left(\sqrt{\hat{r}_{kh}} + \frac{\sigma_0}{2} \Delta B_0^{(k)}\right)^2}{4(1 + \frac{ah}{2})^2} + (1 + \frac{1}{2\epsilon}) \frac{(4ab - \sigma_0^2)h}{8(1 + \frac{ah}{2})} \\ &\le \frac{(2 + \frac{\epsilon}{2} + \frac{1}{2\epsilon})(1 + 2\tau)}{4(1 + \frac{ah}{2})^2} \hat{r}_{kh} + \frac{(2 + \frac{\epsilon}{2} + \frac{1}{2\epsilon})(1 + \frac{1}{2\tau})}{4(1 + \frac{ah}{2})^2} \frac{\sigma_0^2}{4} \{\Delta B_0^{(k)}\}^2 + (1 + \frac{1}{2\epsilon}) \frac{(4ab - \sigma_0^2)h}{8(1 + \frac{ah}{2})}. \end{aligned} $$

Choosing $\epsilon = 1$ and $\tau = \frac{1}{7} < \frac{1}{6}$, we deduce that:

$$ \hat{r}_{(k+1)h} \le \frac{1}{(1 + \frac{ah}{2})^2} \hat{r}_{kh} + \frac{27}{32(1 + \frac{ah}{2})^2} \frac{\sigma_0^2}{4} \{\Delta B_0^{(k)}\}^2 + \frac{3(4ab - \sigma_0^2)h}{16(1 + \frac{ah}{2})} $$

Using the independence and stationarity of the Brownian motion, we set $\xi_{k+1}$ a standard Gaussian random variable such that $\xi_{k+1}\sqrt{h} = \Delta B_0^{(k)}$, therefore:

$$ \hat{r}_{(k+1)h} \le \hat{r}_{kh} + \alpha \xi_{k+1}^2 + \beta $$

where

$$ \alpha = \frac{27}{32(1 + \frac{ah}{2})^2} \frac{\sigma_0^2}{4} h \quad \text{and} \quad \beta = \frac{3(4ab - \sigma_0^2)h}{16(1 + \frac{ah}{2})}. $$

A straightforward recursion yields:

$$ \hat{r}_{(k+1)h} \le \hat{r}_0 + \alpha \sum_{j=0}^{k} \xi_{j+1}^2 + \beta k. $$
In particular, we can use this to upper bound the sum with:

$$ \frac{1}{N} \sum_{k=1}^{N} \hat{r}_{kh} \leq \hat{r}_{0} + \alpha \sum_{j=1}^{N} \xi_{j+1}^{2} + \beta N. $$

We can now bound the Laplace transform, using the independence of the $\xi_j$ and the Laplace transform of a $\chi^2$ distribution $\mathbb{E}(e^{\lambda\xi^2}) = (1 - 2\lambda)^{-1/2}, \forall \lambda < 1/2$.

$$ \begin{aligned} \mathbb{E} \left( e^{\frac{2q}{N} \sum_{k=1}^{N} \hat{r}_{kh}} \right) &\le e^{2q \hat{r}_0 + 2qN\beta} \prod_{j=1}^{N} \mathbb{E} \left[ e^{2q\alpha\xi_j^2} \right] \\ &\le e^{2q \hat{r}_0 + 2qN\beta} \left( \frac{1}{1 - 4q\alpha N} \right)^{1/2} \\ &\le e^{2q \hat{r}_0 + 2qN\beta} e^{2q\alpha N}, \end{aligned} $$

where the last line derives from $1/(1-t) \le e^t$. From our choice of $\alpha$ and $\beta$, we observe that $\alpha N < +\infty$ and $\beta N < +\infty$. We then conclude that $\mathbb{E}(e^{2q\frac{1}{N}\sum_{k=1}^{N} \hat{r}_{kh}})^{\frac{1}{2}} < +\infty$. $\square$

## C.2 Results on the portfolio $Z(t)$

**Proposition 15 (Density of the portfolio).** Assume that $\Sigma$ is positive definite, then the distribution of $S = (Y_1, S_1, \dots, 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, \dots, z_m) \mapsto \frac{M(z_1, \dots, z_m)}{(\sigma_1 \dots \sigma_{m'}) \det(L) z_1 \dots z_m}, $$

where $M$ is a bounded continuous function, whose bound is independent from $m$ and $\Sigma$.

*Proof.* Let us argue first in the case of uncorrelated Brownian motions in equations (16), i.e. when $\Sigma$ is the identity matrix. The correlated case will again be dealt with using Cholesky decomposition. Note that

$$ Y_1(t) = Y_1(0) \exp \left( \int_0^t r_s \right). $$

The density of the integral of the CIR process can be derived from Pitman and Yor [39], 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 [26]. Next, since the assets $S_1(t), \dots, 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, \dots, 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 $(\int_0^1 r_s ds, W_0(1))$, we can once again cite [26]. This density is denoted by $q$ and verifies for any bounded function $\Psi$:

$$ \mathbb{E} \left[ \Psi \left( \int_{0}^{1} r_{s} ds, W_{0}(1) \right) \right] = \int \Psi(I, w_{0}) q(I, w_{0}) dI dw_{0}. $$

Below, we denote by $q_1$ the marginal distribution of $\int_0^1 r_s ds$, whose density is bounded on $\mathbb{R}$ (see among others e.g. [39] and [26]). To describe the other components, we introduce:

$$ c_k = \mu_k - \frac{\{\sigma_k\}^2}{2}, $$
and using a conditional expectation argument, we have:

$$
\begin{align*}
\mathbb{E}[f(Y_1, S_1, \dots, S_{m'})] &= \mathbb{E} \left[ f \left( Y_1(0)e^{\mathcal{I}}, S_1(0)e^{c_1 + \sigma_1 [(\tilde{L}\tilde{W})_1 + \ell_{0,1} W_0]}, \dots, S_{m'}(0)e^{c_{m'} + \sigma_{m'} [(\tilde{L}\tilde{W})_{m'} + \ell_{0,m'} W_0]} \right) \right] \\
&= \int_{\mathbb{R}_+ \times \mathbb{R}} \int_{\mathbb{R}^{m-1}} q(I, w_0) \gamma_{m'}(\tilde{w}) \\
&\qquad \times f \left( Y_1(0)e^{\mathcal{I}}, S_1(0)e^{c_1 + \sigma_1 [(\tilde{L}\tilde{w})_1 + \ell_{0,1} w_0]}, S_{m'}(0)e^{c_{m'} + \sigma_{m'} [(\tilde{L}\tilde{w})_{m'} + \ell_{0,m'} w_0]} \right) d\tilde{w} dI dw_0,
\end{align*}
$$

where $\gamma_{m-1}$ refers to the $m-1$ dimensional standard centered Gaussian distribution. For a fixed value of $w_0$, we consider the change of variables:

$$z_1 = Y_1(0)e^I \quad \text{and} \quad z_{i+1} = S_i(0)e^{c_i + \sigma_i(\tilde{L}\tilde{w})_i + \sigma_i \ell_{0,i} w_0},$$

which is equivalent to:

$$I(z_1) = \log \left( \frac{z_1}{Y_1(0)} \right) \quad \text{and} \quad \tilde{w}(z_2, \dots, z_m, w_0) = \tilde{L}^{-1} \begin{pmatrix} \{\sigma_1\}^{-1} \left[ \log \left( \frac{z_2}{S_1(0)} \right) - c_1 - \ell_{0,1} w_0 \right] \\ \vdots \\ \{\sigma_{m'}\}^{-1} \left[ \log \left( \frac{z_m}{S_{m'}(0)} \right) - c_{m'} - \ell_{0,m'} w_0 \right] \end{pmatrix}.$$

Since $\tilde{L}$ is invertible, the Jacobian of the change of variable is then given by:

$$dId\tilde{w} = \frac{1}{(\sigma_1 \dots \sigma_{m'}) \det \tilde{L}} \frac{dz_1 dz_2 \dots dz_m}{z_1 \dots z_m}.$$

We observe that $\det L = \det \tilde{L}$, which leads to:

$$\mathbb{E}[f(Y_1, S_1, \dots, S_{m'})] = \int_{\{\mathbb{R}_+\}^m} f(z_1, \dots, z_m) \frac{\int_{\mathbb{R}} q\left(\log\left(\frac{z_1}{Y_1(0)}\right), w_0\right) \gamma_{m-1}(\tilde{w}(z_2, \dots, z_m, w_0)) dw_0}{(\sigma_1 \dots \sigma_{m'}) \det L \times (z_1 \dots z_m)} dz_1 \dots dz_m.$$

The previous expression then guarantees that $(Y_1, S_1, \dots, S_{m'})$ has a distribution uniformly continuous with respect to the Lebesgue measure on $\{\mathbb{R}_+\}^m$, whose density is of the form:

$$(z_1, \dots, z_m) \mapsto \frac{M(z_1, \dots, z_m)}{(\sigma_1 \dots \sigma_{m'}) \det(L) z_1 \dots z_m},$$

where

$$M(z_1, \dots, z_m) = \int_{\mathbb{R}} q\left(\log\left(\frac{z_1}{Y_1(0)}\right), w_0\right) \gamma_{m-1}(\tilde{w}(z_2, \dots, z_m, w_0)) dw_0.$$

The Gaussian density $\gamma_{m-1}$ being continuous and bounded by $(2\pi)^{-m'/2}$, we then have:

$$M(z_1, \dots, z_m) \le \int_{\mathbb{R}} q\left(\log\left(\frac{z_1}{Y_1(0)}\right), w_0\right) dw_0 = q_1(I(z_1)) \le \|q_1\|_{\infty} < +\infty.$$

Therefore, $M$ is a bounded continuous function, whose bound is independent from $m$ and $\Sigma$. $\square$

Below, we will need to upper bound the probability of sliced events that are defined as follows. Consider $u = (u_1, v)$ a vector of the simplex $S_m$ and $\epsilon > 0$, we introduce:

$$\Omega_{\epsilon}(\theta, u) = \{\theta \le \langle Z, u \rangle \le \theta + \epsilon\}.$$

From Proposition 15, we deduce the following property.
**Proposition 16** (Sliced events). Consider a portfolio Z defined by Equation (16). Assume $\sigma^+ = \sup_{1\le i\le m} \sigma_i < +\infty$, and that the correlation matrix $\Sigma$ is invertible. Assume that the initialization is given by $Z_i(0) = 1$. A constant K exists such that for any $\epsilon > 0$, for any $\rho \in (0, 1)$:

$$ \mathbb{P}(\Omega_{\epsilon}(\theta, u)) \le K m e^{\frac{(\sigma^{+})^2 m^2 (1-\rho)}{4\rho^2} \epsilon^{\frac{1-\rho}{2}}}, \quad \forall \theta \in \mathbb{R}, \forall u \in S_m. $$

*Proof.* We consider $\eta > 0$ small enough (whose dependency with $\epsilon$ and $m$ will be precised later on) and we observe that:

$$ \begin{aligned} \mathbb{P}(\Omega_{\epsilon}(\theta, u)) &= \mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)}] \\ &= \mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\min Z_i \ge \eta} + \mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\exists i: Z_i \le \eta}] \\ &\le \mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\min Z_i \ge \eta}] + \sum_{i=1}^{m} \mathbb{P}[Z_i \le \eta]. \end{aligned} \quad (50) $$

**Step 1:** Let us first consider the second term in Equation (50). We then consider two separate cases.

* If $i = 1$ and considering $\eta \le Z_1(0)$, we then observe that $Z_1 \ge Z_1(0)$ since $r$ is a non-negative process, so that $\mathbb{P}[Z_1 \le \eta] = 0$.

* Assume now that $i' = i - 1 \in \{1, \dots, m'\}$, using the notations of Proposition 15, we observe that:

$$ Z_{i'} \le \eta \iff Z_{i'}(0) e^{c_{i'} + \sigma_{i'} B_{i'}} \le \eta \iff B_{i'} \le -\{\sigma_{i'}\}^{-1} \left( \log \frac{1}{\eta} + c_{i'} \right), $$

so that

$$ \begin{aligned} \mathbb{P}[Z_{i'} \le \eta] &= \mathbb{P}\left[B_{i'} \le -\{\sigma_{i'}\}^{-1}\left(\log \frac{1}{\eta} + c_{i'}\right)\right] \\ &= \mathbb{P}\left[\mathcal{N}(0, 1) \ge \{\sigma_{i'}\}^{-1}\left(\log \frac{1}{\eta} + c_{i'}\right)\right] \\ &\le \frac{\sigma_{i'} e^{-\frac{(\log \frac{1}{\eta} + c_{i'})^2}{2\{\sigma_{i'}\}^2}}}{\sqrt{2\pi}\left(\log \frac{1}{\eta} + c_{i'}\right)}, \end{aligned} $$

where the last line comes from standard estimation of the Gaussian tail. Choosing $\eta$ small enough, we then deduce that a constant K exists such that:

$$ \forall i \in \{2, \dots, m\} \quad \mathbb{P}[Z_i \le \eta] \le K e^{-\frac{(\log \frac{1}{\eta})^2}{4\{\sigma^i\}^2}}. \quad (51) $$

Defining $\sigma^+ = \sup_{1\le i\le m} \sigma_i$, a union bound then leads to:

$$ \sum_{i=1}^{m} \mathbb{P}[Z_i \le \eta] \le K m e^{-\frac{(\log \frac{1}{\eta})^2}{4\{\sigma^+\}^2}} = K m \left( \eta^{\frac{\log \frac{1}{\eta}}{4\{\sigma^+\}^2}} \right). $$

**Step 2:** We study the first term of (50) and for any vector $q = (q_1, \dots, q_m)$ and for any integer $\ell$, we denote by $q^{-\ell}$ the vector $q = (q_1, \dots, q_{\ell-1}, q_{\ell+1}, \dots, q_m)$. For any $u \in S_m$, we can find $\ell$ such that $u_{\ell} \ge \frac{1}{m}$, so that:

$$ \begin{aligned} \mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\min Z_i \ge \eta}] &= \mathbb{E}[\mathbf{1}_{\theta - (Z^{-\ell}, u^{-\ell}) \le Z_{\ell} u_{\ell} \le \theta + \epsilon - (Z^{-\ell}, u^{-\ell})} \mathbf{1}_{\min Z_i \ge \eta}] \\ &\le \mathbb{E}[\mathbf{1}_{(u_{\ell})^{-1}[\theta - (Z^{-\ell}, u^{-\ell})] \le Z_{\ell} \le (u_{\ell})^{-1}[\theta - (Z^{-\ell}, u^{-\ell})] + (u_{\ell})^{-1}\epsilon} \mathbf{1}_{\min Z_i \ge \eta}]. \end{aligned} $$
Now, we recall the joint density $p(z_1, \dots, z_m)$ defined in Proposition 15 and rewrite the previous inequality as:

$$ 
\begin{aligned}
\mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\min Z_i \ge \eta}] &\le \int_{[\eta, +\infty)^{[m]}} \mathbf{1}_{(u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] \le z_{\ell} \le (u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] + (u_{\ell})^{-1}\epsilon} p(z_1, \dots, z_m) dz_1 \dots dz_m \\
&\le \int_{[\eta, +\infty)^{[m-1]}} \left[ \int_{\eta}^{\infty} \mathbf{1}_{(u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] \le z_{\ell} \le (u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] + (u_{\ell})^{-1}\epsilon} p(z_1, \dots, z_m) dz_{\ell} \right] \prod_{i \ne \ell} dz_i.
\end{aligned}
$$

We use the Cauchy-Schwarz inequality in order to bound the inner integral:

$$ 
\begin{aligned}
\int_{\eta}^{\infty} \mathbf{1}_{(u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] \le z_{\ell} \le (u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] + (u_{\ell})^{-1}\epsilon} p(z_1, \dots, z_m) dz_{\ell} & \\
&\le \sqrt{\int_{\eta}^{\infty} \mathbf{1}_{(u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] \le z_{\ell} \le (u_{\ell})^{-1}[\theta - (z^{-\ell}, u^{-\ell})] + (u_{\ell})^{-1}\epsilon} p(z_1, \dots, z_m)^2 dz_{\ell}} \\
&\le \sqrt{m\epsilon} \sqrt{\int_{\eta}^{\infty} p(z_1, \dots, z_m)^2 dz_{\ell}},
\end{aligned}
$$

where the last line derives from $u_{\ell}^{-1}\epsilon \le m\epsilon$. The previous inequality combined with the Jensen inequality $\int |g|d\mu \le \sqrt{\int g^2 d\mu}$ and the Tonelli theorem yield:

$$ 
\begin{aligned}
\mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\min Z_i \ge \eta}] &\le \sqrt{m\epsilon} \int_{[\eta, +\infty)^{[m-1]}} \left[ \sqrt{\int_{\eta}^{\infty} p(z_1, \dots, z_m)^2 dz_{\ell}} \right] \prod_{i \ne \ell} dz_i \\
&\le \sqrt{m\epsilon} \sqrt{\int_{[\eta, +\infty)^{[m]}} p(z_1, \dots, z_m)^2 dz_1 \dots dz_m}.
\end{aligned}
$$

The expression of $p$ obtained in Proposition 15 yields:

$$ 
\begin{aligned}
\mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\min Z_i \ge \eta}] &\le \sqrt{m\epsilon} \sqrt{\int_{[\eta, +\infty)^{[m]}} \frac{M(z_1, \dots, z_m)^2}{(\sigma_1 \dots \sigma_{m'})^2 \det(L)^2(z_1 \dots z_m)^2}, dz_1 \dots dz_m} \\
&\le \sqrt{m\epsilon} \sqrt{\frac{K}{\eta^m} \int_{[\eta, +\infty)^{[m]}} \frac{M(z_1, \dots, z_m)}{(\sigma_1 \dots \sigma_{m'}) \det(L)(z_1 \dots z_m)}, dz_1 \dots dz_m},
\end{aligned}
$$

where we use the fact that we integrate over $[\eta, +\infty[$ and that the function $M$ is bounded. We finally obtain that for a positive constant $K$:

$$ \mathbb{E}[\mathbf{1}_{\Omega_{\epsilon}(\theta, u)} \mathbf{1}_{\min Z_i \ge \eta}] \le \sqrt{m\epsilon} \frac{K}{\eta^{m/2}}. \quad (52) $$

**Step 3: Final bound** We now choose $\eta$ to balance the size of Equation (51) and Equation (52) in the decomposition (50). Fix $\rho \in (0, 1)$, we choose $\eta$ such that $\eta = \epsilon^{\rho/m}$, so that:

$$ (52) \le K \sqrt{m} \sqrt{\eta^{-m}} \epsilon = K \sqrt{m} \epsilon^{\frac{1-\rho}{2}}. $$

We also verify with this choice of $\eta$ that:

$$ (51) \le K m e^{-\frac{(\log \frac{1}{\eta})^2}{4(\sigma^+)^2}} \le K m e^{-\frac{m^2}{4(\sigma^+)^2 \rho^2} (\log(\epsilon))^2}. $$
To balance both contributions, we aim at finding the smallest constant C such that:

$$-\frac{m^2}{4\{\sigma^+\}^2 \rho^2} (\log(\epsilon))^2 \leq \frac{1-\rho}{2} \log(\epsilon) + C,$$

which is equivalent to

$$\frac{m^2}{4\{\sigma^+\}^2 \rho^2} (\log(\epsilon))^2 + \frac{1-\rho}{2} \log(\epsilon) \geq -C.$$

Now we consider the left hand side as a polynomial of degree 2 in $\log(\epsilon)$ and obtain that the inequality is satisfied for any $\epsilon$ when the constant C is chosen as:

$$C = \frac{m^2 \{\sigma^+\}^2 (1 - \rho)}{4\rho^2}.$$

Finally we obtain that:

$$(51) \leq K m e^{\frac{\{\sigma^+\}^2 m^2 (1-\rho)}{4\rho^2}} \epsilon^{\frac{1-\rho}{2}},$$

which concludes the proof. □

### C.3 Proof of Proposition 7

*Proof.* Let us first rewrite the Wasserstein distance as:

$$W_1(\mathcal{L}(\hat{Y}_1^h), \mathcal{L}(Y_1)) = \sup_{f \in Lip_1} \{\mathbb{E}(f(\hat{Y}_1^h) - \mathbb{E}(f(Y_1)))\}.$$

Since the G.B.M. are exactly simulated and the test function $f$ is chosen among 1-Lipschitz functions it is sufficient to bound $\mathbb{E}(|\hat{Y}_1^h - Y_1|)$. Recall the definition of the approximation $\hat{Y}^h$ in Equation (21), we obtain

$$\begin{aligned} \mathbb{E}(|Y_1 - \hat{Y}_1^h|) &= \mathbb{E}\left[\left|Y_0\left(e^{\int_0^1 r_s ds} - e^{\frac{1}{N}\sum_{k=1}^N \hat{r}_{kh}}\right)\right|\right] \\ &= \mathbb{E}\left[\left|Y_1\left(1 - e^{\Delta_h}\right)\right|\right], \end{aligned}$$

where

$$\Delta_h = \int_0^1 r_s ds - \frac{1}{N} \sum_{k=1}^N \hat{r}_{kh}. \qquad (53)$$

Using $|1 - e^x| \leq |x|e^{|x|}$ and a three way Holder's inequality, we obtain that for $p, q$ and $r$ such that $\frac{1}{p} + \frac{1}{q} + \frac{1}{r} = 1$ then:

$$\mathbb{E}(|\hat{Y}_1^h - Y_1|) \leq \mathbb{E}(|Y_1|^r)^{1/r} \mathbb{E}(|\Delta_h|^p)^{1/p} \mathbb{E}(e^{q\Delta_h})^{1/q}.$$

Using Lemma 13 and 14, we deduce that for any choice of $r, q$ such that $\frac{1}{q} + \frac{1}{r} = 1/2$, we have:

$$\mathbb{E}(|\hat{Y}_1^h - Y_1|) \leq C_q \mathbb{E}(|Y_1|^r)^{1/r} \sqrt{h}, \qquad (54)$$

for $C_q$ large enough (independent from $h$), which concludes the proof. □
### C.4 Upper bound of the bias term

*Proof of Proposition 8.* 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. $\mathbf{S} = (S_1, \dots, S_{m'})$, we have that $Z = (Y, \mathbf{S})$ whereas $\hat{Z} = \hat{Z}^h = (\hat{Y}^h, \mathbf{S})$ where $\hat{Y}^h$ is a biased approximation of $Y_1$ defined in (21). We denote by $w = (u, v) \in \mathbb{R}_+ \times (\mathbb{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$:

$$ \mathcal{B} = \left\| \mathbb{E}[\langle Z, w \rangle \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \langle \hat{Z}, w \rangle \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta}] \right\|_2 = \left\| \mathbb{E}[\langle Z \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \hat{Z} \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta}, w \rangle] \right\|_2. $$

From the re-investment condition, we know that $\|w\|_1 = 1$ and the Cauchy-Schwarz inequality yields:

$$ \mathcal{B} \le \sqrt{m} \mathbb{E} \left[ \left\| Z \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \hat{Z} \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta} \right\|_2 \right]. $$

Since we discretize only the first component, we have that:

$$ Z = (Y_0 e^{\int_0^1 r_s ds}, S_1, \dots, S_{m'}), \text{ and } \hat{Z} = (Y_0 e^{\frac{1}{N} \sum_{k=1}^N \hat{r}_{kh}}, S_1, \dots, S_{m'}), $$

We can rewrite the difference as:

$$ \begin{aligned} \| Z \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \hat{Z} \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta} \|_2 &\le \| Y \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \hat{Y} \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta} \|_2 + \| \mathbf{S}(\mathbf{1}_{\langle Z, w \rangle \ge \theta} - \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta}) \|_2 \\ &\le \| Y (\mathbf{1}_{\langle Z, w \rangle \ge \theta} - \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta}) \|_2 + \| (Y - \hat{Y}) \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta} \|_2 \\ &\quad + \| \mathbf{S}(\mathbf{1}_{\langle Z, w \rangle \ge \theta} - \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta}) \|_2 \\ &\le |Y - \hat{Y}| + (|Y| + \|\mathbf{S}\|_2) |\mathbf{1}_{\langle Z, w \rangle \ge \theta} - \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta}|. \end{aligned} $$

Taking the expectation and applying the Cauchy-Schwarz inequality in the second term we obtain that $\mathcal{B}$ can be bounded by:

$$ \mathcal{B} \le \sqrt{m} \mathbb{E} \left( (|Y| + \|\mathbf{S}\|_2)^2 \right)^{1/2} \mathbb{E} \left( \left| \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta} \right| \right)^{1/2} + \sqrt{m} \mathbb{E} \left[ \left\| Y - \hat{Y} \right\| \right]. \quad (55) $$

We apply Proposition 7 and deduce that the second term is of the order $\sqrt{h}$. We then focus on the difference of indicator functions. Let us remark that:

$$ \langle Z, w \rangle = Y u + \langle \mathbf{S}, v \rangle \quad \text{and} \quad \langle \hat{Z}, w \rangle = \hat{Y} u + \langle \mathbf{S}, v \rangle, $$

thus

$$ \mathbb{E} \left( \left| \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta} \right| \right) = \left| \mathbb{E} \left( \mathbf{1}_{\langle \mathbf{S}, v \rangle \ge \theta - Y u} - \mathbf{1}_{\langle \mathbf{S}, v \rangle \ge \theta - \hat{Y} u} \right) \right| = \mathbb{P}(A_1) + \mathbb{P}(A_2), $$

where

$$ A_1 = \{\theta \le \langle Z, w \rangle \le \theta + u(Y - \hat{Y})\} \text{ and } A_2 = \{\theta + u(Y - \hat{Y}) \le \langle Z, w \rangle \le \theta\}. $$

These two events are handled similarly, we write:

$$ Y - \hat{Y} = Y_1(1 - e^{\Delta_h}), $$

where $\Delta_h$ is given by (53). Let us fix $\epsilon > 0$ and consider the event $\Omega_\epsilon = \{u|Y - \hat{Y}| \le \epsilon\}$, then

$$ \mathbb{P}(A_1) = \mathbb{P}(A_1 \cap \Omega_\epsilon) + \mathbb{P}(A_1 \cap \Omega_\epsilon^c), $$

and now, observe that:

$$ A_1 \cap \Omega_\epsilon \subset \{\theta \le \langle Z, w \rangle \le \theta + \epsilon\}, $$
which lead to control the probability the the portfolio Z belongs to a slice. Consider $w = (u, v)$ a vector of the simplex $\Delta_m$ and $\epsilon > 0$, we introduce:

$$ \Omega_{\epsilon}(\theta, w) = \{\theta \le \langle Z, w \rangle \le \theta + \epsilon\}. $$

An upper bound for this probability was obtained in Proposition 16. A constant $K$ exists such that for any $\epsilon > 0$, for any $\rho \in (0, 1)$:

$$ \mathbb{P}(\Omega_{\epsilon}(\theta, u)) \le K m e^{\frac{(\sigma^+)^2 m^2 (1-\rho)}{4\rho^2}} \epsilon^{\frac{1-\rho}{2}}, \quad \forall \theta \in \mathbb{R}, \forall \rho \in (0,1) \quad \forall u \in \Delta_{m'}. $$

We finally study $A_1 \cap \Omega_{\epsilon}^c$ and write that for any $\zeta > 0$, and using $|1 - e^x| \le |x|e^{|x|}$:

$$ \begin{aligned} \mathbb{P}(A_1 \cap \Omega_{\epsilon}^c) &\le \mathbb{P}[uY(1 - e^{\hat{h}}) \ge \epsilon] \\ &\le \mathbb{P}[Y(1 - e^{\Delta_h}) \ge \epsilon] \\ &\le \mathbb{P}[Y \ge \zeta] + \mathbb{P}[|1 - e^{\Delta_h}| \ge \zeta^{-1}\epsilon] \\ &\le \mathbb{P}[Y \ge \zeta] + \mathbb{P}[|\Delta_h| e^{|\Delta_h|} \ge \zeta^{-1}\epsilon] \\ &\le \frac{\mathbb{E}[Y]}{\zeta} + \zeta \frac{\mathbb{E}[|\Delta_h| e^{|\Delta_h|}]}{\epsilon} \\ &\le \frac{\mathbb{E}[Y^r]}{\zeta^r} + \zeta \frac{(\mathbb{E}[\Delta_h^2])^{1/2} (\mathbb{E}[e^{2|\Delta_h|}])^{1/2}}{\epsilon}, \end{aligned} $$

where we applied first a union bound and then the Markov and the Cauchy-Schwarz inequalities. We finally obtain from Lemma 14 and 13 that:

$$ \mathbb{P}(A_1 \cap \Omega_{\epsilon}^c) \le \frac{\mathbb{E}[Y^r]}{\zeta^r} + C \zeta \frac{\sqrt{h}}{\epsilon}. \quad (56) $$

Combining Proposition 16 and (56), we deduce that:

$$ \mathbb{E} \left( \left| \mathbf{1}_{\langle Z, w \rangle \ge \theta} - \mathbf{1}_{\langle \hat{Z}, w \rangle \ge \theta} \right| \right) \le C \left[ K m e^{\frac{(\sigma^+)^2 m^2 (1-\rho)}{4\rho^2}} \epsilon^{\frac{1-\rho}{2}} + \frac{\mathbb{E}[Y^r]}{\zeta^r} + C \zeta \frac{\sqrt{h}}{\epsilon} \right]. $$

Now using this bound as well as (54) in (55) we obtain that:

$$ \mathcal{B} \le K \sqrt{m} e^{\frac{(\sigma^+)^2 m^2 (1-\rho)}{4\rho^2}} \left[ \epsilon^{\frac{1-\rho}{2}} + \frac{1}{\zeta^r} + \zeta \frac{\sqrt{h}}{\epsilon} + \sqrt{h} \right], \quad (57) $$

where the parameters $r > 1$, $\epsilon > 0$, $\zeta > 0$ and $\rho \in (0, 1)$ have to be chosen. Note that the constant depends heavily on the number of assets $m$ and on the choice of $\rho$. We optimize only on the two parameters $\epsilon > 0$ and $\zeta > 0$ the quantity:

$$ \epsilon^{\frac{1-\rho}{2}} + \frac{1}{\zeta^r} + \zeta \frac{\sqrt{h}}{\epsilon} + \sqrt{h}. $$

It is immediate to verify that $\zeta$ should be chosen such that

$$ \zeta^{-r} = \zeta \frac{\sqrt{h}}{\epsilon} \implies \zeta = \frac{\epsilon^{\frac{1}{1+r}}}{\sqrt{h}^{\frac{1}{1+r}}}. $$

In the meantime, the optimal value of $\epsilon$ satisfies:

$$ \epsilon^{\frac{1-\rho}{2}} = \zeta \frac{\sqrt{h}}{\epsilon} \implies \epsilon^{\frac{3-\rho}{2}} = \epsilon^{\frac{1}{1+r}} \sqrt{h}^{1-\frac{1}{1+r}}. $$
Some straightforward computations show that $\epsilon = h^{\frac{2r}{1-\rho+r(3-\rho)}}$, which in turn implies that a constant $K$ large enough exists such that:

$$ \mathcal{B} \leq K_{r, \rho} \sqrt{me^{\frac{(\sigma^+)^2 m^2 (1-\rho)}{4\rho^2}}} h^{2\left(\frac{3-\rho}{1-\rho} + \frac{1}{r}\right)}. $$

It is then easy to see that the rate may be arbitrarily close to $h^{1/6}$ by picking $r$ large enough and $\rho$ close to 0. More precisely, for any $\epsilon > 0$, we can find $\rho$ and $r$ such that:

$$ \mathcal{B} \leq K_{\epsilon} \sqrt{me^{\frac{(\sigma^+)^2 m^2}{4\epsilon^2}}} h^{\frac{1}{6}-\epsilon}. $$

□

## References

[1] A. Alfonsi. On the discretization schemes for the cir (and bessel squared) processes. *Monte Carlo Methods and Applications*, 11:355-384, 2005.

[2] A. Alfonsi. High order discretization schemes for the cir process: application to affine term structure and heston models. *Mathematics of Computation*, 79(269):209-237, 2010.

[3] A. Alfonsi. Strong order one convergence of a drift implicit euler scheme: Application to the cir process. *Statistics & Probability Letters*, 83(2):602-607, 2013.

[4] A. Alfonsi. *Affine diffusions and related processes: simulation, theory and applications*, volume 6. Springer, 2015.

[5] P. Artzner, F. Delbaen, J.M. Eber, and D. Heath. Thinking coherently. *Risk*, 10, 1997.

[6] P. Artzner, F. Delbaen, J.M. Eber, and D. Heath. Coherent measures of risk. *Mathematical Finance*, 9:203-228, 1999.

[7] Y. Atchadé, G. Fort, and E. Moulines. On perturbed proximal gradient algorithms. *The Journal of Machine Learning Research*, 18(1):310-342, 2017.

[8] V. Bally and D. Talay. The law of the Euler scheme for stochastic differential equations. *Probability theory and related fields*, 104(1):43-60, 1996.

[9] O. Bardou, N. Frikha, and G. Pagès. Recursive computation of value-at-risk and conditional value-at-risk using mc and qmc. In Pierre L'Ecuyer and Art B. Owen, editors, *Monte Carlo and Quasi-Monte Carlo Methods 2008*, pages 193-208, Berlin, Heidelberg, 2009. Springer Berlin Heidelberg.

[10] O. Bardou, N. Frikha, and G. Pagès. Cvar hedging using quantization-based stochastic approximation algorithm. *Mathematical Finance*, 26(1):184-229, 2016.

[11] O. Bardou, N. Frikha, and G. Pagès. Computing var and cvar using stochastic approximation and adaptive unconstrained importance sampling. *Monte Carlo Methods and Applications*, 15(3):173-210, 2009.

[12] Amir Beck and Marc Teboulle. Mirror descent and nonlinear projected subgradient methods for convex optimization. *Operations Research Letters*, 31(3):167-175, 2003.

[13] A. Ben Tal and M. Teboulle. Expected utility, penalty functions and duality in stochastic nonlinear programming. *Manage. Sci.*, 32:1445-1466, 1986.
[14] M. Benaïm. Dynamics of stochastic approximation algorithms. *Séminaire de probabilités de Strasbourg*, 33:1–68, 1999.

[15] M. Benaïm and M. Hirsch. Asymptotic pseudotrajectories and chain recurrent flows, with applications. *Journal of Dynamics and Differential Equations*, 8(1):141–176, 1996.

[16] S. Bubeck. Convex optimization: Algorithms and complexity. *Foundations and Trends in Machine Learning*, (8), 2015.

[17] L. Chen and Q.-M. Shao. Stein's method for normal approximation. In *An introduction to Stein's method*, volume 4 of *Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap.*, pages 1–59. Singapore Univ. Press, Singapore, 2005.

[18] L. Condat. Fast projection onto the simplex and the 11 ball. *Mathematical Programming Series A*, 158:575–585, 2016.

[19] J.C. Cox, J.E. Ingersoll, and S.A. Ross. A theory of the term structure of interest rates. *Econometrica*, 53:385–407, 1985.

[20] A. Cozma and C. Reisinger. Exponential integrability properties of Euler discretization schemes for the cox-ingersoll-ross process. *Discrete and continuous dynamical systems, series B*, 21(10):3359–3377, 2016.

[21] S. Dereich, A. Neuenkirch, and L. Szpruch. An Euler-type method for the strong approximation of the cox ingersoll ross process. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 468(2140):1105–1115, 2012.

[22] J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra. Efficient projections onto the 1 1 -ball for learning in high dimensions. *Int. Conf. on Machine learning (ICML)*, 2008.

[23] N. Frikha and L. Huang. A multi-step richardson-romberg extrapolation method for stochastic approximation. *Stochastic Processes and their Applications*, 125(11):4066–4101, 2015.

[24] P. Glasserman. *Monte Carlo Methods in Financial Engineering*. Stochastic Modelling and Applied Probability. Springer, 2003.

[25] A. Going-Jaeschke. Parameter estimation and bessel processes in financial models and numerical analysis in hamiltonian dynamics. *PhD dissertation of Swiss Federal Institutde of Technology of Zurich*, 1998.

[26] A. Gulisashvili and E. Stein. Asymptotic behavior of the stock price distribution density and implied volatility in stochastic volatility models. *Applied Mathematics and Optimization*, 61(3):287–315, 2010.

[27] K.C. Kiwiel. Breakpoint searching algorithms for the continuous quadratic knapsack problem. *Math. Program., Ser. A*, 112:473–491, 2008.

[28] P. Kloeden and E. Platen. *Numerical solution of SDE*. Springer Science, 1999.

[29] P. Krokhmal, J. Palmquist, and S. Uryasev. Portfolio optimization with conditional value-at-risk objective and constraints. *Journal of Risk*, 4(2):43–68, 2001.

[30] G. Lan, A. Nemirovskij, and A. Shapiro. Validation analysis of mirror descent stochastic approximation method. *Math. Program., Ser. A*, (134):425–458, 2012.

[31] Sophie Laruelle and Gilles Pagès. Stochastic approximation with averaging innovation applied to finance. *Monte Carlo Methods and Applications*, 18(1):1–51, 2012.
[32] A. Fujino M. Blondel and N. Ueda. Large-scale multiclass support vector machine training via euclidean projection onto the simplex. *International Conference on Pattern Recognition*, 2014.

[33] H.M. Markowitz. Portfolio selection. *Journal of Finance*, 7:77-91, 1952.

[34] P. Mertikopoulos and M. Staudigl. On the convergence of gradient-like flows with noisy gradient input. *SIAM Journal on Optimization*, 28(1):163-197, 2018.

[35] A. Nemirovskij and D. Yudin. *Problem complexity and method efficiency in optimization*. Wiley-Interscience, 1983.

[36] Y. Nesterov. Primal-dual subgradient methods for convex problems. *Mathematical Programming, Series B*, 120:221-259, 2009.

[37] [DATA] Python package `panda_datareader`. `get_data_yahoo("nvda", "tsla", "aapl", "csco", "t", "xom", "twtr", "jpm", "bp")`. `start="2016-01-01"`, `end="2018-01-01"`.

[38] G.Ch. Pflug. Some remarks on the value-at-risk and the conditional value-at-risk. In *Probabilistic Constrained Optimization: Methodology and Applications, Ed. S. Uryasev*. Kluwer Academic Publishers, 2000.

[39] M. Pitman, J. Yor. A decomposition of bessel bridges. *Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete*, 59(4):425-457, 1982.

[40] H. Robbins and S. Monro. A stochastic approximation method. *The Annals of Mathematical Statistics*, 22:400-407, 1951.

[41] H. Robbins and D. Siegmund. A convergence theorem for non negative almost supermartingales and some applications. *Optimizing methods in stat.*, pages 233-257, 1971.

[42] R. Rockafellar and J. Royset. Random variables, monotone relations, and convex analysis. *Mathematical Programming*, 148(1):297-331, 2014.

[43] R. T. Rockafellar and S. Uryasev. Optimization of conditional value-at-risk. *The Journal of Risk*, 2(3):21-41, 2000.

[44] R. T. Rockafellar and S. Uryasev. Conditional value-at-risk for general loss distributions. *Journal of Banking and Finance*, 26(7):1443-1471, 2002.

[45] W. F. Sharpe. Mutual fund performance. *Journal of Business*, 39:119-138, 1966.

[46] D. Talay and L. Tubaro. Expansion of the global error for numerical schemes solving stochastic differential equations. *Stochastic analysis and applications*, 8(4):483-509, 1990.

[47] Z. Zhou, P. Mertikopoulos, N. Bambos, S. Boyd, and P. Glynn. Stochastic mirror descent in variationally coherent optimization problems. *Advances in Neural Information Processing Systems*, 30, 2017.