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

HAL
open science

# Stochastic approximation, Markovian perturbation and their applications

Noufel Frikha

► To cite this version:

Noufel Frikha. Stochastic approximation, Markovian perturbation and their applications: Stochastic approximation, Markovian perturbation and their applications. Probability [math.PR]. Université Paris Cité, 2017. tel-04017754

HAL Id: tel-04017754

https://hal.science/tel-04017754

Submitted on 7 Mar 2023

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.

Public Domain
# Habilitation à Diriger des Recherches

Spécialité : Mathématiques Appliquées

présentée
par

Noufel Frikha

## Stochastic approximation, Markovian perturbation of stochastic processes and their applications.

soutenue publiquement le 23 Novembre 2017 devant le jury composé de :

<table><tr><td>Mme</td><td>Annie Millet,</td><td>Professeur à l'université Paris 1,</td><td>Examineur</td></tr><tr><td>M.</td><td>Vlad Bally,</td><td>Professeur à l'université Paris Est,</td><td>Examineur</td></tr><tr><td>M.</td><td>Dan Crisan,</td><td>Professeur à l'Imperial College,</td><td>Rapporteur</td></tr><tr><td>M.</td><td>François Delarue,</td><td>Professeur à l'université Nice-Sophia Antipolis,</td><td>Examineur</td></tr><tr><td>M.</td><td>Gilles Pagès,</td><td>Professeur à l'université Pierre et Marie Curie,</td><td>Examineur</td></tr><tr><td>M.</td><td>Denis Talay,</td><td>Directeur de recherche à l'INRIA,</td><td>Rapporteur</td></tr><tr><td>M.</td><td>Peter Tankov,</td><td>Professeur à l'ENSAE,</td><td>Rapporteur interne</td></tr></table>
Maturité de l'homme:
cela veut dire avoir retrouvé le sérieux qu'enfant, on mettait dans ses jeux.

Friedrich Nietzsche, par delà le bien et le mal, maximes et intermèdes.
Remerciements

En tout premier lieu, je tiens à remercier içi les rapporteurs de mon mémoire d'HDR, Dan Crisan, Denis Talay et Peter
Tankov qui me font le grand honneur d'avoir accepté d'évaluer mes travaux. Je vous suis particulièrement reconnaissant
pour vos remarques et pour l'intérêt dont vous avez témoigné pour mon travail. Je remercie chaleureusement les autres
membres du jury, Annie Millet, Vlad Bally, François Delarue et Gilles Pagès. Votre présence aujourd'hui ainsi que
votre bienveillance appellent ma très grande reconnaissance.

Je remercie particulièrement Annie Millet. La grande clarté de son cours de Master “Méthodes de Monte Carlo”,
son écoute ainsi que sa disponibilité sont pour moi des exemples. Elle a su donner à chacun de nos entretiens une note
de gentillesse qui m’a beaucoup touchée. Je remercie Peter Tankov pour sa constante disponibilité, l’intérêt dont il a
fait preuve à l’égard de l’ensemble de mes travaux et ses précieux conseils. Je remercie également Vlad Bally dont les
travaux sur le calcul de Malliavin appliqué aux équations différentielles stochastiques et l’interprétation probabiliste
de la méthode parametrix m’ont beaucoup apportés ces dernières années. Je tiens à remercier chaleureusement
François Delarue dont la qualité et la profondeur des travaux, notamment ceux sur les jeux à champs moyens, éveillent
constamment ma curiosité. Enfin, en dernier lieu, je souhaite exprimer içi toute ma gratitude à Gilles Pagès, directeur
de thèse exemplaire, pour son enthousiasme et humour communicatifs. Je ferai de mon mieux pour suivre vos traces
avec mes futurs élèves en thèse!

Je tiens également à remercier içi l'ensemble de mes collaborateurs mathématiciens, Olivier Bardou, Vincent Lemaire, Abass Sagna, Stéphane Menozzi, Max Fathi, Lorick Huang, Arturo Kohatsu-Higa et Libo Li. Vous m'avez beaucoup apporté, chacun à votre propre manière, et ce travail vous doit beaucoup. Je suis sûr que nous aurons l'occasion de nous retrouver et de retravailler ensemble!

Je tiens à exprimer ma sincère reconnaissance à l'ensemble de mes collègues du LPMA, et particulièrement l'équipe
MathFiProNum. Je vous remercie pour votre bienveillance et sympathie quotidiennes. C'est aussi grâce à vous qu'il
m'est si agréable de venir travailler jour après jour.

Pour terminer, je souhaite exprimer toute ma gratitude à ma (petite) famille et mes amis ainsi qu'à tous ceux qui
ont su me soutenir ces dernières années.
# Outline of the manuscript

This document presents and put in perspective my research work carried out after my PhD Thesis since my position as an associate professor at the university Paris Diderot (Paris VII) started in September 2012.

This work consists in two topics that interact with one another. The first part concerns stochastic approximation algorithm with some applications in mathematical finance. The second part focuses on Markovian perturbation of stochastic processes with several applications such that probabilistic representation, unbiased Monte Carlo simulation, density estimates, asymptotic expansions of some Markov semigroups, weak approximation error of some discretisation schemes.

The first four chapters deals with stochastic approximation theory. The first chapter is an introduction. We briefly present two problems arising in numerical probability. The first one deals with stochastic approximation algorithm, also referred in the literature as Robbins-Monro algorithms, which are simulation based recursive schemes that are widely used in stochastic optimization. We review results about a.s. convergence, strong and weak convergence rates and present some applications in mathematical finance. The second one concerns several problems related to the discretisation of stochastic processes. We present the main results available in the literature and the different tools to obtain them, with a special emphasis on the weak discretisation error.

Chapter 2 deals with some contributions to risk management in financial and energy markets. It summarises the results contained in [5] and [6]. First, we study theoretical and computational aspects of risk minimisation in financial market models operating in discrete time [6]. To define the risk, we consider a class of convex risk measures defined on $L^p$ in terms of shortfall risk. A dynamic programming principle in a non-Markovian framework is established, under mild assumptions, namely the absence of arbitrage opportunity and the non-degeneracy of the price process. Moreover, we provide a probabilistic algorithm to compute the optimal dynamic strategy using optimal quantization and analyse the global error between the true value function and the approximated one.

As an alternative method to simulation based recursive importance sampling procedure to estimate the optimal change of measure for Monte Carlo simulations, we propose an algorithm which combines (vector and functional) optimal quantization with Newton-Raphson zero search procedure [5]. We analyse the error induced by the proposed procedure and illustrate its efficiency on several examples.

In Chapter 3, we focus on non-asymptotic concentration bounds, deviation estimates and related transport entropy inequalities for stochastic approximation schemes [7], [8]. Our results include the deviation between the expectation of a given function of an Euler like discretization scheme of some diffusion process at a fixed deterministic time and its empirical mean obtained by the Monte-Carlo procedure as well as the deviation between the value at a given time of a stochastic approximation algorithm and its target. Under suitable assumptions both concentration bounds turn out to be Gaussian [7]. In order to obtain different non-Gaussian concentration regimes, one has to quantify the contribution of the diffusion term to the concentration regime [8]. We also derive a general non-asymptotic deviation bound for the difference between a function of the trajectory of a continuous Euler scheme associated to a diffusion process and its mean. Finally, we obtain non-asymptotic bound for stochastic approximation with averaging of trajectories. Notably, it is shown that averaging a stochastic approximation algorithm with a slow decreasing step sequence provide the optimal concentration rate.

In Chapter 4, we present a class of multi-level and multi-step Richardson Romberg extrapolation methods for stochastic approximation algorithms [9], [10]. In a first part, we extend the scope of the multilevel Monte Carlo method recently introduced for the expectation of a given function of some diffusion process to the framework of stochastic optimisation by means of stochastic approximation algorithm. A central limit theorem is proved and the possible optimal choices of step size sequence are provided [9]. In a second part, we also extend the well-known Richardson-Romberg extrapolation method for Monte Carlo linear estimator to the framework of stochastic approximation algorithm [10]. We illustrate the method to the estimation of the quantile of diffusion processes. For both methods, numerical results confirm the theoretical analysis and show a significant reduction in the initial computational cost.

The next chapters are related to the second part of my research work which focuses on Markovian perturbation of stochastic processes and its application.

Chapter 5 is dedicated to two approximation problems [12], [11]. We begin by a brief presentation of the parametrix
method “à la Mc Kean and Singer” in the framework of diffusion processes. This perturbation method is a classical
technique in PDE theory to construct the fundamental solution to an elliptic or parabolic PDE. Our presentation here
is based on probabilistic arguments. We also briefly present its discrete counterpart when one considers discretisation
schemes of stochastic processes.

In the second part of this chapter, we consider a discrete version of the parametrix technique that we apply in order
to study the weak approximation error of a one-dimensional skew diffusion with bounded measurable drift and Hölder
diffusion coefficient by an Euler-type scheme, which consists of iteratively simulating skew Brownian motions with
constant drift. First, we establish two sided Gaussian bounds for the density of this approximation scheme. Then, a
bound for the difference between the densities of the skew diffusion and its Euler approximation is obtained under the
same assumption [12].

In the third part, we study the asymptotic expansion of a Markov semigroup with respect to a small parameter
[11]. The method proposed combines a perturbation approach of Markov semigroup and Malliavin's calculus. We
study two applications. In the first one, heat kernel expansions of some hypo-elliptic diffusion processes are revisited.
In the second one, we obtain an asymptotic expansion of the transition density of a skew diffusion process with small
skew.

Chapter 6 is divided into two parts. In the first part, we focus on the regularity of the law associated to the first hitting time of a threshold by a one-dimensional uniformly elliptic diffusion process and to the associated process stopped at the threshold [13]. We present explicit expressions for the corresponding transition densities and studies their regularity properties up to the boundary under mild assumptions on the coefficients. Some Gaussian upper-estimates for these densities (and their derivatives) are also established. Our construction allows to obtain a probabilistic representation that can be used for the construction of an unbiased Monte Carlo simulation method and also to derive some integration by parts formula.

In the second part, we investigate the weak existence and uniqueness for some SDEs with coefficients depending on
some of its path-functionals. A general method is developed in order to deal with a process whose law is singular with
respect to the Lebesgue measure. We illustrate our approach on two examples: a diffusion process with coefficients
depending on its running local time and a diffusion process with coefficients depending on its running maximum. We
also construct the associated transition density and establish some Gaussian upper-estimates.
List of research works presented in this manuscript:

Publications related to my PhD Thesis:

[1] with O. Bardou, G. Pagès: Recursive Computation of Value-at-Risk and Conditional Value-at-Risk using MC and QMC, *Monte Carlo and Quasi-Monte Carlo Methods 2008*, Part 3, 193-208, 2008.

[2] with O. Bardou, 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.

[3] with O. Bardou, G. Pagès: CVaR hedging using quantization based stochastic approximation algorithm, *Mathematical Finance*, 26(1):184-229, 2016.

[4] with V. Lemaire: Joint modelling of Gas and Electricity spot prices, *Applied Mathematical Finance*, 20(1), 69-93, 2013.

Publications after my PhD Thesis:

[5] with A. Sagna: Quantization based Recursive Importance Sampling, *Monte Carlo Methods and Applications*, 18(4):287-326, 2012.

[6] Shortfall risk minimization in discrete time financial market models, *SIAM Journal on Financial Mathematics*, Volume 5, Issue 1, 384-414, 2012.

[7] with S. Menozzi: Concentration Bounds for Stochastic Approximations, *Electronic Communications in Probability*, Volume 17, no. 47, 1-15, 2012.

[8] with M. Fathi: Transport-Entropy inequalities and deviation estimates for stochastic approximation schemes, *Electronic Journal of Probability*, Volume 18, no. 67, 1-36, 2013.

[9] Multi-level stochastic approximation algorithms, *The Annals of Applied Probability*, Volume 26, Issue 2, 933-985, 2016.

[10] with L. Huang: A multi-step richardson-romberg extrapolation method for stochastic approximation, *Stochastic Processes and their Applications*, Volume 125, Issue 11, 4066-4101, 2015.

[11] with A. Kohatsu-Higa: A Parametrix approach for asymptotic expansion of Markov semigroups with applications to multi-dimensional diffusion processes, 50 pages, submitted, 2015.

[12] On the weak approximation of a skew diffusion by an Euler-type scheme, accepted for publication in *Bernoulli*, 28 pages, 2016.

[13] with A. Kohatsu-Higa, L. Li: On the first hitting times of one dimensional elliptic diffusions, 41 pages, submitted, 2016.

[14] with L. Li: Weak existence and uniqueness for some stochastic differential equations with path-functional coefficients, 26 pages, submitted, 2017.
# Contents

<table><tr><td>Remerciements</td><td>3</td></tr><tr><td>Outline of the manuscript</td><td>5</td></tr><tr><td>Contents</td><td>9</td></tr><tr><td>1 Introduction : On two problems arising in numerical probability</td><td>11</td></tr><tr><td>1.1 Stochastic approximation algorithm</td><td>11</td></tr><tr><td>1.1.1 Introduction</td><td>11</td></tr><tr><td>1.1.2 Standard a.s. convergence results</td><td>13</td></tr><tr><td>1.1.3 Convergence rate</td><td>15</td></tr><tr><td>1.2 Discretisation and simulation of stochastic processes</td><td>17</td></tr><tr><td>1.2.1 Motivation</td><td>17</td></tr><tr><td>1.2.2 Strong error</td><td>18</td></tr><tr><td>1.2.3 Weak error</td><td>19</td></tr><tr><td>1.2.4 Multi-level Monte Carlo simulation</td><td>23</td></tr><tr><td>1.2.5 Unbiased simulation of SDE</td><td>24</td></tr><tr><td>Part I : Stochastic approximation algorithm and applications</td><td>27</td></tr><tr><td>2 Some contributions to risk management in financial and energy markets</td><td>29</td></tr><tr><td>2.1 Risk minimisation in financial or energy markets operating in discrete time</td><td>29</td></tr><tr><td>2.2 Quantisation based recursive importance sampling</td><td>32</td></tr><tr><td>3 Concentration inequalities and deviation estimates for stochastic approximation schemes</td><td>39</td></tr><tr><td>3.1 A brief presentation of transport inequalities and concentration inequalities</td><td>39</td></tr><tr><td>3.2 Gaussian concentration inequalities</td><td>41</td></tr><tr><td>3.3 From Exponential to Gaussian concentration regimes</td><td>44</td></tr><tr><td>3.4 Perspective</td><td>46</td></tr><tr><td>4 Improving the computational efficiency of stochastic approximation algorithms</td><td>49</td></tr><tr><td>4.1 Motivation</td><td>49</td></tr><tr><td>4.2 Multi-level stochastic approximation algorithm</td><td>51</td></tr><tr><td>4.3 Multi-step Richardson-Romberg for stochastic approximation algorithm</td><td>53</td></tr><tr><td>4.4 Perspective</td><td>58</td></tr><tr><td>Part II : Markovian perturbation of stochastic processes and applications</td><td>59</td></tr><tr><td>5 Weak approximation and asymptotic expansion of some Markov semigroups</td><td>61</td></tr><tr><td>5.1 A brief presentation of the parametrix method “à la Mc Kean and Singer”</td><td>61</td></tr><tr><td>5.1.1 The case of diffusion processes</td><td>61</td></tr><tr><td>5.1.2 A discrete version of the parametrix method</td><td>63</td></tr><tr><td>5.2 Weak approximation of a skew diffusion by an Euler-type scheme</td><td>64</td></tr><tr><td>5.3 Asymptotic expansion of Markov semigroups and application to diffusion processes</td><td>66</td></tr><tr><td>6 One-dimensional elliptic diffusions and some of their path-functionals</td><td>71</td></tr><tr><td>6.1 Hitting time and killed process associated to a one-dimensional diffusion</td><td>71</td></tr><tr><td>6.2 Weak uniqueness and density estimates for SDEs with coefficients depending on some path-functionals</td><td>74</td></tr><tr><td>6.3 Perspective</td><td>79</td></tr><tr><td>Bibliography</td><td>81</td></tr></table>
# CONTENTS
# Chapter 1

# Introduction : On two problems arising in numerical probability

In this introductory chapter, we briefly present two problems arising in numerical probability. The first part deals with a brief presentation of the stochastic approximation algorithm theory. Since their introduction in the seminal paper of H. Robbins and S. Monro [117], these simulation based procedures are widely used in various fields such that stochastic optimisation or zero search problems. The stochastic approximation theory provides various theorems that guarantee the a.s., weak and/or strong convergence rates towards the desired target. We present here some results related to the a.s. convergence, known as the Robbins-Siegmund theorem and the ordinary differential equation (ODE) method. We provide the key points for their proofs. Then, we review strong and weak convergence rates. Such algorithm satisfies a central limit theorem (CLT) with a rate that depends on the step sequence. Finally, we present some applications in mathematical finance.

In the second part, we briefly present several problems related to the discretisation of diffusion-jump type stochastic processes. We distinguish between strong and weak approximation since these two problems are of different nature and require different mathematical tools for their study. We first investigate the strong approximation error which is related to the error between the paths of the process and of its approximation scheme. Then, we focus on the weak approximation error which is the error between the law of the process and of its approximation scheme. This error is often estimated on a class of test functions. Since the pioneering paper of Talay and Tubaro [122], we know that there is a strong connection between the weak approximation error and the solution of the underlying parabolic linear partial differential equation (PDE) via the Feynman-Kac representation.

We try to gather major results obtained the past thirty years since they will be useful in order to understand some of the next chapters.

We will focus on the weak approximation and the mathematical tools needed in order to study the resulting error. More precisely, we will distinguish between the two following cases: when the test function and the coefficients of the stochastic differential equation are smooth; when the process is non-degenerated. This will be the topic of Section 5.2 in Chapter 5 which deals with the weak approximation error of a skew diffusion by an Euler-type scheme.

We briefly present the Multi-level Monte Carlo method developed by Giles [48] which allows to significantly increase the computational efficiency of the Monte Carlo method for the computation of the expectation of a non-simulatable random variable. We extend the scope of this method to the framework of stochastic optimisation by means of stochastic approximation algorithm in Chapter 4.

We conclude this chapter by presenting an unbiased simulation technique for stochastic differential equations based on a probabilistic representation of the parametrix technique that has been recently developed Bally and Kohatsu-Higa [7]. In Chapter 6, we will propose a similar probabilistic representation for hitting times and associated killed process of one-dimensional elliptic diffusions.

## 1.1 Stochastic approximation algorithm

### 1.1.1 Introduction

In applied mathematics, one often has to numerically solve some optimisation or zero search problems. The former problem can often be cast into the latter one when one faces a convex optimisation problem. Some commonly encountered examples are the computation of the quantile at level $\alpha \in (0, 1)$ of probability distribution, in mathematical finance the extraction of implicit volatility or the minimisation of a convex function with respect to a parameter (e.g. quantity of assets in a portfolio, minimisation of the variance, etc).
In a deterministic framework, the zero-search recursive procedure has the following dynamics

$$ \theta_{n+1} = \theta_n - \gamma_{n+1}h(\theta_n), \quad n \ge 0 \qquad (1.1) $$

where $h: \mathbb{R}^d \to \mathbb{R}^d$ is a continuous function satisfying a sub-linear growth assumption and $(\gamma_n)_{n \ge 1}$ is a non-negative non-increasing step sequence such that $\sum_{n \ge 1} \gamma_n = +\infty$.

Under an appropriate mean-reverting assumption on the function $h$, one shows that the deterministic sequence $(\theta_n)_{n \ge 1}$ defined by (1.1) is bounded and eventually converges toward a zero $\theta^*$ of $h$. More precisely, assuming that the following mean-reverting assumption

$$ \langle \theta - \theta^*, h(\theta) \rangle \ge 0, \quad \forall \theta \ne \theta^* \qquad (1.2) $$

and that $\{h=0\} = \{\theta^*\}$, then one proves that $(\theta_n)_{n \ge 1}$ converges towards $\theta^*$ as $n \to \infty$. An example of such procedure is the Newton-Raphson zero-search algorithm obtained by choosing $h := Dh^{-1}h$ and $\gamma_n = 1$.

However, in a stochastic optimization framework as described in the examples mentioned above, one does not have straightforward access to the value of $h(\theta)$ because $h$ is given by the expectation of some given random function, say

$$ h(\theta) := \mathbb{E}[H(\theta, U)] $$

where $H: \mathbb{R}^d \times \mathbb{R}^q \to \mathbb{R}^d$ is a measurable function and $U$ is an $\mathbb{R}^q$-valued random variable with law $\mu \in \mathcal{P}(\mathbb{R}^q)$. Here it is implicitly assumed that one can easily simulate according to the law $\mu$ and has direct access to the value of the function $H$. One simple idea could be to simply approximate the above dynamics (1.1) by using at each iterate a Monte Carlo approximation of $h(\theta_n)$. A more interesting idea consists in doing both simultaneously by using on the one hand $H(\theta_n, U^{(n+1)})$ instead of $h(\theta_n)$ and on the other hand by letting the step sequence goes to zero not too fast to take into account the additional variance induced by this randomisation. Hence, one is led to introduce the following dynamics

$$ \theta_{n+1} = \theta_n - \gamma_{n+1} H(\theta_n, U^{(n+1)}), \quad n \ge 0 \qquad (1.3) $$

where $(U^{(n)})_{n \ge 1}$ is an i.i.d. sequence of random vector distributed according to $\mu$ on a given probability space $(\Omega, \mathcal{F}, \mathbb{P})$ and $\theta_0$ is an $\mathbb{R}^d$-valued random vector independent of the sequence $(U^{(n)})_{n \ge 1}$.

### Example:

As a first simple example, consider the empirical average $(\bar{U}^n)_{n \ge 1}$ of i.i.d. and integrable samples $(U^{(n)})_{n \ge 1}$ defined by $\bar{U}^n := \frac{1}{n} \sum_{k=1}^n U^{(k)}$, $\bar{U}^0 = 0$. We observe that it satisfies the following dynamics

$$ \bar{U}^{n+1} = \bar{U}^n - \frac{1}{n+1} (\bar{U}^n - U^{(n+1)}) = \bar{U}^n - \frac{1}{n+1} H(\bar{U}^n, U^{(n+1)}), \quad n \ge 0, \qquad (1.4) $$

where $H(\theta, U) := \theta - U$ and $h(\theta) := \theta - \mathbb{E}[U]$, so that $\theta^* = \mathbb{E}[U]$. Moreover, it is well-known that the sequence $(\bar{U}^n)_{n \ge 1}$ given by (1.4) converges a.s. to the unique zero of $h$. When one knows how to sample from the distribution $\mu$, this procedure is easy to implement contrary to its deterministic counterpart (1.1). We see that stochastic approximation appears as a natural extension of the Monte Carlo method to solve optimisation or zero search problems when the function of interest $h$ writes as the integral of some function against a law $\mu$ from which it is easy to sample.

Let us make an important remark. Writing $H(\theta_n, U^{(n+1)}) = h(\theta_n) + \Delta M_{n+1}$ where $\Delta M_{n+1} := H(\theta_n, U^{(n+1)}) - h(\theta_n)$ is a martingale increment with respect to the natural filtration $\{\mathcal{F}_n\}$, $\mathcal{F}_n := \sigma(\theta_0, U^{(1)}, \dots, U^{(n)})$, $n \ge 1$ induced by the algorithm, we observe that the dynamics (1.3) can be seen as a perturbed Euler scheme (with decreasing step) of the ordinary differential equation (ODE)

$$ \dot{\theta} = -h(\theta) \qquad (1.5) $$

From this connection it may be motivating to use some ODE techniques such that existence of a Lyapunov function in order to study the asymptotic behaviour of the recursive scheme (1.3). A Lyapunov function for the ODE (1.5) is a function $L: \mathbb{R}^d \to \mathbb{R}_+$ such that any solution $t \mapsto \theta(t)$ satisfies $t \mapsto L(\theta(t))$ is non-increasing. Assuming that $L$ is continuously differentiable, it is mainly equivalent to the condition $\langle \nabla L, h \rangle \ge 0$ since $dL(\theta(t)) = \langle \nabla L(\theta(t)), \theta(t) \rangle = -\langle \nabla L, h \rangle(\theta(t))dt$. Finding a Lyapunov function may be a difficult task and there is no result that guarantees its existence! Usually it requires a deep knowledge and understanding of the underlying dynamic (1.5). However, if such function does exist, the system (1.5) is said to be dissipative. We will have in mind two commonly encountered situations:

* when one faces a convex optimisation problem, the recursive scheme (1.3) is usually a stochastic gradient algorithm. The function $L$ is the object of interest that is one has $h = \nabla L$ and one is computing a minimum of $L$.

* Unfortunately, as is so often the case when dealing with stochastic approximation algorithm, the function $L$ is not identified a priori and one has to find a Lyapunov function associated to the function $h$, if such function exists!
1.1.2 Standard a.s. convergence results

We now turn our attention to results related to the a.s. and/or Lp convergence of the recursive procedure (1.3). We first provide a general result known as Robbins-Siegmund Lemma based on the existence of an a priori Lyapunov function from which we will deduce two corollaries that ensures the a.s. (and Lp) convergence of the sequence (θn)n≥1. We refer the interested reader to the monographs of Benveniste and al. [19], Duflo [39] and Kushner and Yin [81] among others for developments and a more complete overview in stochastic approximation theory.

**Theorem 1.1.1** (Robbins-Siegmund Lemma). Assume that there exists a continuously differentiable function L : ℝd → ℝ+ such that ∇L is Lipschitz-continuous and

$$
|\nabla L|^2 \le C(1 + L) \tag{1.6}
$$

and such that the mean-reverting assumption holds

$$
\langle \nabla L, h \rangle \ge 0. \tag{1.7}
$$

Assume furthermore that H satisfies the following sub-linear growth assumption

$$
\forall \theta \in \mathbb{R}^d, \quad \mathbb{E}[|H(\theta, U)|^2]^{\frac{1}{2}} \le C(1 + L(\theta))^{\frac{1}{2}}, \quad (1.8)
$$

for some positive constant C.

Let $(\gamma_n)_{n \ge 1}$ be a positive step sequence satisfying the following decreasing step assumption

$$
\sum_{n \ge 1} \gamma_n = \infty \quad \text{and} \quad \sum_{n \ge 1} \gamma_n^2 < \infty. \tag{1.9}
$$

Finally assume that θ₀ is independent of the innovation sequence (U(n))n≥1 and satisfies E[L(θ₀)] < ∞. Then, the sequence (θn)n≥1 defined by the recursive procedure (1.3) satisfies the following properties:

1. the sequence $(\theta_n - \theta_{n-1})_{n \ge 1}$ converges a.s. and in $L^2(\mathbb{P})$ to zero as $n$ goes infinity,

2. the sequence $(L(\theta_n))_{n \ge 1}$ is $L^1(\mathbb{P})$-bounded,

3. $L(\theta_n) \to L_\infty \in L^1(\mathbb{P})$ a.s. as $n \to \infty$.

4. $\sum_{n \ge 1} \gamma_n \langle \nabla L, h \rangle (\theta_{n-1}) < \infty$ a.s.

Proof. The key ingredient of the proof is the convergence theorem for non-negative super-martingale. We only sketch
the proof. Since $\nabla L$ is Lipschitz-continuous, by the fundamental theorem of calculus, one has

$$
\begin{align}
L(\theta_{n+1}) &\le L(\theta_n) + \gamma_{n+1} \langle \nabla L(\theta_n), \theta_{n+1} - \theta_n \rangle + C |\theta_{n+1} - \theta_n|^2 \nonumber \\
&= L(\theta_n) - \gamma_{n+1} \langle \nabla L, h \rangle(\theta_n) - \gamma_{n+1} \langle \nabla L(\theta_n), \Delta M_{n+1} \rangle + C \gamma_{n+1}^2 |H(\theta_n, U^{(n+1)})|^2 \tag{1.10}
\end{align}
$$

where we recall that $\Delta M_{n+1} = H(\theta_n, U^{(n+1)}) - h(\theta_n)$ is an $\{\mathcal{F}_n\}$-martingale increment. Using the above inequality together with assumptions (1.7), (1.8) and (1.9), one proves that $L(\theta_n) \in L^1(\mathbb{P})$ so that $\Delta M_{n+1} \in L^2(\mathbb{P})$, $n \ge 0$. The key idea is to introduce the following sequence

$$
S_n := \frac{L(\theta_n) + \sum_{k=1}^{n} \gamma_k \langle \nabla L, h \rangle (\theta_{k-1}) + C \sum_{k \ge n+1} \gamma_k^2}{\prod_{k=1}^{n} (1 + C \gamma_k^2)}, \quad n \ge 1
$$

with $S_0 = L(\theta_0)$. From assumption (1.7) and using (1.10), one proves that $(S_n)_{n \ge 0}$ is a non-negative super-martingale so that it converges a.s. toward an integrable random variable $S_\infty$. As a by-product of this convergence, one deduces the $L^1$-boundedness of $(L(\theta_n))_{n \ge 0}$ and that $\sum_{n \ge 1} \gamma_n \langle \nabla L, h \rangle (\theta_{n-1}) < \infty$ a.s. Then one proves that $L(\theta_n) \to L_\infty$, a.s. and $\sum_{n \ge 0} \mathbb{E}[|\theta_{n+1} - \theta_n|^2] < \infty$. $\square$

**Remark 1.1.1.** We also want to mention that when the innovation sequence $(U^{(n)})_{n \ge 1}$ is not i.i.d. but only $\{\mathcal{F}_n\}$-adapted, the recursive procedure (1.3) can be written

$$
\theta_{n+1} = \theta_n - \gamma_{n+1}h(\theta_n) - \gamma_{n+1}(\Delta M_{n+1} + r_{n+1})
$$

where $r_{n+1} := h(\theta_n) - \mathbb{E}[H(\theta_n, U^{(n+1)}|\mathcal{F}_n)]$ is an $\mathcal{F}_{n+1}$-measurable remainder term. In this case, assuming that

$$
\sum_{n \ge 1} \gamma_n |r_n|^2 < \infty, \text{ a.s.}
$$

one still proves that the conclusions of the Robbins-Siegmund Lemma hold.
The next result allows one to conclude that the recursive scheme (1.3) converges toward its target $\theta^*$.

**Corollary 1.1.2** (Robbins-Monro & Stochastic gradient algorithms).

(1) [Robbins-Monro algorithm] Assume that the function of interest $h$ satisfies

$$ \forall \theta \neq \theta^*, \quad \langle \theta - \theta^*, h(\theta) \rangle > 0, \qquad (1.11) $$

and

$$ \forall \theta \in \mathbb{R}^d, \quad \mathbb{E}[|H(\theta, U)|^2]^{\frac{1}{2}} \leq C(1 + |\theta|) \qquad (1.12) $$

with $\theta_0 \in L^2(\mathbb{P})$ independent of the innovation sequence $(U^{(n)})_{n \geq 1}$. Assume that the step sequence $(\gamma_n)_{n \geq 1}$ satisfies (1.9). Then, the sequence $(\theta_n)_{n \geq 1}$ defined by (1.3) converges a.s. to $\theta^*$, which satisfies $\{h = 0\} = \{\theta^*\}$.

(2) [Stochastic Gradient algorithm] Assume that $h = \nabla L$ for a convex and continuously differentiable function $L$ satisfying (1.6), $\lim_{|\theta| \to \infty} L(\theta) = \infty$, $\{\nabla L = 0\} = \{\theta^*\}$. Assume that $\mathbb{E}[|H(\theta, U)|^2] \leq C(1 + L(\theta))$ and that $L(\theta_0) \in L^1(\mathbb{P})$. Assume that the step sequence $(\gamma_n)_{n \geq 1}$ satisfies (1.9). Then, $\theta^* = \arg \min_{\theta \in \mathbb{R}^d} h(\theta)$ and $(\theta_n)_{n \geq 1}$ converges a.s. toward $\theta^*$.

**Remark 1.1.2.** There is also an alternative to the stochastic gradient algorithm known as the Kiefer-Wolfowitz (K.W.) procedure. In this approach, the function of interest $h$ is given by $h(\theta) = \mathbb{E}[\hat{\partial}_{\theta}H(\theta, U)]$. In some practical implementations, the local gradient $\hat{\partial}_{\theta}H(\theta, U)$ may be difficult to simulate because of Malliavin weights, log-likelihood ratio, etc, whereas the computation of $H(\theta, U)$ is easy. The K.W. algorithm combines the recursive stochastic approximation algorithm with a finite difference approach to differentiation. Based on the following approximation of the local gradient

$$ \partial_{\theta_i} H(\theta, U) \approx \frac{H(\theta + \eta^i e_i, U) - H(\theta - \eta^i e_i, U)}{2\eta^i} $$

where $(e_i)_{1 \leq i \leq d}$ denotes the canonical basis of $\mathbb{R}^d$ and $|\eta| \ll 1$, $(\eta^i)_{1 \leq i \leq d} \in \mathbb{R}^d$, one introduces the following recursive scheme

$$ \theta_{n+1}^i = \theta_n^i - \gamma_{n+1} \frac{H(\theta_n + \eta_{n+1}^i, U^{(n+1)}) - H(\theta_n - \eta_{n+1}^i, U^{(n+1)})}{2\eta_{n+1}^i}, \quad 1 \le i \le d, \qquad (1.13) $$

where $(\gamma_n)_{n \geq 1}$ is a step sequence satisfying (1.9) and $(\eta_n^i)_{n \geq 1}$, $1 \leq i \leq d$, are $d$ non-negative step sequences satisfying for $1 \leq i \leq d$,

$$ \eta_n^i \to 0, \quad \sum_{n \geq 1} \eta_n^i = \infty \quad \text{and} \quad \sum_{n \geq 1} (\gamma_n / \eta_n^i) < \infty. $$

We refer to the monograph Benveniste and al. [19] for a result about the a.s. convergence of (1.13).

We conclude by the ODE method which establishes a link between the recursive procedure (1.3) and the underlying ODE (1.5) as mentioned above. The main improvement provided by the ODE method is to study the asymptotic behaviour of the sequence $(\theta_n)_{n \geq 0}$ (assumed to be bounded a priori) from a sequence function $(\theta^{(k)})_{k \geq n}$, $n \geq 1$ defined by linear interpolation. Importantly, we will not require the mean reverting assumption (1.11) or the existence of a Lyapunov function that separates the target $\theta^*$. However, the convergence of the procedure will be related to the asymptotic behaviour of the underlying ODE.

We set $\Gamma_0 = 0$ and $\Gamma_n = \sum_{k=1}^n \gamma_k$. Then we define the cádlàg step function $(\theta_t^{(0)})_{t \geq 0}$ by

$$ \theta_t^{(0)} := \theta_n, \quad t \in [\Gamma_n, \Gamma_{n+1}). $$

and the sequence of time shifted function

$$ \theta_t^{(n)} := \theta_{\Gamma_n+t}^{(0)}, \quad t \in \mathbb{R}_+. $$

For every $t \in \mathbb{R}_+$, we also set $N(t) := \max\{k : \Gamma_k \leq t\}$. For every $t \in [\Gamma_n, \Gamma_{n+1})$, the dynamics (1.3) can be written as follows

$$ \theta_t^{(0)} = \theta_0^{(0)} - \int_0^{\Gamma_n} h(\theta_s^{(0)}) ds - \sum_{k=1}^n \gamma_k \Delta M_k $$

which in turn satisfies

$$ \begin{aligned} \forall t \in \mathbb{R}_+, \quad & \theta_t^{(0)} = \theta_0^{(0)} - \int_0^{\Gamma_{N(t)}} h(\theta_s^{(0)}) ds - \sum_{k=1}^{N(t)} \gamma_k \Delta M_k \\ & = \theta_0^{(0)} - \int_0^t h(\theta_s^{(0)}) ds + \int_{\Gamma_{N(t)}}^t h(\theta_s^{(0)}) ds - \sum_{k=1}^{N(t)} \gamma_k \Delta M_k. \end{aligned} $$
Now by the definition of the shifted function $\theta^{(n)}$, one gets

$$
\begin{align*}
\theta_t^{(n)} &= \theta_{\Gamma_n}^{(0)} + \theta_{\Gamma_n+t}^{(0)} - \theta_{\Gamma_n}^{(0)} \\
&= \theta_{\Gamma_n}^{(0)} - \int_{\Gamma_n}^{\Gamma_n+t} h(\theta_s^{(0)})ds + \mathcal{R}^n(t)
\end{align*}
$$

where $\mathcal{R}^n(t) := \int_{\Gamma_{N(\Gamma_n+t)}}^{\Gamma_n+t} h(\theta_s^{(0)})ds - \sum_{k=n+1}^{N(\Gamma_n+t)} \gamma_k \Delta M_k$ is a remainder term.

The next result establishes a connection between the asymptotic behaviour of $(\theta_n)_{n\ge 0}$ and the asymptotic of the sequence of functions $(\theta^{(n)})_{n\ge 0}$.

**Proposition 1.1.3.** Let $L: \mathbb{R}^d \to \mathbb{R}_+$ be a continuously differentiable function such that $\nabla L$ is Lipschitz-continuous, $\lim_{|\theta| \to +\infty} L(\theta) = +\infty$ and satisfying (1.6). If $|h| \le C(1+L)^{\frac{1}{2}}$ and $\langle \nabla L, h \rangle$ is a non-negative lower semi-continuous function then the limiting values of $(\theta_n)_{n \ge 1}$ is a connected component of $\{\langle \nabla L, h \rangle = 0\}$.

1.1.3 Convergence rate

We now present some results related to the convergence rate of the recursive scheme (1.3) in the i.i.d. framework. In standard settings, a stochastic algorithm converges to its target at rate $\sqrt{\gamma_n}$, which suggests to choose $\gamma_n = \gamma/n$, $\gamma > 0$, under assumption (1.9). Under a strong mean-reverting assumption, we will see that this is a sharp non-asymptotic $L^2(\mathbb{P})$, which is consistent with the Monte Carlo method. In Chapter 3, we will see that this rate also appears when one investigates concentration inequalities, that is non-asymptotic deviation estimates between $(\theta_n)_{n\ge 1}$ and its target. Then, we present the standard central limit theorem (CLT) which establishes that the renormalised sequence $\sqrt{\gamma_n^{-1}}(\theta_n - \theta^*)$ converges in distribution to a normal distribution with mean 0 and a covariance matrix based on $\mathbb{E}[H(\theta^*, U)H(\theta^*, U)^t]$.

**Proposition 1.1.4.** Assume that the assumptions of Corollary 1.1.2 (1) (Robbins-Monro algoirthm) are satisfied.
*Assume that there exists *α* > 0 such that the strong mean reverting assumption holds*

$$
\forall \theta \in \mathbb{R}^d, \quad \langle \theta - \theta^*, h(\theta) \rangle \ge \alpha |\theta - \theta^*|. \tag{1.14}
$$

Then, there exists a positive constant C, such that one has

$$
\mathbb{E}[|\theta_n - \theta^*|^2] \le C\gamma_n, \quad n \ge 1,
$$

for the two following step sequences:

• $\gamma_n = \gamma/n^\beta$, $\beta \in (1/2, 1)$, $n \ge 1$.

• $\gamma_n = \gamma/n$, with $\gamma > \frac{1}{2\alpha}$. $n \ge 1$.

The above result suggests a weak rate given by $\sqrt{\gamma_n^{-1}}$. The study of the asymptotic distribution of the renormalised sequence $\sqrt{\gamma_n^{-1}}(\theta_n - \theta^*)$ led to a considerable and expanding literature starting from some pioneering (independent) works by Bouton and Kushner in the early 1980 until recently, see e.g. Fort [43] or Lelong [90]. We give here a result obtained by Pelletier [113] which has as the main advantage to be local in the sense that the CLT holds on the set of convergence of the algorithm to an equilibrium which makes possible a straightforward application to multi-target algorithms.

**Theorem 1.1.5 (Central limit theorem).** Consider the sequence $(\theta_n)_{n\ge 1}$ defined by (1.3) and let $\theta^* \in \{h=0\}$ be one possible target. Assume that the following assumptions hold:

(1) Attractivity: $\theta^*$ is a strong attractor for the underlying ODE, that is h is twice continuously differentiable at $\theta^*$ and Dh($\theta^*$) is an Hurwitz matrix, i.e. all its eigenvalues have positive real parts. We denote by $\lambda_m$ the eigenvalue with the lowest real part.

(2) Regularity of H: The function H satisfies the following regularity and growth control property

$$
\theta \mapsto \mathbb{E}[H(\theta, U)H(\theta, U)^t] \text{ is continuous at } \theta^* \text{ and } \theta \mapsto \mathbb{E}[|H(\theta, U)|^{2+\delta}] \text{ is locally bounded on } \mathbb{R}^d
$$

for some $\delta > 0$.

(3) Non-degenerate asymptotic variance: The covariance matrix

$$
\Sigma^* := \mathbb{E}[H(\theta^*, U)(H(\theta^*, U))^t]
$$

is a positive definite matrix.
(4) Choice of the step sequence: The step sequence is given by

$$
\gamma_n = \frac{\gamma}{a + n^\beta}, \quad a, \gamma > 0
$$

for some $\beta \in (\frac{1}{2}, 1]$, with the additional constraint, when $\beta = 1$,

$$
\gamma > \frac{1}{2\mathcal{R}e(\lambda_m)}. \qquad (1.15)
$$

Then, the a.s. convergence is ruled on the event $\{\theta_n \to \theta^*\}$ by the following stable central limit theorem

$$
\gamma_n^{-1/2} (\theta_n - \theta^*) \Rightarrow \mathcal{N}(0, \Gamma), \qquad (1.16)
$$

where

$$
\Gamma := \int_{0}^{\infty} \left( e^{-(Dh(\theta^*) - \frac{I_d}{2c})s} \right)^t \Sigma^* e^{-(Dh(\theta^*) - \frac{I_d}{2c})s} ds
$$

with $c = \infty$ if $\beta \neq 1$ and $c = \gamma$ if $\beta = 1$. Moreover, the convergence is stable, that is, for every bounded continuous function $f$ and every $A \in \mathcal{F}_\infty$, one has

$$
\mathbb{E}[\mathbf{I}_{\{\theta_n \to \theta^*\} \cap A} f(\sqrt{\gamma_n^{-1}(\theta_n - \theta^*)})] \to \mathbb{E}[\mathbf{I}_{\{\theta_n \to \theta^*\} \cap A} f(\Gamma^{\frac{1}{2}} \zeta)], \quad n \to \infty
$$

where $\zeta \sim \mathcal{N}(0, I_d)$.

We see that the best rate of convergence is achieve for $\beta = 1$ but we have an additional constraint (1.15) on the constant $\gamma$ depending on the unknown matrix $Dh(\theta^*)$. In practical implementations, one does not “observe” such constraint but we remark that the algorithm (1.3) spends most of its time exploring the state space. The CLT regime occurs when the algorithm gets trapped after the mentioned search phase.

In order to circumvent this (theoretical) constraint on the step sequence, one can use the *averaging principle*, also known as Ruppert & Polyak principle in stochastic approximation community, in order to achieve (for free!) the best convergence rate. The original motivation was to smoothen the behaviour of the original procedure by considering the empirical mean of the past values up to the nth iteration rather than the computed value at the nth iteration. If one devices this averaging procedure with a slow decreasing step sequence $(\gamma_n)_{n \ge 1}$ then one reaches for free the best possible rate of convergence.

More precisely, we consider the step sequence $\gamma_n = (\gamma/(a+n))^\beta$, with $\beta \in (1/2, 1)$ and $n \ge 1$ and define the empirical mean

$$
\bar{\theta}_n := \frac{\theta_0 + \dots + \theta_n}{n+1} = \bar{\theta}_{n-1} - \frac{1}{n}(\bar{\theta}_{n-1} - \theta_n), \quad n \ge 1,
$$

of the sequence $(\theta_n)_{n \ge 0}$ defined by (1.3). We observe that by Cesàro's lemma, one has $\bar{\theta}_n \to \theta^*$ a.s. if $\theta_n \to \theta^*$ a.s. as $n \to \infty$. Then, under the assumptions (1), (2), (3) of Theorem 1.1.5, one has

$$
\sqrt{n}(\bar{\theta}_n - \theta^*) \Rightarrow \mathcal{N}(0, \Gamma^*), \quad n \to +\infty
$$

where

$$
\Gamma^* := Dh(\theta^*)^{-1} \Sigma^* (Dh(\theta^*)^t)^{-1}
$$

is the optimal asymptotic covariance matrix.

From a practical point of view, one should not start averaging the initial procedure at the very beginning but wait
until the exploration phase is finished.

We have seen some basic results related to a.s. convergence and convergence rate of stochastic approximation, whilst other interesting aspects have purposely not been cited. For instance, we will mention *traps*, that is parasitic equilibrium such saddle points local maxima, see e.g. Brandière and Duflo [25], Pemantle [114] for some aspects on how an excited enough parasitic equilibrium point is a.s. not a possible limit point for a stochastic approximation algorithm. We also mention stochastic approximation with *averaging innovation* where the innovation sequence $(U^{(n)})_{n \ge 1}$ is no longer i.i.d. but satisfies some mild ergodic or simply averaging assumption, thus allowing for deterministic sequence *Quasi Monte Carlo stochastic approximation* or to plug exogenous data, like market data, see Laruelle & Pagès [87] and the references therein.

*Example: A Robbins-Monro algorithm for quantile estimation*
We conclude this section by a simple illustration of the mentioned results to the recursive computation of the α-quantile (or VaR$_{α}$) of a real-valued random variable X, with law µ, defined as a solution θ* to the equation:

$$
\mathbb{P}(X \le \theta) = \alpha \Leftrightarrow h(\theta) := \mathbb{E}[H(\theta, X)] = 0, \quad \alpha \in (0, 1)
$$

with $H(\theta, x) := \mathbf{I}_{\{x \le \theta\}} - \alpha$ and $h(\theta) = \mathbb{P}(X \le \theta) - \alpha$. Since $X$ is finite, a solution to the above equation always exists. If the probability distribution function of $X$ is increasing, this solution is unique. For sake of simplicity, we will assume that this is case. Let $(\gamma_n)_{n \ge 1}$ be a deterministic positive step sequence satisfying (1.9). The dynamics of the Robbins-Monro algorithm is given by

$$
\theta_{n+1} = \theta_n - \gamma_{n+1} H(\theta_n, X^{(n+1)}), \quad n \ge 0
$$

where $\theta_0 \in L^2(\mathbb{P})$ and is independent of the i.i.d. innovation sequence $(X^{(n)})_{n \ge 1}$ (with common law $\mu$). We introduce the Lyapunov function $L(\theta) := \frac{1}{2}|\theta - \theta^*|^2$ and remark that

$$
\langle \nabla L, h \rangle(\theta) = (\theta - \theta^*)(\mathbb{P}(X \le \theta) - \mathbb{P}(X \le \theta^*)) > 0, \quad \forall \theta \neq \theta^*
$$

since $\theta \mapsto \mathbb{P}(X \le \theta)$ is increasing. Moreover, the function $H$ being bounded assumption (1.8) is clearly satisfied.
Hence, by Corollary 1.1.2 the sequence $(\theta_n)_{n \ge 1}$ converges a.s. toward its target $\theta^*$.

In order to investigate the weak rate of convergence, we assume that $X$ admits a continuously differentiable density $p_X$ satisfying $p_X(\theta^*) > 0$. Our aim is to apply Theorem 1.1.5. Assumptions (1) and (2) are clearly satisfied. Moreover, by the very definition of $\theta^*$, one has $\Sigma^* = \mathbb{E}[H(\theta^*, X)^2] = \alpha(1-\alpha) > 0$. In order to achieve the optimal asymptotic rate we set $\gamma_n = \gamma/(a+n)$, with $a>0$ and $\gamma > 1/(2p_X(\theta^*))$. From Theorem 1.1.5, one gets

$$
\sqrt{n}(\theta_n - \theta^*) \Rightarrow \mathcal{N}(0, \Gamma), \quad n \to \infty
$$

with

$$
\Gamma := \alpha(1 - \alpha) \frac{\gamma^2}{2\gamma p_X(\theta^*) - 1}.
$$

This asymptotic variance is minimal for $\gamma = \frac{1}{p_X(\theta^*)}$ with an optimal value $\Gamma = \alpha(1-\alpha)/p_X^2(\theta^*)$. However, in practical implementations, one does not have access to the value of $\theta^*$ or the density of $X$. As a consequence, one is naturally led to consider the averaged form of the algorithm although it is completely satisfactory for a practical point of view. When one tries to implement the stochastic approximation scheme $(\theta_n)_{n \ge 1}$, one observes a chaotic behaviour with a slow convergence, specially when $\alpha \approx 1$, say $\alpha > 95\%$. In [1] and [2], we proposed a R.M. algorithm to compute the quantile and the expected shortfall and improved its performance by combining it with a stochastic approximation based importance sampling technique.

In Chapter 2, we propose an application to the recursive computation of the shortfall $L^p$-risk of a financial portfolio and the recursive computation of risk minimising strategies.

We refer the interested reader to the Phd dissertation of Laruelle [86] for some applications of the stochastic approximation theory to the optimal splitting of orders across liquidity pools such that dark pools, to multi-arm clinical trials by urn models. We also mention the monograph of Pagès [110] for stochastic approximation based procedures for parameter implicitation and calibration of financial models, or the computation of optimal quantizers by e.g. the Competitive Learning Vector Quantization (CLVQ) method.

1.2 Discretisation and simulation of stochastic processes

1.2.1 Motivation

In this section, we consider an $\mathbb{R}^d$-valued stochastic differential equation (SDE in short) with dynamics

$$
X_t = x + \int_0^t b(s, X_{s-}) ds + \int_0^t \sigma(s, X_{s-}) dZ_s, \quad t \ge 0, \qquad (1.17)
$$

where the coefficients $b: \mathbb{R}_+ \times \mathbb{R}^d \to \mathbb{R}^d$, $\sigma: \mathbb{R}_+ \times \mathbb{R}^d \to \mathbb{R}^d \otimes \mathbb{R}^q$ are bounded measurable in time, Lipschitz-continuous in space (uniformly in time) and where $(Z_t)_{t \ge 0}$ is a $q$-dimensional Lévy process with Lévy-Khintchine decomposition given by

$$
\mathbb{E}[\exp(i\langle u, Z_t \rangle)] = \exp(t\psi(u)), \quad \psi(u) = \exp\left(-\frac{|u|^2}{2}\right) + \int_{\mathbb{R}^q} (\exp(i\langle u, z \rangle) - 1 - i\langle u, z \rangle \mathbf{I}_{\{|z| \le 1\}})\nu(dz). \quad (1.18)
$$
In other words, the SDE (1.17) is driven by a process Z that is given by the sum of a q-dimensional Brownian motion and pure jump process N of Lévy measure ν.

Such dynamics appears in various fields. This is so because of the strong connection between SDEs and PDEs (see e.g. Bass [14], Friedman [45]) via Feynman-Kac representation formulas. In various situations, the unique classical solution of a linear parabolic or elliptic PDE can be expressed as the expectation of a (path) functional of the SDE (1.17).

When Z is a q-Brownian motion, the diffusion process (1.17) is commonly used to provide probabilistic representation of the heat equation in non-homogeneous media. It is also commonly used to model the dynamics of assets in financial and energy markets. In [4], we proposed a joint model of Gas and electricity spot prices based on Ornstein processes with parameterized diffusion coefficients that captures the most salient statistical properties observed on data. When Z is a pure jump process, the dynamics (1.17) allows to reproduce some stylised features observed in energy markets such that spikes in the spot prices, see e.g. Benth and al. [18]. In mathematical finance, the price of an option is usually seen as the solution of a parabolic PDE in which the terminal condition is given by the payoff.

In order to compute a quantity related to the law of the process with dynamics (1.17), which is not known in most applications, one usually considers a discrete approximation scheme. A natural scheme that can be implemented very easily is the Euler-Maruyama scheme that we now describe. On the interval $[0, T]$, $T > 0$, for a given time step $h = T/N > 0$, $N \in \mathbb{N}^*$, we introduce the time grid $t_i = ih$, $i = 0, \dots, N$ and the continuous Euler scheme $(X_t^h)_{t \in [0,T]}$ defined for $t \in [0, T]$ by

$$X_t^h = x + \int_0^t b(\phi(s), X_{\phi(s)}^h)ds + \int_0^t \sigma(\phi(s), X_{\phi(s)}^h)dZ_s, \quad t \ge 0, \quad (1.19)$$

where $\phi(s) = \sup\{t_i : t_i \le s\}$. In other words, the Euler scheme is obtained by freezing the coefficients to the current discretisation point in each small interval $[t_i, t_{i+1}]$. Let us mention that in the case where Z has Brownian and jump parts, it is not always possible to simulate the increments $(Z_t - Z_s)_{0 \le s \le t \le T}$ of the innovation process so that it is not always possible to implement the scheme (1.19). However, it is still possible to provide an analysis of the error, see Section 1.2.3. In many important cases, one can implement this scheme for instance when Z is a stable process or Gamma process, etc.

We will present two different results related to the approximation of $(X_t)_{t \in [0,T]}$ solution to the SDE (1.17) by the scheme $(X_t^h)_{t \in [0,T]}$ given by (1.19). Approximation of stochastic processes goes back to the works of Stroock and Varadhan. They investigated approximation of the distribution of solutions to martingale problems by Markov chains. We refer the interested reader to Chapter 11 in [120] and to [119] for Dirichlet boundary problems.

The pionnering works of Milstein [103], [104] investigated the strong approximation error, that is the $L^p$-error between the path of $(X_t)_{t \in [0,T]}$ and the path of its approximation $(X_t^h)_{t \ge 0}$. Then, Milstein [105] and Talay [121] proved that the weak error, that is the error between the law of the random variable $X_T$ and the law of $X_T^h$, is of order h under some regularity assumptions.

### 1.2.2 Strong error

The strong error associated to the process (1.17) and its approximation scheme (1.19) is the quantity $E[\sup_{0 \le t \le T} |X_t - X_t^h|^p]^{1/p}$ for $p > 0$. When the underlying innovation process Z is an $\mathbb{R}^q$-valued Brownian motion, one has the following result:

**Proposition 1.2.1.** Suppose that the coefficients b and σ of the SDE (1.17) satisfy the following regularity condition:

$$\forall (s, t) \in [0, T]^2, \quad \forall (x, y) \in \mathbb{R}^d, \quad |b(t, x) - b(s, y)| + \|\sigma(s, x) - \sigma(t, y)\| \le C(|t - s|^\alpha + |x - y|) \quad (1.20)$$

where $\alpha \in [0, 1]$ and for some positive constant C. Then, for all $p > 0$, there exists a positive constant $C_p$ such that

$$E[\sup_{0 \le t \le T} |X_t - \bar{X}_t^h|^p]^{1/p} \le C_p(1 + |x|) h^{\alpha \wedge \frac{1}{2}}. \quad (1.21)$$

As a consequence, we remark that the Hölder regularity time exponent $\alpha$ rules the strong convergence rate of the Euler approximation scheme as soon as $\alpha < 1/2$. The key ingredient for the proof of Proposition 1.2.1 is to use a Gronwall Lemma to upper-bound $E[\sup_{0 \le t \le T} |X_t - \bar{X}_t^h|^p]^{1/p}$ by the $L^p$-norm of the increments $X_s - X_{\phi(s)}$ of the diffusion. As a straightforward corollary, we deduce that the Euler approximation scheme (1.19) converges a.s. toward the diffusion process at a rate $\beta \in [0, \alpha \wedge \frac{1}{2})$, namely

$$\forall \beta \in [0, \alpha \wedge \frac{1}{2}), \quad h^{-\beta} \sup_{0 \le t \le T} |X_t - \bar{X}_t^h| \to 0, \quad \text{a.s.}$$

as $h \downarrow 0$. The proof follows from the $L^p$-convergence Proposition 1.2.1 and an argument using Borel-Cantelli Lemma.
In order to compensate this slow convergence rate Milstein introduced in [104] a new scheme in order to achieve the usual strong error of order $h$ that appears in the deterministic framework. In dimension $d = q = 1$ its expression is simple and it can be easily implemented provided that $b$ and $\sigma$ have enough regularity. In higher dimension, some theoretical and simulation problems make its use more questionable. The key idea is to use an Itô-Taylor formula in the stochastic integral between to discrete times of the grid (which behaves like $\sqrt{h}$ and is the leading term in the strong approximation error) in order to introduce higher order terms in the dynamics of the scheme (1.19). We refer to the book of Kloeden and Platten [71] for an overview about higher order schemes.

Similar arguments can be used in order to quantify the strong error with the same rate of convergence when the innovation process $Z$ is a square integrable Lévy process, see e.g. Kohatsu-Higa and Protter [72]. We also mention the work of Bruti Liberati and Platen [27] where the authors investigate the strong approximation error for pure jump SDEs, propose jump-adapted time discretizations and higher order scheme (in the strong sense) when the intensity of the jump process is finite $\nu(\mathbb{R}^q) < \infty$.

### 1.2.3 Weak error

In many situations, like the pricing of European options in mathematical finance, an Euler discretisation scheme $(X_t^h)_{t \in [0,T]}$ such that (1.19) is introduced to evaluate quantities writing as $\mathbb{E}[F((X_t)_{t \in [0,T]})]$ where $F: \mathcal{C}(\mathbb{R}_+, \mathbb{R}^d) \to \mathbb{R}$ is a path functional. In order to evaluate such quantities, one is led to use the Monte Carlo estimator

$$E_M^h := \frac{1}{M} \sum_{i=1}^{M} F((X_t^{h,i})_{t \in [0,T]}), \quad M \ge 1$$

where $((X_t^{h,i})_{t \in [0,T]})_{1 \le i \le M}$ are i.i.d. samples of the path of the Euler scheme $(X_t^h)_{t \in [0,T]}$ defined by (1.19).

The global error between $\mathbb{E}[F((X_t)_{t \in [0,T]})]$, the quantity to estimate, and its implementable approximation $E_M^h$ can be decomposed as follows:

$$\begin{align}
\mathcal{E}_{glob} &:= (\mathbb{E}[F((X_t)_{t \in [0,T]})] - \mathbb{E}[F((X_t^h)_{t \in [0,T]})]) + (\mathbb{E}[f((X_t^h)_{t \in [0,T]})] - E_M^h) \tag{1.22} \\
&:= \mathcal{E}_W^h + \mathcal{E}_S^M. \nonumber
\end{align}$$

The term $\mathcal{E}_S^M$ corresponds to the statistical error. If $F((X_t^h)_{t \in [0,T]}) \in L^2(\mathbb{P})$, the statistical error is of order $\sigma^h/\sqrt{M}$ in $L^2(\mathbb{P})$ and the renormalised statistical error satisfies a CLT. One may also be interested in non-asymptotic deviation estimates of this quantity in order to device non-asymptotic confidence interval. We will come back to this particular topic in Chapter 3.

The term $\mathcal{E}_W^h$ corresponds to the weak error and has been widely investigated in the literature since the seminal work of Talay and Tubaro [122] for the case $F((X_t)_{t \in [0,T]}) = F(X_T)$. One may be interested in quantifying this error for other type of functionals. We will briefly present some results related to the two functionals:

$$1. F((X_t)_{t \in [0,T]}) = F(X_T)$$

$$2. F((X_t)_{t \in [0,T]}) = F(\tau \wedge T, X_{\tau \wedge T}), \text{ where } \tau := \inf \{t \ge 0 : X_t \notin D\} \text{ for a domain } D \subset \mathbb{R}^d.$$

The first case corresponds to the standard Cauchy problem since the quantity to estimate is related to the unique classical solution of linear parabolic PDE with terminal function $F$.

The second case corresponds to the Cauchy-Dirichlet problem since the quantity to estimate is related to the solution of linear parabolic PDE with Dirichlet conditions.

Other type of functional can be discussed such that $F(\max_{0 \le t \le T} X_t, X_T)$ or $F(\int_0^T \mathbf{I}_{\{X_t \ge 0\}} dt, X_T)$ both in dimension $d = q = 1$.

Still in the one-dimensional diffusion case $d = q = 1$, we mention the works of Alfonsi, Jourdain and Kohatsu-Higa [2], [3] where the Wasserstein distance between the laws of $(X_t)_{t \in [0,T]}$ and $(X_t^h)_{t \in [0,T]}$ defined by

$$W_1(\mathcal{L}(X), \mathcal{L}(X^h)) := \sup_{F: \mathcal{C}([0,T], \mathbb{R}) \to \mathbb{R}, \text{ Lip}(F) \le 1} |\mathbb{E}[F(X^h)] - \mathbb{E}[F(X)]|$$

is investigated by means of a pathwise optimal transport technique. Here $\text{Lip}(F)$ denotes the Lipschitz constant of $F$ (with respect to the sup norm on $\mathcal{C}([0,T], \mathbb{R}))$. It is proved that this quantity is upper-bounded by $h^{\frac{2}{3}-\epsilon}$, $\forall \epsilon > 0$, which is an intermediate rate between the strong error and the usual weak error rates. In the multi-dimensional diffusion case, when the coefficients $b$ and $\sigma$ are time-homogeneous $\mathcal{C}^2$, bounded together with their derivatives up to the order 2 and uniform ellipticity holds, the $\rho$-Wasserstein distance is shown to be of order $\sqrt{\log(N)/N}$.
The Cauchy problem for diffusion processes

For the standard continuous Euler scheme (1.19), the weak discretisation error is of order h, provided the coefficients b, σ and the test function F are C 4 with the same polynomial growth assumption. The key idea [121], [122] consists in observing that the function u(t, x) = E[F(X_T)|X_t = x] = E[F(X_T^{t,x})] (note that (X_t)_{t \in [0,T]} is a Strong Markov process under this assumption) is a classical solution to the parabolic linear PDE

$$
\left\{
\begin{array}{l}
(\partial_t + \mathcal{L}_t)u = 0, \quad (t,x) \in [0, T) \times \mathbb{R}^d \\
u(T, x) = F(x),
\end{array}
\right.
\tag{1.23}
$$

where $\mathcal{L}_t$ is the infinitesimal generator associated to (1.17), namely $\mathcal{L}_t F(x) = \langle b(t,x), \nabla F(x) \rangle + \frac{1}{2} \mathrm{Tr}(a(t,x) D_x^2 F(x))$, for $F \in C_0^\infty(\mathbb{R}^d)$, where $a(t,x) = (\sigma\sigma^*)(t,x)$. Under the mentioned assumption on the coefficients and test function $F$, one gets that $u \in C^{2,4}([0,T] \times \mathbb{R}^d)$ with the following controls on the derivatives: $\forall (t,x) \in [0,T] \times \mathbb{R}^d$, $|\partial_t^i \partial_x^j u(t,x)| \le C_T(1+|x|^{p_{i,j}})$. We remark that the last controls allow to deduce that $|\partial_t^i \partial_x^j u(t_i, X_{t_i}^h)| \in L^p(\mathbb{P})$ for any $p \ge 1$. The next step is to introduce the solution of the PDE (1.23) as follows

$$
\mathcal{E}_W^h := \mathbb{E}[F(X_T)] - \mathbb{E}[F(X_T^h)] = -(\mathbb{E}[u(T, X_T^h)] - u(0, x)) = - \sum_{i=0}^{N-1} \mathcal{E}_W^{h,i}
$$

with $\mathcal{E}_W^{h,i} := \mathbb{E}[u(t_{i+1}, X_{t_{i+1}}^h) - u(t_i, X_{t_i}^h)]$. From these representation, one applies Itô-Taylor formula around the point $(t_i, X_{t_i}^h)$ for $i = 0, \dots, N-1$ using the regularity of $u$. In order to prove that $\mathcal{E}_W^{h,i} = \mathcal{O}(h^2)$, it is necessary to obtain an expansion of order two in time and 4 in space, the term of order 3 in space vanishes because of the third moment of the underlying noise is zero.

Let us also mention here that when the Gaussian increments of the scheme are replaced by more general (possibly
discrete) random variables having the same covariance matrix and odd moments up to order 5 as the standard Gaussian
vector of $\mathbb{R}^q$, it can be checked that the error expansion at order $h$ still holds for $b, \sigma, F$ smooth enough.

If we keep on that direction by doing an expansion at order 3 in time and 6 in space, observing that the odd moments of the underlying Gaussian vector vanishes, one gets

$$
\mathcal{E}_W^h = C_T^1 h + \mathcal{O}(h^2),
$$

where $C_T^1$ is an explicit constant independent of $h$. According to Talay and Tubaro [122], when $b, \sigma$ are $C^\infty$ with bounded derivatives of all orders and $F$ is $C^\infty$ with polynomial growth together with its derivatives, then for each integer $L \ge 1$, the expansion

$$
\mathcal{E}_W^h = \sum_{\ell=1}^{L} C_T^\ell h^\ell + \mathcal{O}(h^{L+1}) \quad (1.24)
$$

holds.

On can take advantage of the previous expansion in order to perform a *Richardson-Romberg extrapolation*. This
idea was originally introduced in [122] and then deeply investigated by Pagès in [109]. One considers the strong
solution $X^2$ of a “copy” of equation (1.17), driven by a second Brownian motion $W^2$ defined on the same probability
space $(\Omega, \mathcal{A}, \mathbb{P})$. One may always choose such a Brownian motion by enlarging the space $\Omega$ if necessary.

Then one considers the continuous Euler scheme with a twice smaller step h/2, denoted by $(X_t^{2,h/2})_{t \in [0,T]}$ of the diffusion $X^2$. Then combining the two time discretisation error expansions associated to $X^{1,h}$ (related to the first Brownian motion $W^1$) and $X^{2,h/2}$, we get

$$
\mathbb{E}[F(X_T)] = \mathbb{E}[2F(X_T^{2,h/2}) - F(X_T^{1,h})] - \frac{C_T^2}{2} h^2 + \mathcal{O}(h^3)
$$

so that the global quadratic error becomes

$$
\mathbb{E}\left[\left(\mathbb{E}[F(X_T)] - \frac{1}{M} \sum_{i=1}^{M} 2F(X_T^{2,h/2,(i)}) - F(X_T^{1,h,(i)})\right)^2\right] = \left(\frac{C_T^2}{2} h^2\right)^2 + \frac{\mathrm{Var}(2F(X_T^{2,h/2}) - F(X_T^{1,h}))}{M} + \mathcal{O}(h^5)
$$

If one follows a naive approach by simulating $X^{2,h/2}$ and $X^{1,h/2}$ independently, this corresponds to two independent Brownian motion $W^1$, $W^2$, then one increases the variance by a factor 5 (asymptotically) since $\mathrm{Var}(2F(X_T^{2,h/2}) - F(X_T^{1,h})) = 4\mathrm{Var}(2F(X_T^{2,h/2}) + \mathrm{Var}(F(X_T^{1,h}))) \to 5\mathrm{Var}(F(X_T))$ as $h \downarrow 0$. It is shown in [109] that the consistent choice $W^1 \equiv W^2$ leading to consistent Brownian increments for the two schemes is asymptotically optimal among all possible choice of Brownian motions $W^1$ and $W^2$. We will see in Chapter 4 an extension of the Richardson-Romberg extrapolation method to stochastic optimisation by means of stochastic approximation algorithm.
Anyway in order to analyse the weak discretisation error we will always proceed this way. As mentioned above the key ingredient is to have controls on $u$ and its derivatives and to proceed using Itô-Taylor's formula. This is for that reason that we require regularity on $b$, $\sigma$ and $F$. However, in practical situations the test function $F$ is not regular. For instance, in mathematical finance, one often deals with payoff $F$ that are not regular, e.g. Call option $F(x) = (x - K)_+$ or digital option $F(x) = \mathbf{I}_{\{x \ge K\}}$.

One can weaken the regularity assumption on $F$ if one can benefit from the smoothness of the underlying transition density. More precisely, if the process $(X_t)_{t \in [0,T]}$ solution to (1.17) admits a smooth transition density, one can relax the smoothness assumption on $F$. Indeed, in this case one has $u(t,x) = \int_{\mathbb{R}^q} F(y)p(t, T, x, y)dy$ where $y \mapsto p(t, T, x, y)$ is the density function of the random vector $X_T^{t,x}$. When the coefficients are time homogeneous, the density is given by $p(T - t, x, y)$.

However, in order to benefit from this regularising property one has to assume that the operator $\mathcal{L}$ satisfies some non-degeneracy assumption. One usually considers two type of assumptions: the Hörmander assumption (H) and uniform ellipticity assumption (UE) of $a = \sigma\sigma^t$.

The Hörmander condition is a sufficient (and is known to be nearly necessary) condition for a second-order linear Kolmogorov PDE with smooth coefficients to be hypoelliptic, that is, $P_t C_b(\mathbb{R}^d, \mathbb{R}) \subset C_b^\infty(\mathbb{R}^d, \mathbb{R})$, for all $t \in (0, \infty)$, where $(P_t)_{t \ge 0}$ denotes the associated Markov semigroup. Note that in the non-hypoelliptic regime, even in the case of smooth coefficients, there exist counterexamples to regularity preservation. In Hairer, Hutzenthaler, Jentzen [59], the roughening effect of the noise for non-hypoelliptic linear PDE is exhibited. From the perspective of numerical probability, this has the consequence that numerical approximations may converge (in the strong and weak sense) without any arbitrarily small polynomial rate.

**Assumption (H):** The coefficients $b$ and $\sigma$ are $C^\infty$ and time homogeneous, with bounded derivatives of order larger than 1 (importantly the functions $b$ and $\sigma$ are not supposed to be bounded). Let $A_0 = b - \frac{1}{2}\sigma \nabla\sigma$ and $A_i = \sigma_i$, $i = 1, \dots, q$ where $\sigma = (\sigma_1, \dots, \sigma_q)$. For a multi-index $\alpha = (\alpha_1, \dots, \alpha_r) \in \{0, \dots, q\}^r$ we define the vector fields $A_i^\alpha$ by $A_i^\emptyset = A_i$ for $1 \le i \le q$ and $A_i^{(\alpha,j)} = [A_j, A_i^\alpha]$ for $0 \le j \le q$ where for two smooth $\mathbb{R}^d$-valued vector fields $V, W$, $[V, W] = \nabla W V - \nabla V W$, $\nabla V = (\partial_j V^i)_{1 \le i,j \le d}$, stands for their Lie bracket. The Hörmander condition holds if the vector space spanned by $A_i^\alpha$, $1 \le i \le q$, $\alpha$ being a multi-index, at the point $x$ is $\mathbb{R}^d$.

Under (H), the density function $y \mapsto p(T - t, x, y)$ of $X_T^{t,x}$ exists and is smooth. Moreover, if the vector fields $A_0, \dots, A_q$ satisfies the uniform Hörmander assumption (UH), that is if there exists an integer $L$ and $\lambda > 0$ such that

$$ \forall x \in \mathbb{R}^d, \quad \inf_{u \in S^{d-1}} \sum_{|\alpha| \le L} \langle A^\alpha(x), u \rangle^2 \ge \lambda \quad (1.25) $$

Kusuoka and Stroock [83] provided the following controls on the density. For any integers $m, k$ and any multi-index $\alpha$ and $\beta$ such that $2m + |\alpha| + |\beta| \le k$, there exist an integer $M$ and a non-decreasing function $T \mapsto C(T)$ such that the following inequality holds:

$$ |\hat{c}_t^m \hat{c}_x^\alpha \hat{c}_y^\beta p(t, x, y)| \le \frac{C(T)(1 + |x|^M)}{t^q(1 + |y - x|^2)^k} \exp(-c \frac{|y - x|^2}{t(1 + |x|^2)}). \quad (1.26) $$

Moreover, the inverse of the Malliavin covariance matrix denoted by $\Gamma_t(x)$ is bounded in any $L^p(\mathbb{P})$. More precisely, for any $p \ge 1$, there exists a non-decreasing function $C(T)$ and an integer $k$ such that

$$ \mathbb{E}[\|\Gamma_t\|^p]^{\frac{1}{p}} \le C(T) \frac{1 + |x|^k}{t^{dL}} \quad (1.27) $$

In [10], Bally and Talay take benefit from (1.26) and (1.27) in order to extend (1.24) for hypo-elliptic diffusions without any regularity on the test function $F$, which is supposed to be measurable with exponential growth at infinity. The main quantities that appears in the expansion of the weak discretisation error are $Q_i^h = \mathbb{E}[\hat{c}_x^\alpha u(t_i, X_{t_i}^h)G(X_{t_i}^h)]$ for a smooth function $G$ with polynomial growth. Roughly speaking, when $t_i$ is close to 0, say $i \le N/2$, then one uses (1.26) in order to control $Q_i^h$. When $t_i$ is close to $T$, one has to use Malliavin's integration by parts formula with respect to the random variable $X_{t_i}^h$ which is expected to be non-degenerated with high-probability. We are led to compare $\Gamma_t^h$ and $\Gamma_t$, where $\Gamma_t^h$ is the inverse of the Malliavin covariance matrix of $X^h$ at time $t$. This can be done using a localisation argument.

Let us also mention that in general the Euler scheme does not admit a density under assumption (UH). In [11], Bally and Talay introduced a perturbed Euler scheme which admits a transition density $p^h$ and proved an expansion of the error $p - p^h$.

**Assumption (UE):** The diffusion matrix $a$ is bounded and uniformly elliptic, that is there exists $\bar{a}$ such that $\forall (t, x) \in [0, T] \times \mathbb{R}^d, \forall u \in \mathbb{R}^d, \bar{a}^{-1}|u|^2 \le \langle a(t, x)u, u \rangle \le \bar{a}|u|^2$.

When $b$ and $\sigma$ are Hölder continuous in space (uniformly in time) and bounded, the unique (weak) solution to (1.17) (driven by a Brownian motion $W$) satisfies the following Gaussian bound $\forall T > 0, \exists C > 1, \forall (t, x, y) \in$
$[0, T) \times \mathbb{R}^d \times \mathbb{R}^d, |p(t, T, x, y)| \le C \exp(-C^{-1}|y - x|^2/t)$. The latter bound follows from the *parametrix method*, see Friedman [44] or McKean and Singer [95]. The parametrix technique is a classical perturbation method used in PDE theory that was originally formulated in order to give an expansion in infinite series of iterated kernels of the fundamental solution of an elliptic or parabolic PDE of order $2p, p \in \mathbb{N}^*$, with time inhomogeneous coefficients. More precisely, one has

$$p(t, T, x, y) = \sum_{n \ge 0} (\tilde{p} \otimes H^n)(t, T, x, y) \qquad (1.28)$$

where $\tilde{p}(t, T, x, y)$ is the Gaussian density of the system (1.17) where we removed the drift $b$ and fixed the diffusion coefficient to the terminal point $y$. This system is often referred to the *frozen* or *parametrix* process. The space-time convolution operator $\otimes$ is defined by $(f \otimes g)(t, T, x, y) := \int_t^T ds \int_{\mathbb{R}^d} dz f(t, s, x, z)g(s, T, z, y)$ and $H$ is a kernel. The keystone of this approach is that the underlying parametrix kernel $H$ enjoys a regularising property in the sense that $(\tilde{p} \otimes H^n)(t, T, x, y)$ is explicitly controlled by the term of an absolutely (and uniformly) convergent series times a Gaussian density. We will come back to this point latter on in the second part of this manuscript in Chapter 5.

This method has been successfully applied to many equations and various situations. Its merit and success are due to its flexibility as it can be invoked for a wide variety of PDEs both for theoretical goals such as density estimates [34], martingale problems [97], strong existence and uniqueness for SDE [29] and for numerical approximations [75]. Let us be more precise on the latter topic since it will be used in Chapter 5. If we replace in the dynamics of the continuous Euler scheme (1.19) the Brownian increments by a sequence of *i.i.d.* random variables with common law $\mu$ which admits a density with respect to the Lebesgue measure (with the same moments as the Gaussian density up to a certain order), then the new scheme admits a transition density $p^N$ which admits a representation in infinite series

$$p^N(t, T, x, y) = \sum_{n \ge 1} (\tilde{p} \otimes_N H^{N,n})(t, T, x, y) \qquad (1.29)$$

where $y \mapsto \tilde{p}^N(t, T, x, y)$ is the density of the recursive scheme where we removed the drift $b$ and fixed the diffusion coefficient to the terminal point $y$. The kernel $H^N$ and the operator $\otimes_N$ are discrete version of the kernel $H$ and the space time convolution operator $\otimes$. In [75], Konakov and Mammen followed these approach. By comparing the two series (1.28) and (1.29) they obtained local limit theorems for the difference $(p - p^N)(t, T, x, y)$ at a rate $h^{1/2}$.

When one considers the true continuous Euler scheme with Brownian increments, an expansion of the error for $p - p^N$ can be obtained. More precisely, Konakov and Mammen in [76] proved the following expansion

$$\forall L \in \mathbb{N}^*, \quad (p - p^N)(t, T, x, y) = \sum_{\ell=1}^{L} \frac{\pi_{\ell}^N(t, T, x, y)}{N^{\ell}} + \frac{1}{N^{L+1}} \mathcal{R}_N(t, T, x, y) \qquad (1.30)$$

with the following Gaussian controls

$$\exists C_{T-t} > 1, \forall (x, y) \in (\mathbb{R}^d)^2, |\pi_{\ell}^N(t, T, x, y)| + |\mathcal{R}_N(t, T, x, y)| \le C_{T-t} \exp(-c \frac{|y-x|^2}{T-t})$$

under the assumptions (UE) and assuming that the coefficients $b, \sigma$ are $C^\infty$ with bounded derivatives of all order. We will present in Chapter 5 some results related to the weak approximation error of a skew diffusion with bounded measurable drift and Hölder diffusion coefficient by an Euler-type scheme. We obtained a similar error bound for the difference between the densities of the skew diffusion and its Euler approximation under mild smoothness assumption on the coefficients.

Under stronger smoothness assumption on the coefficients, Gobet and Labart [53] quantified the small time asymptotics of the difference $p - p^N$. More precisely, using Malliavin Calculus techniques, they established the bound $|(p - p^N)(t, T, x, y)| \le C h(T - t)^{-\frac{d+1}{2}} \exp(-c|y - x|^2/(T - t))$, where $C, c > 1$ are constants independent of $t, T$.

The main advantage of (1.30) compared to the previous approach is that one can considerably weaken the assumption on the test function $F$. Guyon [58] proved that (1.24) holds if $F$ is a tempered distributions, e.g. a Dirac mass or derivatives of a Dirac mass. We will see in Chapter 4 that this point is very important in order to obtain an expansion of the implicit discretisation error for solutions to inverse problems. It is also important to quantify this kind of error when one considers the computation of greeks in option hedging.

Recent research works have shown a rising interest in weak higher order approximation schemes starting with the pioneering work of Kusuoka [82]. The method relies on the splitting of the infinitesimal generator in order to use composition techniques of ODEs flow see Ninomiya and Victoir [107] and very recently Bally and Rey [9]. In the latter work, the authors obtained estimates of the error in total variation distance for the Ninomiya Victoir scheme using Malliavin calculus techniques under a Doeblin's condition on the underlying noise.

An interesting problem, which has not been investigated (to the best of our knowledge), will be to study the existence of a transition density for higher order scheme like the Ninomiya scheme (as well as Gaussian lower and upper-bounds) and to see if expansions similar to (1.30) holds. For instance, one may try to follow the approach initiated by Konakov and Mammen [76].
The Cauchy problem in the case of SDE driven by jump process

The approach of Talay and Tubaro [122] still works when one adds jumps in the innovation process Z appearing in the dynamics (1.17). The key idea still consists in controlling the derivatives of the function u, the unique classical solution to linear integro-differential PDE

$$
\left\{
\begin{array}{l}
(\partial_t + \mathcal{L}_t)u(t,x) + \\
\quad \int_{\mathbb{R}^q} (u(t,x + \sigma(t,x)z) - u(t,x)) - \langle \nabla u(t,x), \sigma(t,x)z \rangle \mathbf{I}_{\{|z| \le 1\}} = 0, \quad (t,x) \in [0,T) \times \mathbb{R}^d \\
u(T,x) = F(x),
\end{array}
\right.
\tag{1.31}
$$

Like in the Brownian case, under regularity assumptions on the coefficients and terminal condition F and under integrability condition of the Lévy measure, Protter and Talay [115] derived an expansion of the weak discretisation error similar to (1.24).

One important issue that does not appear in the Brownian framework is the simulation of the increments of the Lévy process Z. In many cases such that stable, Gamma, compound poisson process, one knows how simulate the increments but in general one does not know how to do. Jacod and al. [66] proposed an implementable Euler scheme where the increments $(Z_{t_{i+1}} - Z_{t_i})_{0 \le i \le N-1}$ are approximated by a sequence $(U_i^h)_{0 \le i \le N-1}$ of i.i.d random variables satisfying $|\mathbb{E}[g(Z_h) - g(U_1^h)]| \le C h \epsilon_h$ for $g \in C^4(\mathbb{R}^q)$. Under some integrability conditions on $\nu$, they obtained the bound $|\mathbb{E}[F(X_T) - F(X_T^h)]| \le C(\epsilon_h \vee h)$ and derived an expansion of the error similar to (1.24) under stronger smoothness assumption on the coefficients b and $\sigma$.

When the innovation is a non-degenerated stable process, Konakov and Menozzi in [77] extended the parametrix
machinery developed in [76] and successfully derived an error expansion similar to (1.30).

Let us also mention the work of Kohatsu-Higa and Tankov [74] where the authors study the weak approximation
for a jump-adapted discretisation scheme of Lévy driven SDEs. The approach consists in building upon adaptive
non-uniform discretisation based on the jump times of the driving process coupled with suitable approximations of
the solutions between these jump times.

The Cauchy-Dirichlet problem for diffusion processes

For a Brownian driven SDE with dynamics (1.17), let us consider the first exit time $\tau$ of the domain $D \subset \mathbb{R}^d$ defined by $\tau = \inf \{t \ge 0 : X_t \notin D\}$. The computation of $\mathbb{E}[F(\tau \wedge T, X_{\tau \wedge T})]$ is an important issue that appears in various fields such as mathematical finance for the pricing and hedging of look-back options. More generally, it is the probabilistic representation of the classical solution of the following Cauchy-Dirichlet problem

$$
\left\{
\begin{array}{ll}
(\partial_t + \mathcal{L}_t)u(t,x) = 0, & (t,x) \in [0,T) \times D, \\
u(t,x) = F(t,x), & (t,x) \in [0,T) \times \partial D \cup \{T\} \times \bar{D}
\end{array}
\right.
\tag{1.32}
$$

One is naturally led to consider the approximation $\mathbb{E}[F(\tau^h \wedge T, \Pi_{\bar{D}}(X_{\tau^h \wedge T}^h))]$ where $\tau^h := \inf \{t_i : X_{t_i}^h \notin D\}$ and $\Pi_{\bar{D}}$ is the projection on $\bar{D}$. Then, using the same decomposition as in the standard Cauchy problem, one has

$$
\mathcal{E}_W^h = \mathbb{E}[F(\tau^h \wedge T, \Pi_{\bar{D}}(X_{\tau^h \wedge T}^h))] - \mathbb{E}[F(\tau \wedge T, X_{\tau \wedge T})] = \mathbb{E}[u(\tau^h \wedge T, \Pi_{\bar{D}}(X_{\tau^h \wedge T}))] - u(0, x).
$$

Then, following a similar strategy as in the standard Cauchy problem, one has to apply Itô-Taylor formula and
use the regularity of the solution u to (1.32).

If the domain $D$ is sufficiently smooth and under non-degeneracy assumption (namely (UE) or (H)), one is able to prove that $\mathcal{E}_W^h$ is of order $h^{1/2}$. More precisely, Gobet [50], [51] proved that $\mathcal{E}_W^h = \mathcal{O}(h^{1/2})$ under some regularity assumption on $F$. Then, Gobet and Menozzi [54] sharpened this result by proving a lower bound of order $h^{1/2}$ and a first order expansion that writes $\mathcal{E}_W^h = C_1 h^{1/2} + o(h^{1/2})$.

1.2.4 Multi-level Monte Carlo simulation

The multi-level Monte Carlo method introduced by Heinrich [61] for parametric integration and recently developed by Giles [48] for Monte Carlo simulation as an extension of the two level method of Kebaier [70] allows to minimise the simulation cost of the quantity $\mathbb{E}[Y]$ where $Y$ is an $\mathbb{R}$-valued random variable that can only be approximated by a sequence $(Y^n)_{n \ge 1}$ of easily simulatable random variables as the *bias parameter* $n$ goes to infinity with a *weak discretisation error* or *bias* $\mathbb{E}[Y] - \mathbb{E}[Y^n]$ of order $n^{-\alpha}$, $\alpha > 0$.

This method is now very popular, notably in mathematical finance where as mentioned in the previous sections the
random variable of interest is given by a functional of an SDE with dynamics (1.17). However, it has been widely applied
to various problems of numerical probability, see Giles [47], Dereich [36], Giles, Higham and Mao [49] among others.
We refer the interested reader to the webpage: http://people.maths.ox.ac.uk/gilesm/mlmc_community.html for
further developments.
Let us be more specific. In this context, the standard Monte Carlo method uses the statistical estimator $M^{-1} \times \sum_{j=1}^{M} Y^{n,(j)}$ where the $(Y^{n,(j)})_{1 \le j \le M}$ are $M$ independent copies of $Y^n$. Given the order of the weak error, a natural question is to find the optimal choice of the sample size $M$ to achieve a global error. If the weak discretisation error is of order $n^{-\alpha}$ (for instance $\alpha = 1$ for the standard Cauchy problem and $\alpha = \frac{1}{2}$ for the Cauchy Dirichlet problem) then for a total error of order $n^{-\alpha}$ ($\alpha \in [1/2, 1]$), the minimal computation necessary for the standard Monte Carlo algorithm is obtained for $M = n^{2\alpha}$, see Duffie and Glynn [38]. So, if the computational cost required to simulate one sample of $Y^n$ is of order $n$ (as it is the case for the Euler scheme (1.19)) then the optimal computational cost of the Monte Carlo method is $C_{MC} = C \times n^{2\alpha+1}$, for a positive constant $C > 0$.

In order to reduce the complexity of the computation, the principle of the multi-level Monte Carlo method introduced by Giles [48] consists in using the telescopic sum

$$ \mathbb{E}[Y^{m^L}] = \mathbb{E}[Y^1] + \sum_{\ell=1}^{L} \mathbb{E}[Y^{m^\ell} - Y^{m^{\ell-1}}], $$

for an integer $m > 1$ satisfying $m^L = n^\alpha$. For each level $\ell \in \{1, \dots, L\}$ the numerical computation of $\mathbb{E}[Y^{m^\ell} - Y^{m^{\ell-1}}]$ is achieved by the standard Monte Carlo method using $N_\ell$ independent samples of $(Y^{m^{\ell-1}}, Y^{m^\ell})$. An important point is that the random samples $Y^{m^\ell}$ and $Y^{m^{\ell-1}}$ are perfectly correlated. Then the expectation $\mathbb{E}[Y^n]$ is approximated by the following multi-level estimator

$$ E_{MLMC} := \frac{1}{N_0} \sum_{j=1}^{N_0} Y^{1,(j)} + \sum_{\ell=1}^{L} \frac{1}{N_\ell} \sum_{j=1}^{N_\ell} (Y^{m^\ell,(j)} - Y^{m^{\ell-1),(j)}), $$

where for each level $\ell$, $(Y^{m^\ell,(j)})_{1 \le j \le N_\ell}$ is a sequence of i.i.d. random variables with the same law as $Y^{m^\ell}$. The variance of the new estimator is the sum of the local variance:

$$ \mathrm{Var}(E_{MLMC}) = \sum_{\ell=0}^{L} N_{\ell}^{-1} V_{\ell}, \quad V_{\ell} := \mathrm{Var}(Y^{m^{\ell}} - Y^{m^{\ell-1}}), \quad \ell = 0, \dots, L $$

with the convention $Y^{m^{-1}} = 0$. The computational cost is proportional to $\sum_{\ell=0}^{L} N_\ell m^{-\ell}$.

A simple optimisation problem shows that the variance of the multi-level estimator is minimised for a fixed computational cost by choosing $N_\ell$ proportional to $\sqrt{V_\ell m^{-\ell}}$. In the case we have in mind corresponding to the Euler discretisation scheme with Lipschitz-continuous coefficients, one has $Y^{m^\ell} = F(X_T^{m^\ell})$ for a Lipschitz continuous function $F$ so that $V_\ell = \mathcal{O}(m^{-\ell})$ and $N_\ell$ is proportional to the variance $V_\ell$ (at least asymptotically). This leads us to choose $N_\ell$ satisfying $N_\ell^{-1} V_\ell = \mathcal{O}(n^{-2\alpha}(L+1)^{-1})$. Then the variance of $E_{MLMC}$ is $\mathcal{O}(n^{-2\alpha})$ and the biais is still $n^{-\alpha}!$ Moreover, in order to achieve a global error of order $n^{-\alpha}$, the computational cost of the new estimator is $\mathcal{O}(n^{-2\alpha}L^2) = \mathcal{O}(n^{-2\alpha}(\log(n))^2)$. Based on this heuristic, we clearly see that the new multi-level Monte Carlo estimator outperforms the basic Monte Carlo estimator in terms of computational complexity.

More recently, Ben Alaya and Kebaier [1] proposed a different analysis to obtain the optimal choice of the parameters that relies on a Lindeberg-Feller CLT for the multi-level Monte Carlo algorithm.

More importantly, Lemaire and Pagès [93] proposed a multi-level Richardson-Romberg estimator which combines the higher order bias cancellation of the multi-step Richardson Romberg method investigated in [109] and the variance control of the multi-level Monte Carlo method briefly exposed above. Notably, in standard frameworks like the discretisation of diffusion processes, for a global error of order $n^{-\alpha}$, the computational cost is of order $\mathcal{O}(n^{-2\alpha} \log(n))$.

In Chapter 4, we extend the scope of the multi-level Monte Carlo method to the framework of stochastic optimisation by means of stochastic approximation algorithm. The optimal choice of the parameters is done by proving a CLT for the new estimators. We also present some numerical results that confirm the theoretical analysis and show a significant reduction in the initial computational cost.

### 1.2.5 Unbiased simulation of SDE

Unbiased simulation methods of SDEs offers an alternative to the Multi level Monte Carlo for the computation of $\mathbb{E}[F(X_T)]$ where $(X_t)_{t \in [0,T]}$ is solution to (1.17). We differentiate here between exact methods and unbiased methods. Exact methods allows to sample a path, at a finite set of points, with the exact distribution of the SDE, while unbiased means that we can estimate $\mathbb{E}[F(X_T)]$ without any bias, which is enough for many applications such that the pricing and hedging of financial claims.

An example of an exact simulation method is given in Beskos and al. [20] in dimension $d = q = 1$ based on the Lamperti transform of the original SDE. For some extensions to the SABR and the Heston models we refer the reader to Broadie and Kaya [26], Chen and al. [30]. However, all these methods cannot be easily extended to a general SDE.
The unbiased simulation technique for the multi-dimensional SDE has been introduced by Bally and Kohatsu-Higa [7] and developed by Anderson and Kohatsu-Higa [4] for multi-dimensional diffusions with hölder continuous coefficients. This method is based on a probabilistic interpretation of the parametrix method introduced by E. Lévi to construct the fundamental solution of elliptic and parabolic PDE. The proposed method is also based on the Euler scheme although in the proposed simulation scheme the time partition is random.

We now briefly expose the main guidelines for multi-dimensional diffusion processes with infinitesimal general $\mathcal{L} = \frac{1}{2} \sum_{i,j} a_{i,j} \partial_{i,j} + \sum_i b_i \partial_i$ where $a$ is assumed to be uniformly elliptic. The coefficients $a$ and $b$ are smooth, say $a_{i,j} \in C_b^2(\mathbb{R}^d)$ and $b_i \in C_b^1(\mathbb{R}^d)$. We know that the unique solution to (1.17) defines a strong Markov semigroup $(P_t)_{t \ge 0}$ defined for all bounded measurable maps by $P_t f(x) = \mathbb{E}[f(X_t)]$. Under these assumptions it is also known that $P_t f \in C_b^2(\mathbb{R}^d)$ for all $t > 0$, see e.g. Friedmann [44].

We introduce the approximation process $\bar{X}$ defined by $\bar{X}_t = x + b(x)t + \sigma(x)W_t$, that is the process obtained from (1.17) by freezing the drift and the diffusion coefficients at the initial point $x$. We introduce accordingly the collection of linear maps $(\bar{P}_t)_{t \ge 0}$ defined for all bounded measurable functions $f$ by $\bar{P}_t f(x) = \mathbb{E}[f(\bar{X}_t)] = \int f(y)g(a(x)t, y-b(x)t-x)$, where, for a symmetric positive definite matrix $a$, $g(a, x)$ is the probability density function of $N(0, a)$ taken at $x$. In particular, one has $\partial_t \bar{P}_t f(x) = \bar{P}_t \bar{\mathcal{L}} f$ for every $t > 0$. As a consequence, using an integration by parts formula the following one step formula holds

$$P_t f(x) = \bar{P}_t f(x) + \int_0^t \partial_s (\bar{P}_{t-s} P_s f(x)) ds = \bar{P}_t f(x) + \int_0^t \bar{P}_{t-s} (\mathcal{L} - \bar{\mathcal{L}}) P_s f(x) ds = \bar{P}_t f(x) + \int_0^t S_{t-s} P_s f(x) ds \quad (1.33)$$

where we introduced the operator $S_t f(x) = \int_{\mathbb{R}^d} f(y) \theta_t(x, y) g(a(x)t, y - b(x)t - x) dy$ with

$$\theta_t(x, y) = \frac{1}{2} \sum_{i,j=1}^{d} \kappa_t^{i,j}(x, y) - \sum_{i=1}^{d} \rho_t^i(x, y),$$

$$\kappa_t^{i,j}(x, y) = \hat{c}_{i,j}^2 a_{i,j}(y) + \hat{\partial}_j a_{i,j}(y) H^i(a(x)t, y - x - b(x)t) + \hat{\partial}_i a^{i,j}(y) H^j(a(x)t, y - x - b(x)t) + (a_{i,j}(y) - a_{i,j}(x)) H^{i,j}(a(x)t, y - x - b(x)t),$$

$$\rho_t^i(x, y) = \hat{\partial}_i b^i(y) + (b^i(y) - b^i(x)) H^i(a(x)t; y - x - b(x)t).$$

Here $H^i(a, x) = -(a^{-1}x)^i$ and $H^{i,j}(a, x) = (a^{-1}x)^i(a^{-1}x)^j - (a^{-1})^{i,j}$ are the Hermite polynomials of order 1 and 2 associated to the Gaussian density $g(a, x)$. From the space-time inequality $\forall x \in \mathbb{R}, |x|^{p_e q x^2} \le ((p/(2qe))^{p/2})$, valid for any $p, q > 0$, one easily gets $|S_{t-s} f|_{\infty} \le C(t-s)^{-\frac{1}{2}}|f|_{\infty}$, for a constant $C$ depending only on the coefficients $a, b$ and their derivatives.

Iterating the first step expansion (1.33), one proves that the following expansion holds

$$P_t f(x) = \sum_{n \ge 0} \int_{\Delta_n^*(t)} ds_n \prod_{j=0}^{n-1} S_{s_{j+1}-s_j} \bar{P}_{t-s_n} f(x), \quad (1.34)$$

where $\Delta_n^*(t) := \{s_n = (s_1, \dots, s_n) \in [0, t]^n : 0 < s_1 < \dots < s_n < t\}$ and $s_0 = 0$.

In order to provide a probabilistic representation of the above series, for a time partition $\pi : 0 = s_0 < s_1 < \dots < s_n < s_{n+1} = t$, we introduce the following discrete time process

$$X_{s_{i+1}}^{\pi} = X_{s_i}^{\pi} + b(X_{s_i}^{\pi})(s_{i+1} - s_i) + \sigma(X_{s_i}^{\pi})(W_{s_{i+1}} - W_{s_i}), \quad i = 0, \dots, n.$$

Then, the following intermediary probabilistic representation holds

$$P_t f(x) = \sum_{n \ge 0} \int_{\Delta_n(t)} ds_n \mathbb{E}[f(X_t^\pi) \prod_{i=0}^{n-1} \theta_{s_{j+1}-s_j}(X_{s_j}^\pi, X_{s_{j+1}}^\pi)], \quad (1.35)$$

and it remains to provide a probabilistic interpretation of the time integrals. To this end, we let $N(t)$ be a Poisson process with intensity parameter $\lambda > 0$ and define $N \equiv N(t)$. Let $\tau_1, \dots, \tau_N$ be the event times of the Poisson process and set $\tau_0 = 0$ and $\tau_{N+1} = t$. Since, conditional on $N$, the event times are distributed as a uniform statistic, $\mathbb{P}(N = n, \tau_1 \in ds_1, \dots, \tau_n \in ds_n) = \lambda^n e^{-\lambda t}$, for $s_n \in \Delta_n^*(t)$ and the following probabilistic representation holds

$$P_t f(x) = \mathbb{E}[f(X_t)] = e^{\lambda t} \mathbb{E}[f(X_t^\pi) \prod_{i=0}^{N-1} \theta_{\tau_{j+1}-\tau_j}(X_{\tau_j}^\pi, X_{\tau_{j+1}}^\pi)]. \quad (1.36)$$

Let us note that the right-hand side of (1.36) can be sampled exactly, this gives us an unbiased simulation method for the computation of $\mathbb{E}[f(X_T)]$. The simulation algorithm consists in independently simulating the event times
of a Poisson process and an Euler scheme with a time partition given by the simulated random times. However,
the method suffers from a poor convergence rate since the variance of the associated estimator may be infinite. An
importance sampling on the time steps using a Beta or Gamma distribution is proposed in [4] which allows to achieve
finite variance.

Another method is proposed in [7] and [4] which relies on a backward Euler scheme. The backward method can be applied for Hölder continuous coefficients. In this case, the Euler scheme is ran backward on the random time grid.

In Chapter 6, we will extend the above probabilistic representation to the one dimensional Cauchy-Dirichlet problem. More precisely, we obtain probabilistic representation of the process $(\tau \wedge t, X_{\tau \wedge t})_{t \ge 0}$ and obtain partial integration by parts formula that can be explicitly simulated without any bias.

We also mention the work of Henry-Labordère and al. [62] where the authors proposed a probabilistic representation for a class of semi-linear parabolic PDEs using a branching process. Usually this class of PDEs is solved using probabilistic numerical schemes of Backward SDEs. Here the algorithm is forward and backward regressions are no longer needed.
# Part I : Stochastic approximation algorithm and applications
# Chapter 2

## Some contributions to risk management in financial and energy markets

In this chapter we present the two research works [5] and [6].

In [6], which can be seen as a continuation of [3], we study the theoretical and computational aspects of risk minimisation in financial and energy market models operating in discrete time. The risk is defined by considering a class of convex risk measures defined on $L^p(\mathbb{P})$ in terms of shortfall risk [6]. Existence of an optimal strategy is proved under the absence of arbitrage and the non-degeneracy of the price process in a non-Markovian framework. We propose a Robbins-Monro algorithm to compute the shortfall risk of a given real-valued random variable. Then, we propose and analyse a Newton-Raphson optimisation algorithm to compute the optimal static risk minimising strategy. Finally, in a dynamic framework a quantisation based stochastic approximation algorithms to compute the corresponding risk measures as well as the optimal dynamic strategies is proposed.

In the research work [5], in collaboration with A. Sagna of University Evry Val d'Essonne, we study a quantisation based Newton-Raphson algorithm to compute the optimal importance sampling (IS) parameter. The approach proposed here can be seen as a robust and automatic deterministic counterpart of recursive IS (by translation of the mean) by means of stochastic approximation algorithm, as proposed by Lemaire and Pagès [92] which may require tuning of the step sequence and a good knowledge of the payoff function in practice.

### 2.1 Risk minimisation in financial or energy markets operating in discrete time

In [6], we investigate the risk minimisation problem of a portfolio loss with maturity $T$ described by an $\mathbb{R}$-valued random variable $L$ which contains an observable but non-tradable source of risk. More precisely we consider an energy market operating at discrete trading dates $t_0 = 0 < t_1 < \dots < t_N = T$. We have $d$ assets available for trade with price process $X = (X^1, \dots, X^d)$ and $X^i = (X^i_\ell)_{0 \le \ell \le M}$ for $i = 1, \dots, d$.

In our framework, the source of market incompleteness comes from the presence in $L$ of a state process $Z$ that is observable but not available for trade. Thus, it induces a source of risk that is not completely hedgeable. Typically, in the electricity market, the loss $L$ suffered by an energy company may be due to an abnormal annual electricity (or gas) consumption. This consumption depends on the temperature, which is an observable but non tradable source of risk. In this example the process $(Z_\ell)_{1 \le \ell \le M}$ can be considered as the temperature which may influence not only the loss but the assets available for trade, i.e. electricity prices of spot and forward contracts (which are in this example the only assets available for hedge). The probability space is equipped with a filtration $\mathcal{F} = (\mathcal{F}_\ell)_{0 \le \ell \le M}$. Intuitively, $\mathcal{F}_\ell$ represents the observable information at time $t_\ell$ by all investors, so that $\mathcal{F}_\ell = \sigma\{X_i, Z_i; 0 \le i \le \ell\}$.

The gains from a trading strategy $\theta$ with an initial investment of 0 are described by the discrete stochastic integral $\sum_{\ell=1}^N \theta_{\ell-1} \Delta X_\ell$, where we denote by $\Delta X_\ell$ the increments $X_\ell - X_{\ell-1}$ and $\theta = (\theta_\ell)_{0 \le \ell \le M} \in \mathcal{A}_\mathcal{F}$, where $\mathcal{A}_\mathcal{F} = \{\theta = (\theta_\ell)_{0 \le \ell \le M-1} \mid \theta_\ell \in L^0_{\mathbb{R}^d}(\mathcal{F}_\ell, \mathbb{P}), \ell = 0, \dots, M-1\}$ is the set of admissible strategies, $L^0_{\mathbb{R}^d}(\mathcal{F}_\ell, \mathbb{P})$ denoting the space of all $\mathcal{F}$-mesurable and $\mathbb{P}$-a.s. finite random variables with values in $\mathbb{R}^d$. In order to measure the risk, we consider a class of risk measures defined on $L^p(\mathbb{P})$ in terms of shortfall risk. It is a standard class of convex risk measure as mentioned by Föllmer and Schied [41] and Kaina and Rüschendorf [69].

To be more precise about this risk measure, let $\ell: \mathbb{R} \to \mathbb{R}_+$ be a non-decreasing convex loss function, not identically constant, hence satisfying $\lim_{x \to +\infty} \ell(x) = +\infty$. Assume that $\lim_{x \to -\infty} \ell(x) = \underline{l} < l(x_0)$. For some $x_0 \in \mathbb{R}$, we define the acceptance set

$$ \mathcal{A} := \{L \in L^p : \mathbb{E}[\ell(-L)] \le \ell(x_0)\}, $$
and the corresponding shortfall convex risk measure $\rho: L^p(\mathbb{P}) \to \mathbb{R} \cup \{+\infty\}$

$$ \rho(L) := \inf \{ \xi \in \mathbb{R} : L + \xi \in \mathcal{A} \} = \inf \{ \xi \in \mathbb{R} : \mathbb{E}[\ell(-L - \xi)] \le l(x_0) \}, \quad (2.1) $$

Under the assumption on the function $\ell$, the infimum in (2.1) exists. Assuming that $\forall M > 0$, $\mathbb{E}[l(-L + M)] < +\infty$, one has $\rho(L) \in \mathbb{R}$. Moreover, it is the lowest solution of the equation:

$$ \mathbb{E}[\ell(-L - \xi)] = l(x_0). \quad (2.2) $$

Classical examples are the entropic risk measure defined by setting $\ell(x) = e^{\lambda x}$, $x \in \mathbb{R}$, $\lambda > 0$, $x_0 \in \mathbb{R}$, $\rho(L) = \frac{1}{\lambda} \log \mathbb{E}[e^{-\lambda(L+x_0)}]$ and the $L^p$-shortfall risk measure, defined by setting $\ell(x) = x^p \mathbf{I}_{\{x \ge x_0\}}$, $x \in \mathbb{R}$, $p \ge 1$, $x_0 > 0$.

The basic problem for the investor with a given portfolio loss L is to minimise the risk

$$ \inf_{\theta \in \mathcal{A}_F} \rho \left( L - \sum_{n=1}^{N} \theta_{n-1} \cdot \Delta X_n \right) = \inf_{\theta \in \mathcal{A}_F} \inf \left\{ \xi \in \mathbb{R} : \mathbb{E} \left[ \ell(-L + \sum_{n=1}^{N} \theta_{n-1} \cdot \Delta X_n - \xi) \right] \le \ell(x_0) \right\}. \quad (2.3) $$

The main difference between the CVaR risk minimisation problem investigated in [3] and the shortfall risk minimisation problem (2.3) is that a shortfall risk measure does not write as an expectation contrary to the CVaR but appears as the level of a function which can be written as the expectation of a loss. Hence the arguments developed in [3] do not work.

Under the absence of arbitrage, we prove that the following two steps strategy works:

* the first step consists in solving for every $\xi \in \mathbb{R}$, the following stochastic control problem,

$$ \inf_{\theta \in \mathcal{A}_F} \mathbb{E} \left[ \ell(-L + \sum_{n=1}^{N} \theta_{n-1} \cdot \Delta X_n - \xi) \right] \quad (2.4) $$

By performing a non-trivial dynamic programming argument and using a measurable selection theorem, see e.g. p. 266 Lemma 5 in Jacod and Shirayev [67] or Dellacherie and Meyer [35] Theorem 82, p. 252, we prove the existence of optimal strategies for the risk minimisation problem (2.4) under the absence of arbitrage property (NA). $^1$The optimal strategy is computed by a backward induction

$$ \varphi_N(\xi, \theta_1, \dots, \theta_{N-1}) := \ell(-L + \sum_{n=1}^{N} \theta_{n-1} \cdot \Delta X_n - \xi), \quad \text{a.s.} $$

for all $\xi \in \mathbb{R}$, $\theta \in \mathcal{A}_F$; for $0 \le t < N$ and for all $\xi \in \mathbb{R}$, $(\theta_1, \dots, \theta_t) \in L_{\mathcal{F}_0}^0(\mathbb{P}) \times \dots \times L_{\mathcal{F}_{t-1}}^0(\mathbb{P})$,

$$ \varphi_t(\xi, \theta_1, \dots, \theta_t) := \inf_{\theta_{t+1} \in \Theta_t} \mathbb{E}[\varphi_{t+1}(\xi, \theta_1, \dots, \theta_t) | \mathcal{F}_t], \quad \text{a.s.} \quad (2.5) $$

* once an optimal solution to (2.4) is computed (if any), say $\theta^*(\xi)$, $\xi \in \mathbb{R}$, the second step consists in finding the lowest solution $\xi^*$ of the following equation

$$ \inf_{\theta \in \mathcal{A}_F} \mathbb{E} \left[ \ell(-L + \sum_{n=1}^{N} \theta_{n-1} \cdot \Delta X_n - \xi) \right] = \mathbb{E} \left[ \ell(-L + \sum_{n=1}^{N} \theta_{n-1}^*(\xi) \cdot \Delta X_n - \xi) \right] = \ell(x_0). \quad (2.6) $$

Then we are concerned by the problem of computing the shortfall risk (2.1) of a random variable X and the optimal strategy for the risk minimization problem (2.3). In order to estimate the shortfall risk, which is a solution to equation (2.2), we propose a Robbins-Monro stochastic algorithm. Assuming that there exists a positive function $\varphi$ such that $\mathbb{E}[\ell^2(-L - \xi)] \le C(1 + \varphi(\xi))$, $\xi \in \mathbb{R}$, we implement

$$ \xi_{n+1} = \xi_n + \gamma_{n+1} H(\xi_n, L^{(n+1)}), \quad n \ge 0, \quad H(\xi, x) := \frac{\ell(-x - \xi) - l(x_0)}{(1 + \varphi(\xi))^{\frac{1}{2}}} \quad (2.7) $$

where $(L^{(n)})_{n \ge 1}$ is an i.i.d. sequence of random variables with the same distribution as L, independent of $\xi_0 \in L^2(\mathbb{P})$ and $(\gamma_n)_{n \ge 1}$ is a step sequence satisfying the classical condition (1.9). We prove that the sequence $(\xi_n)_{n \ge 0}$ a.s. converges to a random variable $\rho(L)$ taking values in the set $\mathcal{T}^* := \{\xi : \mathbb{E}[\ell(-L - \xi)] = -\ell(x_0)\}$ as $n \to \infty$.

<footnote>Denoting $V_t^{v,\theta} = v + \sum_{\ell=1}^t \theta_\ell \cdot \Delta X_\ell$, one has (NA) : $\forall \theta \in \mathcal{A}_G$, $(V_M^{0,\theta} \ge 0 \text{ a.s.} \implies V_M^{0,\theta} = 0 \text{ a.s.})$</footnote>
Moreover, if one implements the recursive scheme (2.7) with a step sequence $\gamma_n = \frac{a}{(b+n)^\beta}$, $a, b > 0$, $\frac{1}{2} < \beta < 1$ and computes the empirical mean $\bar{\xi}_{n+1} = \frac{1}{n+1} \sum_{k=0}^{n} \xi_k$, $n \ge 0$ then according to the Ruppert and Polyak averaging principle one has

$$ \sqrt{n}(\bar{\xi}_n - \rho(X)) \stackrel{\mathcal{L}}{\Rightarrow} \mathcal{N}(0, \sigma_*^2), \quad \sigma_*^2 := \frac{\mathbb{E}\left[(\ell(-L - \rho(L))^2\right] - \ell^2(x_0)}{\mathbb{E}\left[\ell'(\ell(-L - \rho(L)))^2\right]}. $$

Similarly to CVaR hedging we first investigate the case of static hedging, that is, the holder of the portfolio only uses one step strategy. In this case, the problem (2.4) and (2.6) can be written

$$ \inf_{\theta \in \mathbb{R}^d} \mathbb{E}[\ell(-L + \theta \cdot X - \xi)] = \ell(x_0), \quad \text{with } X := X_N - X_0. \qquad (2.8) $$

which amounts to the computation of $(\xi^*, \theta^*) \in \mathbb{R} \times \mathbb{R}^d$ of the unique solution (for simplicity $\ell$ is assumed to be increasing and strictly convex) to the following systems of $d+1$ non-linear equations of $d+1$ unknowns

$$ \mathbb{E}[\ell(-H + \theta \cdot X - \xi)] = \ell(x_0), \quad \mathbb{E}[X \ell'(-H + \theta \cdot X - \xi)] = 0. \qquad (2.9) $$

In order to solve this system, we implement the following Newton-Raphson's optimisation algorithm

$$ z_{k+1} = z_k - Dh_n(z_k)^{-1}h_n(z_k), \quad k \ge 0, \quad z_0 := (\xi_0, \theta_0) \qquad (2.10) $$

where we approximated $h$ and $Dh$ using Monte Carlo simulation with $n$ samples

$$ h_n(\xi, \theta) = \left( \frac{1}{n} \sum_{k=1}^{n} \ell(-L^{(k)} + \theta \cdot X^{(k)} - \xi) - l(x_0), \frac{1}{n} \sum_{k=1}^{n} L^{(k)} \ell'(-L^{(k)} + \theta \cdot X^{(k)} - \xi) \right) $$

and

$$ Dh_n(\xi, \theta) = \left( \frac{-\frac{1}{n} \sum_{k=1}^{n} \ell'(-L^{(k)} + \theta \cdot X^{(k)} - \xi)}{-\frac{1}{n} \sum_{k=1}^{n} X^{(k)} \ell''(-L^{(k)} + \theta \cdot X^{(k)} - \xi)} \frac{\frac{1}{n} \sum_{k=1}^{n} X^{(k)} \ell'(-L^{(k)} + \theta \cdot X^{(k)} - \xi)}{\frac{1}{n} \sum_{k=1}^{n} X^{(k)} (X^{(k)})^T \ell''(-L^{(k)} + \theta \cdot X^{(k)} - \xi)} \right) $$

$(L^{(k)}, X^{(k)})_{1 \le k \le n}$ being an i.i.d. sequence of random vectors having the distribution of $(H, X)$. The recursive algorithm (2.10) is known to converge to $(\xi_n^*, \theta_n^*)$ unique solution to $h_n(\xi, \theta) = 0$ if $z_0 := (\xi_0, \theta_0)$ is sufficiently close to $(\xi_n^*, \theta_n^*)$. Moreover, we prove that $(\xi_n^*, \theta_n^*)_{n \ge 1}$ a.s. converges to $(\xi^*, \theta^*)$ as $n \to \infty$.

Finally, in the dynamic framework we propose a probabilistic algorithm to compute the optimal strategy using optimal Markovian quantization. As argued in Bally and Pagès [8] and Pagès and al. [111], Markovian quantization algorithms are suited for stochastic control problem in high dimensions. It consists in approximating the underlying Markov process by a process taking its values in optimized finite grids which take into account the fine structure of the Markov dynamics. We suppose that $\mathcal{F}_t := \sigma((X_s, Z_s), s \le t)$, $t = 1, \dots, N$ and restrict consideration to cases of strategies taking values in a bounded domain $A \subset \mathbb{R}^d$. Like in the CVaR hedging framework, the process $Z$ is observable but not available for trade. We also assume that the loss is given by $L = \phi(x)$, $x \in \mathbb{R}^{d+q}$. We introduce the notation $V_t = \sum_{n=1}^t \theta_{n-1} \Delta X_n$, $t \ge 1$ and $V_0 = -\xi$.

The dynamics of the discretized process $X := (X_t, Z_t)_{0 \le t \le N}$ is given by a Gaussian Euler scheme with step $h$ that is

$$ X_{t+1} = X_t + b(X_t)h + \sigma(X_t)h^{1/2}\varepsilon_{t+1} := G_h(X_t, \varepsilon_{t+1}), \quad t = 0, \dots, N-1, \qquad (2.11) $$

where $\varepsilon_{t+1} = h^{-1/2}(W_{t+h} - W_t) \sim \mathcal{N}(0, I_{d+q})$ independent of $\mathcal{F}_t$. Given a strategy $\theta \in \mathcal{A}_{\mathcal{F}}$, the dynamics of $(V_t)_{t \ge 0}$ is given by:

$$ V_{t+1} = V_t + \theta_t \cdot (X_{t+1} - X_t) = V_t + \theta_t \cdot B(X_{t+1} - X_t) := H_h(X_t, V_t, \theta_t, \varepsilon_{t+1}), \qquad (2.12) $$

where $B = (I_d, 0_{d \times q})$. The dynamic programming principle that we established can be solved recursively by the backward algorithm:

$$ \bar{\varphi}_N(x, v) = l(-\phi(x) + v), \quad (x, v) \in \mathbb{R}^{d+q} \times \mathbb{R}, \qquad (2.13) $$

$$ \bar{\varphi}_t(x, v) = \inf_{\theta \in A} \mathbb{E} \left[ \bar{\varphi}_{t+1}(X_{t+1}^x, V_{t+1}^{x,v,\theta}) \middle| (X_t, V_t) = (x, v) \right], (x, v) \in \mathbb{R}^{d+q} \times \mathbb{R} \qquad (2.14) $$

for $t = 0, \dots, N-1$, where $X_{t+1}^x := G_h(x, \varepsilon_{t+1})$ and $V_{t+1}^{x,v,\theta} := H_h(x, v, \theta, \varepsilon_{t+1})$. At step 0, one has to solve $\bar{\varphi}_0(X_0, -\xi) = l(x_0)$. In the spirit of the (multi-dimensional) Markov chain approximation method for solving non-linear problems (see e.g. [8], [111]), we consider an approximating dynamic of $(X_t, V_t)_{0 \le t \le N}$ by a process $(\tilde{X}_t, \tilde{V}_t)_{0 \le t \le N}$ taking its value in a finite grid.

Let $\Gamma_V := (2\delta)\mathbb{Z} \cap \bar{B}(0, R) = \{v \in \mathbb{R} : v = 2\delta z, \text{ for some } z \in \mathbb{Z}, |v| \le R\}$ be a fixed bounded lattice grid on $\mathbb{R}$. The projection on $\Gamma_V$ according to the closest neighbour rule is denoted $\pi_V$. At each time $t$, we consider a finite grid
<table><caption>Table 2.1: Shortfall risk estimation and One Step risk hedging of Spark Spread option</caption><thead><tr><th>l</th><th>No hedging</th><th>One step hedging</th><th>θ*</th></tr></thead><tbody><tr><td>l<sub>exp</sub></td><td>-18.60 ±0.02119</td><td>-21.71 ±0.01809</td><td>15.98</td></tr><tr><td>l<sub>2</sub></td><td>-2.79 ±0.00349</td><td>-5.83 ±0.01440</td><td>4.68</td></tr><tr><td>l<sub>3</sub></td><td>-2.05 ±0.00197</td><td>-3.69 ±0.01253</td><td>2.56</td></tr><tr><td>l<sub>4</sub></td><td>-1.73 ±0.00132</td><td>-2.71 ±0.00983</td><td>1.54</td></tr></tbody></table>

$\Gamma_t = \{x_t^1, \dots, x_t^{N_t}\}$ on the state space $\mathbb{R}^{d+q}$, for the uncontrolled process $X$. The projection on the grid $\Gamma_t$ will be denoted $\pi_t$. We consider a quantized Euler scheme $\hat{X}$ and a finite state space Markov chain $\hat{V}$ both defined by

$$ \hat{X}_{t+1} = \pi_{t+1}(G_h(\hat{X}_t, \epsilon_{t+1})), \quad \hat{V}_{t+1} = \pi_V(H_h(\hat{X}_t, \hat{V}_t, \theta_{t+1}, \hat{X}_{t+1}), \quad t=0, \dots, N-1 $$

and following [111], we approximate (2.13) and (2.14), by the following stochastic control problem

$$ \hat{\phi}_N(x, v) = l(-\phi(x) + v), \quad (x, v) \in \Gamma_N \times \Gamma_V, \qquad (2.15) $$

$$ \hat{\phi}_t(x, v) = \inf_{\theta \in A} \mathbb{E} \left[ \hat{\phi}_{t+1}(\hat{X}_{t+1}, \hat{V}_{t+1}^{\theta}) \mid (\hat{X}_t, \hat{V}_t) = (x, v) \right] = \inf_{\theta \in A} \mathbb{E} \left[ \hat{\phi}_{t+1}(\hat{X}_{t+1}^x, \hat{V}_{t+1}^{x,v,\theta}) \right], (x, v) \in \Gamma_t \times \Gamma_V \qquad (2.16) $$

where $\hat{X}_{t+1}^x = \pi_{t+1}(G_h(x, \epsilon_{t+1}))$ and $\hat{V}_{t+1}^{x,v,\theta} = \pi_V(H_h(x, v, \theta, \hat{X}_{t+1}^x)$, $t=0, \dots, N-1$.

Now, using the optimal dispatching of the grid sizes $(N_t)_{1 \le t \le N}$ for the Markovian quantization of the process $X$ and Theorem 3.1 in Pagès and al. [111], under some mild assumptions on the coefficients $b$ and $\sigma$, for all $v \in \Gamma_V$ one has

$$ |\bar{\phi}_0(v) - \hat{\phi}_0(v)| \le C_1 \frac{N^{2-\frac{1}{p}}}{(M/N)^{1/(d+q)}} + C_2 N \delta + C_3(\bar{p}) \frac{N}{R^{\bar{p}-1}}, \qquad (2.17) $$

where $M := N_1 + \dots + N_N$, $\hat{p} := \max(2, p_1+1)$, $\bar{p} > 1$, $C_1, C_2$ and $C_3(\bar{p})$ are positive constants. Therefore, we deduce the convergence of the value function approximated by optimal quantization toward the value function of the original control problem and the convergence of $\hat{\xi}^*$ toward $\xi^*$ as $1/M$, $N\delta$ and $N/R^{\bar{p}-1}$ go to zero.

As a simple numerical illustration of our approach we consider a portfolio composed of an exchange option between gas and electricity (called Spark Spread) of maturity $T = 1$ year. The payoof of such option is defined by $L = (S_T^e - h_R S_T^g - C)_+$. Since Electricity has very limited storage possibilities, the holder of this option hedges by trading only a gas forward contract of maturity $T$. We choose to model the price $(S_t^e)_{0 \le t \le T}$ of the electricity spot contract and the price $(F_{t,T}^g)_{0 \le t \le T}$ of the gas forward contract of maturity $T$ by two correlated geometric Brownian motions, the correlation factor between the two Brownian motions is $\rho = 0.5$. The parameters are: heat rate $h_R = 4$BTU/kWh (BTU: British Thermal Unit), the generation costs $C = 3$/MWh, the two volatilities $\sigma_g = 0.4$, $\sigma_e = 0.8$, the correlation factor $\rho = 0.5$ and the electricity and gas initial spot prices are $S_0^e = 40$/MWh, $S_0^g = 3$/MMBTU. In this section, we only consider the one-step hedging and we refer the reader to [6] for numerical results in the dynamic framework. Results about shortfall risk estimations and one step hedging are summarised in Table 2.1. The first column corresponds to the loss function. Four loss functions are considered: $l_{exp}(x) = e^{\lambda x}$, $\lambda = 1/50$ and $x_0 = 0$; $l_p(x) = x^p \mathbf{I}_{\{x \ge 0\}}$, $x_0 = 1$. The second column corresponds to the risk without any hedging using the Robbins-Monro algorithm (2.7) with its confidence interval at level 95%. The third column is the estimate of the shortfall risk with a one step risk minimization strategy and the fourth is the estimate of the optimal strategy $\theta^*$.

## 2.2 Quantisation based recursive importance sampling

In this section, we present [5] related to a quantisation based recursive importance sampling (IS) technique. We follow the idea of Pagès and Lemaire [92] who proposed a recursive IS stochastic gradient algorithm for the computation of $\mathbb{E}[F(X)]$ where $X: (\Omega, \mathcal{A}, \mathbb{P}) \to (E, |\cdot|_E)$ is a random vector taking values in a Banach space $E$ and $F: E \to \mathbb{R}$ is a Borel function such that $F(X) \in L^2(\mathbb{P})$. We present briefly the procedure proposed in [92]. For sake of simplicity we confine the presentation to the case $E = \mathbb{R}^d$. Note that one can also consider a functional space of paths of a process $X = (X_t)_{t \in [0,T]}$ solution to the SDE (1.17).

The IS is done by translating the distribution of the underlying random vector. The main idea of importance sampling by translation applied to the computation of $\mathbb{E}[F(X)]$ is to use the invariance by translation of the Lebesgue measure, for every $\theta \in \mathbb{R}^d$,

$$ \mathbb{E}[F(X)] = \mathbb{E} \left[ F(X + \theta) \frac{p(X + \theta)}{p(X)} \right], \qquad (2.18) $$
where $p: \mathbb{R}^d \to \mathbb{R}$ is probability density function of $X$. Among all these random variables with the same expectation we want to select the one with the lowest variance, i.e. the one with lowest quadratic norm

$$Q(\theta) := \mathbb{E} \left[ F^2(X + \theta) \frac{p^2(X + \theta)}{p^2(X)} \right] = \mathbb{E} \left[ F^2(X) \frac{p(X)}{p(X - \theta)} \right], \quad \theta \in \mathbb{R}^d. \qquad (2.19)$$

Assuming that $\theta \mapsto \log p(x - \theta)$ is log-concave, $x \in \mathbb{R}^d$, and $\lim_{|\theta| \to \infty} p(x - \theta) = 0$ then one shows that $Q$ is convex and goes to infinity at infinity so that $\arg\min Q$ is non-empty. Assume now that $p$ is differentiable on $\mathbb{R}^d$. Under additional growth and integrability conditions that we do not detail here, $Q$ is differentiable on $\mathbb{R}^d$ with a gradient given by

$$\nabla Q(\theta) := \mathbb{E} \left[ F(X - \theta)^2 \underbrace{\frac{p^2(X - \theta)}{p(X)p(X - 2\theta)}}_{W(\theta, X)} \frac{\nabla p(X - 2\theta)}{p(X - 2\theta)} \right]. \qquad (2.20)$$

Now if the growth of the function $F$ is controlled by $G$ (say $G$ grows exponentially), under some integrability assumption, setting

$$\widetilde{W}(\theta, X) := \frac{e^{-2\rho|\theta|^b}}{1 + G(-\theta)^{2c}} W(\theta, X), \quad K(\theta, X) := F(X - \theta)^2 \widetilde{W}(\theta, X), \qquad (2.21)$$

where $\rho$ and $c$ are some positive constants, one proves that the function $K$ satisfies the linear growth assumption (1.12) of the Robbins-Monro theorem, see Corollary 1.1.2.

For instance, when $X \sim \mathcal{N}(0, I_d)$, taking $\rho = \frac{1}{2}$, $b = 2$, one gets

$$K(\theta, x) = \frac{1}{1 + G(-\theta)^{2c}} F^2(x - \theta)(2\theta - x).$$

Coming back to the general case, the mean function of $K$ has the same zeroes as $\nabla Q$, namely one has

$$\{\theta \in \mathbb{R}^d | \mathbb{E}[K(\theta, X)] = 0\} = \{\theta \in \mathbb{R}^d | \nabla Q(\theta) = 0\}.$$

Since $Q$ is convex $h := \nabla Q$ satisfies (1.11) and the same holds for the function $h: \theta \mapsto h(\theta) = \mathbb{E}[K(\theta, X)]$. Hence, according to the Robbins-Monro algorithm the sequence $(\theta_n)_{n \ge 0}$ defined by

$$\theta_{n+1} = \theta_n - \gamma_{n+1} K(\theta_n, X^{(n+1)}), \quad \theta_0 \in \mathbb{R}^d, \qquad (2.22)$$

a.s. converges to an argmin $Q$-valued (square integrable) random variable $\theta^*$.

In the infinite dimensional setting when $E = C([0, T], \mathbb{R}^d)$ (equipped with the supremum norm) in order to deal with the case where $X$ is a continuous path-dependent diffusion process, a recursive scheme similar to (2.22) is proposed in [92]. However, one relies on the Girsanov theorem instead of the invariance by translation of the Lebesgue measure in order to derive a representation similar to (2.18).

From a numerical point of view, the procedure (2.22) suffers from several drawbacks. First the tuning of the algorithm needs a good knowledge of the behaviour of $F$ at infinity. Moreover, in practical implementations, recursive importance sampling methods using stochastic approximation algorithm need specific tuning of the step sequence and the convergence may be slow even when the dimension of the underlying random vector $X$ is low. As a consequence, sometimes one may prefer deterministic optimisation procedures when it is possible.

In [5], we propose a deterministic version of the IS algorithm investigated in [92]. The main idea can be summed up as follows. Assume that $X$ with probability density function $p$, supposed to be two times continuously differentiable, can be approximated by a discrete random variable taking $N$ values and whose distribution is known. Let us denote by $(x_k)_{1 \le k \le N}$ its possible values and by $(p_k)_{1 \le k \le N}$ the associated probabilities. For every $\theta \in \mathbb{R}^d$, the function $Q$ given by (2.19) as well as $\nabla Q$ and $\nabla^2 Q$ can be easily approximated by

$$Q_N(\theta) = \sum_{k=1}^{N} p_k F^2(x_k) \frac{p(x_k)}{p(x_k - \theta)}, \quad \nabla Q_N(\theta) = \sum_{k=1}^{N} p_k F^2(x_k) \frac{p(x_k)}{p^2(x_k - \theta)} \nabla p(x_k - \theta), \\ \nabla^2 Q_N(\theta) = \sum_{k=1}^{N} p_k F^2(x_k) \frac{p(x_k)}{p(x_k - \theta)} \left( 2 \frac{\nabla p(x_k - \theta) \nabla p(x_k - \theta)^t}{p^2(x_k - \theta)} - \frac{\nabla^2 p(x_k - \theta)}{p(x_k - \theta)} \right). \qquad (2.23)$$

Hence, one may devise a classical Newton-Raphson algorithm, namely

$$\theta_{k+1}^N = \theta_k^N - \nabla^2 Q_N(\theta_k^N)^{-1} \nabla Q_N(\theta_k^N), \quad \theta_0 \in \mathbb{R}^d, \quad k \ge 0 \qquad (2.24)$$
to approximate the optimal (variance) minimiser $\theta^*$ of $Q$. Starting from that idea, Jourdain and Lelong [68] proposes to replace $X$, where $X$ is assumed to be a $d$-dimensional standard normal random vector, by a close discrete random variable obtained by a Monte Carlo simulation. The (unique) minimum of $Q$ is approximated by the (unique) zero $\theta^N$ of $\nabla Q_N$ which can be computed using the deterministic Newton-Raphson algorithm (2.24) involving $\nabla Q_N$ and $\nabla^2 Q_N$. Moreover, several asymptotic properties are addressed.

Here, we propose to replace the Monte Carlo phase by a vector or functional quantization phase of $X$. Optimal quadratic quantisation of the random vector $X$ consists basically of finding the best approximation (in $L^2(\mathbb{P})$) of $X$ by a discrete random variable $\hat{X}^N$ (the quantization of $X$) taking at most $N$ values. When $X$ takes its values in $\mathbb{R}^d$, one speaks about vector quantization. When $X$ takes its values in an infinite dimensional Hilbert space like $L^2([0, T], dt)$ endowed with the usual norm $|u|_2 = (\int_0^T u(t)^2 dt)^{\frac{1}{2}}$, one speaks about functional quantization.

For $N \ge 1$ and $x := (x_1, \dots, x_N) \in E^N$ ($E$ will either be $\mathbb{R}^d$ or $L^2([0, T], dt)$ in what follows) an $N$-tuple referred as an $N$-quantizer, we let $\text{Proj}_x : E \to \{x_1, \dots, x_N\}$ be a Borel projection following the nearest neighbour rule. Then, the Borel partition $C = \{C_1, \dots, C_N\}$ of $E$ defined by $C_k := \text{Proj}_x^{-1}(x_k)$, $k = 1, \dots, N$ and satisfying

$$ \text{Proj}_x^{-1}(\{x_k\}) \subset \{y \in E, |x_k - y|_E = \min_{1 \le j \le N} |x_j - y|_E\}, \quad 1 \le k \le N, $$

is called a Voronoi tessellation of $E$ induced by $x$ and one defines accordingly the Voronoi quantization of $X$ induced by $x$ as follows

$$ \hat{X} = \text{Proj}_x (X) = \sum_{k=1}^{N} x_k \mathbf{I}_{\{X \in C_k\}}. $$

Note that the discrete random variable $\hat{X}$ is the best $L^p(\mathbb{P})$-approximation of $X$ among all measurable random variable taking values in $x := \{x_1, \dots, x_N\}$. For any fixed $N$-quantizer $x$, we associate the $L^p(\mathbb{P})$-mean error $\mathbb{E}[|X - \hat{X}|^p]^{1/p}$ induced by $x$ and computing an optimal quantization of $X$ consists in finding an $N$-tuple $x \in E^N$ which minimises the $L^p$-mean error over $E^N$. It amounts to minimising the function $\Lambda_N^X : (x_1, \dots, x_N) \mapsto \mathbb{E}[|X - \hat{X}|^p]^{1/p} = \mathbb{E}[(\min_{1 \le i \le N} |X - x_i|_E)^p]^{1/p}$.

One shows that, if $X \in L_E^p(\mathbb{P})$, the function $Q_N^X$ reaches its minimum (at least) at one $N$-tuple $x^*$ called an optimal $N$-quantizer. This minimum is in general not unique, except in some cases, e.g. when $d = 1$ and the density of $X$ is log-concave. Moreover, the $L^p$-mean quantization error $e_{N,p}^X := \min_{E^N} Q_N^X$ converges toward 0 as $N \uparrow \infty$ and for “non-singular” $\mathbb{R}^d$-valued random vectors, the rate of convergence of $e_{N,p}^X$ is ruled by the so-called Zador Theorem (see Graf and Luschgy [57]) which asserts that $\lim_{N \to +\infty} N^{1/d} e_{N,p}^X = q_p(d) (\int f(x)^{d/(d+p)} dx)^{(d+p)/dp}$ if $X \in L_{\mathbb{R}^d}^{p+\delta}(\mathbb{P})$ where $f$ denotes the density of the absolutely continuous part of $\mathbb{P}_X$ and $q_p(d)$ a positive constant depending on $p$ and $d$.

The main idea now is that we know that $\hat{X}$ is close to $X$ in distribution and if one has a numerical access to the optimal $N$-quantizer $x$ with the associated weights sequence $(p_k = \mathbb{P}(X \in C_k(x)))_{1 \le k \le N}$ of the quantisation $\hat{X}$ then the computation of the functions $Q_N$, $\nabla Q_N$ and $\nabla^2 Q_N$ are dramatically faster compared to any other simulation based approach. For instance, a sharply optimised database of quantizers of standard univariate and multivariate Gaussian distributions is available on the web site www.quantize.maths-fi.com [112] for download. To an $N$-quantizer $x$ we associate the unique minimum $\hat{\theta}^N$ of $Q_N$, that is, the unique solution to the equation $\nabla Q_N(\theta) = 0$. One expects that $\hat{\theta}^N \to \theta^*$ as the size $N$ of its associated $N$-quantizer $x$ goes to infinity.

Under the assumptions that the payoff function $F$ is Lipschitz-continuous, the density $p$ is strictly log-concave with $\lim_{|x| \to +\infty} p(x) = 0$ then under some growth condition on the first and second derivatives of $p$ that we do not detail here, one gets

$$ \hat{\theta}^N \to \theta^*, \quad Q_N(\hat{\theta}^N) \to Q(\theta^*) \text{ as } N \to +\infty $$

as soon as $\mathbb{E}[|X - \hat{X}^N|^2] \to 0$ as $N \to +\infty$.

In order to obtain a convergence rate, we assume that $\nabla p/p$ is $\eta$-Hölder. For instance, for the multi-dimensional Gaussian law $\eta = 1$. Under these additional assumption, one obtains the following convergence rate $|\hat{\theta}^N - \theta^*| = \mathcal{O}(\|X - \hat{X}^N\|_2^\alpha)$ as $N \to \infty$. For instance, when $X$ is a $d$-dimensional Gaussian vector, one has $|\hat{\theta}^N - \theta^*| \le CN^{-1/d}$.

We also extend our approach to the infinite dimensional setting, i.e. the case of path-dependent diffusion. Instead of the invariance by translation of the Lebesgue measure, we rely on the Girsanov transform to play the role of mean translator. To be more precise, we consider a $d$-dimensional Itô process $X$ solution to the SDE

$$ X_t = x + \int_0^t b(s, X^s) ds + \int_0^t \sigma(s, X^s) dW_s \quad (2.25) $$

where $W = (W_t)_{t \in [0,T]}$ is a $q$-dimensional Brownian motion and where $X^t := (X_{t \wedge s})_{s \in [0,T]}$ (with $X = X^T$) is the stopped process at time $t$. We assume that a strong solution to the SDE (2.25) exists. In order to devise a quantized
IS procedure for $\mathbb{E}[F(X)]$, $F$ being a Borel map defined on $C([0, T], \mathbb{R}^d)$ we consider a translation functional given by $\theta \in L_{T,q}^2 := L^2([0, T], \mathbb{R}^q)$ which is slightly less general than the one used in [92]. From the Girsanov Theorem it follows that for every $\theta \in L_{T,q}^2$

$$ \mathbb{E}[F(X)] = \mathbb{E}\left[F(X^{(\theta)})e^{-\int_{0}^{T}\langle\theta_s, dW_s\rangle - \frac{1}{2}\|\theta\|_{L_{T,q}^2}^2}\right] $$

where $X^{(\theta)}$ denotes the unique strong solution to (2.25) with drift $b + \sigma\theta$. Among all these estimators we want to select the one with the lowest quadratic norm so that we want to solve the following minimization problem

$$ \min_{\theta \in L_{T,q}^2} Q(\theta) \quad \text{where} \quad Q(\theta) := \mathbb{E}\left[F^2(X^{(\theta)})e^{-2\int_0^T \langle\theta_s, dW_s\rangle - \|\theta\|_{L_{T,q}^2}^2}\right] = \mathbb{E}\left[F^2(X)e^{-\int_0^T \langle\theta_s, dW_s\rangle + \frac{1}{2}\|\theta\|_{L_{T,q}^2}^2}\right]. $$

Note that $Q$ is strictly convex and $\lim_{\|\theta\|_{L_{T,q}^2} \to \infty} Q(\theta) = +\infty$ so that it admits a unique minimum. Moreover, $Q$ is twice differentiable at any $\theta \in L_{T,q}^2$ and $\langle DQ(\theta), \psi \rangle_{L_{T,q}^2}$, $(D^2Q(\theta)\psi, \phi)$, for $\psi, \phi \in L_{T,q}^2$ admit a representation as an expectation.

In view of a numerical implementation of Newton-Raphson's algorithm to estimate a minimum of $Q$, we are led to consider a (non-trivial) finite dimensional subspace $E$ of $L_{T,p}^2$ spanned by an orthonormal basis $(e_1, \dots, e_m)$. Like for the finite dimensional framework, our procedure will be based on the representation (as an expectation) of the first and second differential of $Q$ on $E$ combined with functional quantization of the solution to (2.25).

The functional quantization of (2.25) is based on the Karhunen-Loève (K-L) decomposition $W \stackrel{L_T^2}{=} \sum_{n \ge 1} \sqrt{\lambda_n} \xi_n e_n$ where $(e_n)_{n \ge 1}$ is the K-L orthonormal basis, $\lambda_n$ are the corresponding eigenvalues of the covariance operator associated to the one-dimensional Brownian motion $W$ and $\{\xi_n := \langle W, e_n \rangle / \sqrt{\lambda_n}, n \ge 1\}$, is an i.i.d. sequence of random variables with standard normal distribution. One considers an optimal $N_n$-quantization $(N_n \ge 1)$ of $\xi_n$, denoted $\hat{\xi}_n$ where $\hat{\xi}_n := \text{Proj}_{x_n}(\xi_n)$, $x_n := (x_1^{N_n}, \dots, x_{N_n}^{N_n})$ is the unique optimal $N_n$-quantizer of the normal distribution and $N_1 \times \dots \times N_n \le N$, with $N_1, \dots, N_n \ge 1$. For large enough $n$, we set $N_n = 1$, $\hat{\xi}_n = 0$ (which is the optimal 1-quantization of $\mathcal{N}(0, 1)$) and we define the product quantizer by (the finite sum) $\hat{W}_t = \sum_{n \ge 1} \sqrt{\lambda_n} \hat{\xi}_n e_n(t)$. The product quantizer $\chi$ that produces the above Voronoi quantization $\hat{W}$ is defined by $\chi_i(t) = \sum_{n \ge 1} \sqrt{\lambda_n} x_{i_n} e_n(t)$, $i = (i_1, \dots, i_n, \dots) \in \prod_{n \ge 1} \{1, \dots, N_n\}$ in the sense that $\mathbb{P}(\hat{W} = \chi_i) = \prod_{n \ge 1} \mathbb{P}(\xi_n \in C_{i_n}(x_n))$ with an approximation rate $\|W - \hat{W}^N\|_2 \le C_T(\log N)^{-\frac{1}{2}}$. Then, for each $\underline{i}$, one solves the following ODEs

$$ dx_i(t) = \left( b(x_i(t)) - \frac{1}{2}\sigma\sigma'(x_i(t)) \right) dt + \sigma(x_i(t))d\chi_i(t), \quad (2.26) $$

and set $\tilde{X}_t^N = \sum_{i} x_i(t) \mathbf{I}_{\{C_{\underline{i}}(\chi)\}}(W)$, $N \ge 1$. Assuming that $b, \sigma$ are continuously differentiable, this non-Voronoi $N$-functional quantization of (2.25) produces an error rate $\|X - \tilde{X}^N\|_2 \le C_T(\log N)^{-\frac{1}{2}}$. We approximate $Q$ by $\tilde{Q}_N$, the $m$-tuple $\langle DQ(\theta), e_i \rangle_{L_{T,q}^2}$ by $J(\theta)_i = \langle D\tilde{Q}_N(\theta), e_i \rangle_{L_{T,q}^2}$, $1 \le i \le m$ and the $m \times m$ symmetric positive definite matrix $(D^2Q(\theta)e_i, e_j)$ by $H(\theta)_{i,j} := (D^2\tilde{Q}_N(\theta)e_i, e_j)$, $1 \le i \le m$, $1 \le j \le m$, for every $\theta \in E$ where we replaced the random variables $(X, W)$ by their quantized version $(\tilde{X}^N, \hat{W}^N)$ in the representation as an expectation of $Q$, $DQ$ and $D^2Q$. Hence, we compute the minimum $\hat{\theta}^N$ of $\tilde{Q}_N$ by devising the following Newton-Raphson algorithm

$$ \hat{\theta}_{k+1}^N = \hat{\theta}_k^N - H(\hat{\theta}_k^N)^{-1} J(\hat{\theta}_k^N), \quad k \ge 0, \quad \hat{\theta}_0^N \in E \text{ given.} \quad (2.27) $$

We briefly illustrate the performance of the algorithm in the infinite dimensional setting by considering three different basis of $L^2([0, 1], \mathbb{R})$

* a polynomial basis composed of the shifted Legendre polynomials $(\tilde{P}_n)_{n \ge 0}$ defined by

$$ \forall n \ge 0, \forall t \in [0, 1], \quad \tilde{P}_n(t) = P_n(2t - 1) \quad \text{where} \quad P_n(t) = \frac{1}{2^n n!} \frac{d^n}{dt^n}((t^2 - 1)). \quad (\text{ShLeg}) $$

* the Karhunen-Loève basis defined by

$$ \forall n \ge 0, \forall t \in [0, 1], \quad e_n(t) = \sqrt{2} \sin\left(\left(n + \frac{1}{2}\right)\pi t\right). \quad (\text{KL}) $$
• the Haar basis which is defined by

$$
\forall n \ge 0, \forall k = 0, \dots, 2^n - 1, \forall t \in [0, 1], \psi_{n,k}(t) = 2^{\frac{k}{2}} \psi(2^k t - n), \quad \psi(t) = \begin{cases} 1 & \text{if } t \in [0, \frac{1}{2}) \\ -1 & \text{if } t \in (\frac{1}{2}, 1) \\ 0 & \text{otherwise.} \end{cases} \quad (\text{Haar})
$$

We consider a Down&In Call option of strike K and barrier L. This option is activated when the underlying process X moves down and hits the barrier L. The payoff function at maturity T is defined by

$$
F(X) = (X_T - K)_+ \mathbf{1}_{\{\min_0 \le t \le T, X_t \le L\}}.
$$

A standard approach to price the option is to consider the continuous Euler scheme $\bar{X}$ of step $t_k = k \frac{T}{M}$ obtained by extrapolation of the Brownian Motion between two instants of discretization. For every $t \in [t_k, t_{k+1}]$, we can write

$$
\bar{X}_t = \bar{X}_{t_k} + b(\bar{X}_{t_k})(t - t_k) + \sigma(\bar{X}_{t_k})(W_t - W_{t_k}), \quad \bar{X}_0 = x_0 \in \mathbb{R}.
$$

By pre-conditioning,

$$
\mathbb{E}[(\bar{X}_T - K)_+ \mathbf{1}_{\{\min_{0 \le t \le T} \bar{X}_t \le L\}}] = \mathbb{E}\left[(\bar{X}_T - K)_+ \left(1 - \prod_{p=0}^{M-1} p(\bar{X}_{t_p}, \bar{X}_{t_{p+1}})\right)\right], \quad (2.28)
$$

where $p(x_p, x_{p+1}) = \mathbb{P}(\min_{t_p \le t \le t_{p+1}} \bar{X}_t \ge L | (\bar{X}_{t_p}, \bar{X}_{t_{p+1}}) = (x_p, x_{p+1}))$ is the probability of non-exit from $[L, +\infty)$ of some Brownian bridge. Using the law of the Brownian bridge (see e.g. Gobet [52]), we can write

$$
p(x_p, x_{p+1}) = 1 - \mathbb{P} \left( \min_{t \in [0, t_1]} W_t \le \frac{L - x_p}{\sigma(x_p)} \mid W_{t_1} = \frac{x_{p+1} - x_p}{\sigma(x_p)} \right) = \begin{cases} 1 - e^{-\frac{2(L-x_p)(L-x_{p+1})}{(t_{p+1}-t_p)\sigma^2(x_p)}}, & L \le \min(x_p, x_{p+1}), \\ 0, & \text{otherwise.} \end{cases} \quad (2.29)
$$

Hence we run our algorithm with the probabilistic representation appearing in the right-hand side of (2.28). We consider the local volatility model $b(x) = rx$, $\sigma(x) = \sigma x^{1+\beta}/(1+x^2)^{\frac{1}{2}}$, with $r = 4\%$, $\sigma = 5$, $x = 100$, $\beta = 0.5$. The number of Monte Carlo simulations $n$ is $5.10^4$ in every case. We set the optimal product quantizer at level $d_N = 966$ which corresponds to the optimal decomposition $N_1 = 23$, $N_2 = 7$, $N_3 = 3$, $N_4 = 2$ for the product quantization $(\hat{\xi}_1, \hat{\xi}_2, \hat{\xi}_3, \hat{\xi}_4)$. The numerical results are summarised in Table 2.2. In Figure 2.1, the optimal variance reducers for the local volatility model are depicted

<table><caption>Table 2.2: Down&In Call option in Local volatility model with $X_0 = 100$, $K = 115$, $L = 65$, $T = 1$, $\sigma = 5$, $r = 4\%$, $\beta = 0.5$, $d_N = 966$, $n = 5.10^4$, $M = 100$.</caption>
  <thead>
    <tr>
      <th>Basis</th>
      <th>m</th>
      <th>Price MC</th>
      <th>Variance MC</th>
      <th>Price QIS</th>
      <th>Variance QIS</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Constant</td>
      <td>1</td>
      <td>0.684</td>
      <td>25.80</td>
      <td>0.673</td>
      <td>15.33</td>
    </tr>
    <tr>
      <td>Legendre</td>
      <td>2</td>
      <td>0.711</td>
      <td>27.92</td>
      <td>0.662</td>
      <td>4.58</td>
    </tr>
    <tr>
      <td>(ShLeg)</td>
      <td>4</td>
      <td>0.683</td>
      <td>25.66</td>
      <td>0.684</td>
      <td>3.46</td>
    </tr>
    <tr>
      <td></td>
      <td>8</td>
      <td>0.686</td>
      <td>26.59</td>
      <td>0.685</td>
      <td>3.35</td>
    </tr>
    <tr>
      <td>Karhunen-Loève</td>
      <td>2</td>
      <td>0.680</td>
      <td>25.38</td>
      <td>0.696</td>
      <td>5.22</td>
    </tr>
    <tr>
      <td>(KL)</td>
      <td>4</td>
      <td>0.702</td>
      <td>26.39</td>
      <td>0.683</td>
      <td>6.53</td>
    </tr>
    <tr>
      <td></td>
      <td>8</td>
      <td>0.687</td>
      <td>26.39</td>
      <td>0.688</td>
      <td>5.80</td>
    </tr>
    <tr>
      <td>Haar</td>
      <td>2</td>
      <td>0.648</td>
      <td>24.90</td>
      <td>0.673</td>
      <td>8.15</td>
    </tr>
    <tr>
      <td>(Haar)</td>
      <td>4</td>
      <td>0.671</td>
      <td>25.15</td>
      <td>0.692</td>
      <td>5.46</td>
    </tr>
    <tr>
      <td></td>
      <td>8</td>
      <td>0.709</td>
      <td>30.17</td>
      <td>0.700</td>
      <td>5.29</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_2.png" /><figcaption>Figure 2.1: Down & In Call option: Optimal $\theta^N$ obtained by our algorithm in the case of the Local Volatility model for different basis and several values of $m$ ($m = 2$ for the left upper curves, $m = 4$ for the right upper curves and $m = 8$ for the lower curves).</figcaption></figure>
# Chapter 3

## Concentration inequalities and deviation estimates for stochastic approximation schemes

In this chapter, we present [7] which is a joint work with S. Menozzi from the university of Evry Val d'Esonne and [8] which is a joint work with M. Fathi from the university Paul Sabatier of Toulouse. Both works are related to concentration inequalities and deviation estimates for stochastic approximation schemes. Here we are interested in two kinds of discrete stochastic approximation schemes. The first one refers to the law of an Euler like discretisation scheme of a diffusion process at a fixed deterministic date and the second one concerns the law of a stochastic approximation algorithm at a given time-step.

Concentration inequalities consists in establishing a *non-asymptotic* Gaussian [7] or sub-Gaussian [8] control of the probability that the law of the scheme deviates from its target. More precisely, we consider the deviations between the expectation of a given function of an Euler like discretisation scheme of some diffusion process at a fixed deterministic time and its empirical mean obtained by the Monte Carlo procedure. We then consider some estimates concerning the deviation between the value at a given time-step of a stochastic approximation algorithm and its target.

Much work on non asymptotic deviation estimates for the law of the Euler discretisation scheme has now been conducted. In the ergodic framework and for a constant diffusion coefficient Gaussian controls have been obtained by Malrieu and Talay [94]. For the standard Euler scheme, a first attempt to establish two-sided Gaussian bounds for the statistical error can be found in Lemaire and Menozzi [91] under some non-degeneracy conditions and up to a systematic bias. Very recently, Honoré, Menozzi and Pagès [63] obtained non-asymptotic Gaussian concentration bounds for the difference between the invariant measure $\nu$ of an ergodic Brownian diffusion process and the empirical distribution of an approximating scheme with decreasing time step along a suitable class of smooth enough test functions $f$ such that $f - \nu(f)$ is a *coboundary* of the infinitesimal generator of the process.

Such concentration inequalities are strongly connected to transport-entropy inequalities. In Section 3.1, we present this connection and recall basic definitions and properties. For a complete overview and recent developments in the theory of transport inequalities, the reader may refer to the recent survey of Gozlan and Leonard [56] and the book of Villani [126]. We establish Gaussian concentration inequalities in [7] for the two specific problems mentioned above under suitable assumptions. The key tool consists in exploiting accurately the concentration properties of the increments of the schemes and proving that this concentration transfers to the law of the scheme. We will explain this point in Section 3.2. Finally, in Section 3.3 we explain the results obtained in [8] which improve and complete those obtained in [7]. The key point is to properly quantify the contribution of the diffusion term to the concentration regime. We also derive a non-asymptotic bound for stochastic approximation with averaging of trajectories, in particular we prove that averaging a stochastic approximation algorithm with a slow decreasing step sequence gives rise to optimal concentration rate.

### 3.1 A brief presentation of transport inequalities and concentration inequalities

In this section, we present the connection between transport inequalities and concentration inequalities. We also recall some basic definitions and properties that will be needed in the two next sections. We denote by $P(\mathbb{R}^d)$ the set of probability measures on $\mathbb{R}^d$.

For $p \ge 1$, we consider the set $P_p(\mathbb{R}^d)$ of probability measures with finite moment of order $p$. The Wasserstein
metric $W_p(\mu, \nu)$ of order $p$ between two probability measures $\mu, \nu \in \mathcal{P}_p(\mathbb{R}^d)$ is defined by

$$W_p^P(\mu, \nu) = \inf \left\{ \int_{\mathbb{R}^d \times \mathbb{R}^d} |x - y|^p \pi(dx, dy) : \pi \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d), \pi_0 = \mu, \pi_1 = \nu \right\}$$

where $\pi_0$ and $\pi_1$ are two probability measures standing for the first and second marginals of $\pi \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d)$. For $\mu \in \mathcal{P}(\mathbb{R}^d)$, we define the relative entropy w.r.t $\nu \in \mathcal{P}(\mathbb{R}^d)$ as

$$H(\mu, \nu) = \int_{\mathbb{R}^d} \log \left( \frac{d\mu}{d\nu} \right) d\mu$$

if $\mu \ll \nu$ and $H(\mu, \nu) = +\infty$ otherwise. We are now in position to define the notion of transport-entropy inequality. Here as below, $\Phi: \mathbb{R}_+ \to \mathbb{R}_+$ is a convex, increasing function with $\Phi(0) = 0$.

**Definition 3.1.1.** A probability measure $\mu$ on $\mathbb{R}^d$ satisfies a transport-entropy inequality with function $\Phi$ if for all $\nu \in \mathcal{P}(\mathbb{R}^d)$, one has

$$\Phi(W_1(\nu, \mu)) \le H(\nu, \mu).$$

For the sake of simplicity, we will write that $\mu$ satisfies $T_\Phi$.

The following proposition comes from Corollary 3.4. in Gozlan and Leonard [56].

**Proposition 3.1.1.** The following propositions are equivalent:

* The probability measure $\mu$ satisfies $T_\Phi$.

* For all 1-Lipschitz function $f: \mathbb{R}^d \to \mathbb{R}$, one has

$$\forall \lambda \ge 0, \quad \int \exp(\lambda f) d\mu \le \exp \left( \lambda \int f d\mu + \Phi^*(\lambda) \right),$$

where $\Phi^*$ is the monotone conjugate of $\Phi$ defined on $\mathbb{R}_+$ as $\Phi^*(\lambda) = \sup_{\rho \ge 0} \{\lambda \rho - \Phi(\rho)\}$. Note that since $\Phi$ is a non-negative convex function which is non-decreasing with $\Phi(0) < \infty$, the monotone conjugate of $\Phi^*$ is $\Phi$.

As an important example we will say that $\mu$ satisfies the *Gaussian concentration property* if there exists $\beta > 0$ such that $\Phi(\lambda) = \lambda^2/\beta$, $\lambda \ge 0$ so that $\Phi^*(\lambda) = \lambda^2\beta/4$. In this chapter we will simply say that $\mu$ satisfies $GC(\beta)$.

From the Markov exponential inequality and $GC(\beta)$, one derives $\mathbb{P}[f(Y_1) - \mathbb{E}[f(Y_1)] \ge r] \le \exp(-\lambda r + \frac{\beta\lambda^2}{4})$, $\forall \lambda, r \ge 0$. An optimisation over $\lambda$ gives that this probability is bounded by $\exp(-\frac{r^2}{\beta})$.

A practical criterion for $GC(\beta)$ to hold is given by Bolley and Villani [24]. If there exists $\varepsilon > 0$ s.t. $\int \exp(\varepsilon|y|^2)\mu(dy) < \infty$, then $\mu$ satisfies $GC(\beta)$ for some $\beta := \beta(\varepsilon)$. The two claims are actually equivalent.

In [7] and [8], $GC(\beta)$ is the only crucial property we require on the underlying innovations of the approximation schemes.

Such transport-entropy inequalities are very attractive especially from a numerical point of view since they are related to the concentration of measure phenomenon which allows to establish non-asymptotic deviation estimates. In particular, some *non-asymptotic confidence intervals* can be explicitly computed. This can be very useful in practice when the computational resources are constrained.

The three next results put an emphasis on this point. Suppose that $(X^{(n)})_{n \ge 1}$ is a sequence of i.i.d. $\mathbb{R}^d$-valued random variables with common law $\mu$.

**Corollary 3.1.2.** If $\mu$ satisfies $T_\Phi$ then for all 1-Lipschitz function $f: \mathbb{R}^d \to \mathbb{R}$ and for all $r \ge 0$, for all $M \ge 1$, one has

$$\mathbb{P} \left( \left| \frac{1}{M} \sum_{k=1}^{M} f(X^{(k)}) - \mathbb{E}[f(X^{(1)})] \right| \ge r \right) \le 2 \exp(-M\Phi(r))$$

Deriving non-asymptotic deviation bounds for $W_1(\mu_M, \mu)$ is of interest for many applications in the fields of numerical probability and statistic. In its present form, next result is due to Gozlan and Leonard [55], Theorem 12.

**Proposition 3.1.3.** If $\mu$ satisfies $T_\Phi$ then the empirical measure $\mu_M$ defined as $\mu_M = \frac{1}{M} \sum_{k=1}^M \delta_{X^{(k)}}$ satisfies the following concentration bound

$$\mathbb{P}(W_1(\mu_M, \mu) \ge \mathbb{E}[W_1(\mu_M, \mu)] + r) \le \exp(-M\Phi(r)).$$

where for $x \in \mathbb{R}^d$, $\delta_x$ stands for the Dirac mass at point $x$.
The quantity $\mathbb{E}[W_1(\mu_M, \mu)]$ will go to zero as $M$ goes to infinity, by convergence of empirical measures, but we still need quantitative bounds. The next result is an adaptation of Theorem 10.2.1 in [116] on similar bounds but for the distance $W_2$. A proof can be found in the appendix of [8].

**Proposition 3.1.4.** Assume that $\mu$ has a finite moment of order $d+3$. Then, one has

$$ \mathbb{E}[W_1(\mu_M, \mu)] \le C(d, \mu)M^{-1/(d+2)} $$

where

$$ C(d, \mu) := 4\sqrt{d} + 2\sqrt{\int_{\mathbb{R}^d} (1 + |x|^{d+1})^{-1} dx} \sqrt{2^{-2d} + 2^{3-d} \int |y|^{d+3}\mu(dy) + 2^{3-d}d(d+3)!} $$

This bound is not optimal in general, but has the advantage of having very explicit constants. In the case of a distribution with compact support, it has been shown in [12], Section 7, that $\mathbb{E}[W_1(\mu_M, \mu)]$ is of order $O(M^{-1/d})$, and that this is the optimal exponent in $d$ when $d \ge 3$. In view of Kantorovich-Rubinstein duality formula, namely

$$ W_1(\mu, \nu) = \sup \left\{ \int fd\mu - \int fd\nu : [f]_1 \le 1 \right\} $$

where $[f]_1$ denotes the Lipschitz-modulus of $f$, the latter result provides the following concentration bounds

$$ \forall r \ge 0, \forall M \ge 1, \quad \mathbb{P} \left( \sup_{f: [f]_1 \le 1} \left( \frac{1}{M} \sum_{k=1}^{M} f(X^{(k)}) - \mathbb{E}[f(X^{(1)})] \right) \ge C(d, \mu) M^{-1/(d+2)} + r \right) \le \exp(-M\Phi(r)). $$

Similar results were first obtained for different concentration regimes by Bolley, Guillin and Villani [23] relying on a non-asymptotic version of Sanov's Theorem. Some of these results have also been derived by Boissard [22] using concentration inequalities, and were also extended to ergodic Markov chains up to some contractivity assumptions in the Wasserstein metric on the transition kernel.

Some applications are proposed in Bolley, Guillin and Villani [23]. Such results can indeed provide non-asymptotic deviation bounds for the estimation of the density of the invariant measure of a Markov chain. Let us note that the (possibly large) constant $C(d, \mu)$ appears as a trade-off to obtain uniform deviations over all Lipschitz functions.

In the two next sections we will present some transport-entropy inequalities and, as a consequence, some non-asymptotic deviation estimates for the laws at a given time step of two kinds of discrete-time and $d$-dimensional stochastic evolution scheme of the form

$$ X_{n+1} = X_n + \gamma_{n+1} H(n, X_n, U_{n+1}), \quad n \ge 0, X_0 = x \in \mathbb{R}^d, \qquad (3.1) $$

where $(\gamma_n)_{n \ge 1}$ is a deterministic positive sequence of time steps, $(U_i)_{i \in \mathbb{N}^*}$ is an i.i.d. $\mathbb{R}^q$-valued sequence of random variables defined on some probability space $(\Omega, \mathcal{F}, \mathbb{P})$ with law $\mu$ and the function $H: \mathbb{N} \times \mathbb{R}^d \times \mathbb{R}^q \to \mathbb{R}^d$ is a measurable function satisfying for all $x \in \mathbb{R}^d$, for all $n \in \mathbb{N}$, $H(n, x, \cdot) \in \mathcal{L}^1(\mu)$, and $\mu(du)$-a.s., $H(n, \cdot, u)$ is continuous. We will also assume that $\mu$ satisfies the *Gaussian concentration property* $GC(\beta)$.

As a consequence of the transport-entropy inequalities obtained for the laws at a given time step of Euler like schemes and stochastic approximation algorithm, we will derive in Section 3.3 non-asymptotic deviation bounds in the Wasserstein metric.

## 3.2 Gaussian concentration inequalities

In this section, we present the results obtained in [7] in collaboration with S. Menozzi from the university of Evry Val d'Essonne related to Gaussian concentration inequalities for two kinds of recursive scheme with dynamics (3.1). In Section 3.3, we will extend those results to other regimes ranging from exponential to Gaussian.

We consider the unique strong solution $X$ to the SDE (1.17) where the underlying innovation process $Z$ is a $q$-dimensional Brownian motion $(W_t)_{t \ge 0}$. In [7], we assume that the coefficients satisfy the mild smoothness condition:

(A) The coefficients $b, \sigma$ are uniformly Lipschitz continuous in space uniformly in time, $\sigma$ is bounded.

Note that we do not assume any non-degeneracy condition on $\sigma$ in (A).

We do not consider the associated continuous Euler Maruyama scheme (1.19) but rather define an Euler like scheme. For a fixed time step $\Delta = T/N$, $N \in \mathbb{N}^*$, and the time grid $t_i := i\Delta$, for all $i = 0, \dots, N$, we consider the following dynamics

$$ X_0^\Delta = x, \quad X_{t_{i+1}}^\Delta = X_{t_i}^\Delta + b(t_i, X_{t_i}^\Delta)\Delta + \sigma(t_i, X_{t_i}^\Delta)\Delta^{1/2}U^{(i+1)}, \quad i = 0, \dots, N-1 \qquad (3.2) $$
where $(U^{(i)})_{i \in \mathbb{N}^*}$ is a sequence of $\mathbb{R}^q$-valued i.i.d. random variables with law $\mu$ satisfying: $\mathbb{E}[U^{(1)}] = 0_q$, $\mathbb{E}[U^{(1)}(U^{(1)})^*] = I_q$, where $(U^{(1)})^*$ denotes the transpose of the column vector $U^{(1)}$ and $0_q, I_q$ respectively denote the zero vector of $\mathbb{R}^q$ and the identity matrix of $\mathbb{R}^q \otimes \mathbb{R}^q$. We also assume that $\mu$ satisfies $\mathrm{GC}(\beta)$ for some $\beta > 0$. The main advantage of such a situation is that it includes the case of the standard Euler scheme where $U^{(1)} \stackrel{d}{=} \mathcal{N}(0, I_q)$ and the case of the Bernoulli law where $U^{(1)} \stackrel{d}{=} (B_1, \dots, B_q)$, $(B_k)_{1 \le k \le q}$ are i.i.d random variables with law $\mu = \frac{1}{2}(\delta_{-1} + \delta_1)$, both satisfying $\mathrm{GC}(\beta)$ with $\beta = 2$.

As already mentioned, for practical purposes, we are interested in non-asymptotical controls of the statistical error, which is given by $\frac{1}{M} \sum_{i=1}^M f(X_T^{\Delta,(i)}) - \mathbb{E}[f(X_T^{\Delta})]$ where the $((X_T^{\Delta,(i)})^i)_{1 \le i \le M}$ are independent copies of the scheme (3.2) starting at $x$ at time 0 and evaluated at time $T$. Precisely, for a fixed $M$ and a given threshold $r > 0$, one would like to give bounds on the quantity $\mathbb{P}[\frac{1}{M} \sum_{i=1}^M f(X_T^{\Delta,(i)}) - \mathbb{E}[f(X_T^{\Delta})] > r]$.

Our first main result obtained in [7] shows that when the innovations satisfy $\mathrm{GC}(\beta)$, the Gaussian concentration property transfers to the statistical error.

**Theorem 3.2.1 (Concentration Bounds for an Euler-like scheme).** Denote by $X_T^\Delta$ the value at time $T$ of the scheme (3.2). Assume that the common law of the i.i.d. sequence $(U_i)_{i \in \mathbb{N}^*}$ in (3.2) satisfies $\mathrm{GC}(\beta)$ for some $\beta > 0$ and that the coefficients $b, \sigma$ satisfy (A). Let $f$ be a real valued uniformly Lipschitz continuous function on $\mathbb{R}^d$. For all $M \in \mathbb{N}^*$ and all $r \ge 0$, one has

$$ \mathbb{P}[\left|\frac{1}{M}\sum_{i=1}^{M} f(X_T^{\Delta,(i)}) - \mathbb{E}[f(X_T^{\Delta})]\right| \ge r] \le 2 \exp\left(-\frac{r^2 M}{T \Psi(T, f, b, \sigma, q)}\right), \\ \Psi(T, f, b, \sigma, q) := 4\beta[f]_1^2 |\sigma|_\infty^2 \exp(2([b]_1 + c[\sigma])_1(1 + c[\sigma])_1 T), $$

where $q$ is the dimension of the underlying Brownian motion $W$ and $c := c(q)$.

Note that contrary to Lemaire and Menozzi [91] we do not need any non-degeneracy condition on the diffusion coefficient and we got rid off the systematic bias.

The strategy to prove the above concentration inequality is similar to the first step of the analysis of the discretisation error of Talay and Tubaro [122] presented in Section 1.2.3. Indeed, introducing the functions $v^\Delta(t_i, x) := \mathbb{E}[f(X_T^\Delta)|X_{t_i}^\Delta = x]$, $0 \le i \le N$, $x \in \mathbb{R}^d$, we use the Markov property satisfied by the scheme and the following telescopic sum

$$ \begin{aligned} f(X_T^\Delta) - \mathbb{E}[f(X_T^\Delta)] &= \sum_{i=1}^{N} v^\Delta(t_i, X_{t_i}^\Delta) - \mathbb{E}[v^\Delta(t_i, X_{t_i}^\Delta)|X_{t_{i-1}}^\Delta] \\ &= \sum_{i=1}^{N} f_i^\Delta(X_{t_{i-1}}^\Delta, \sqrt{\Delta}U^{(i)}) - \mathbb{E}[f_i^\Delta(X_{t_{i-1}}^\Delta, \sqrt{\Delta}U^{(i)})|X_{t_{i-1}}^\Delta] \end{aligned} \quad (3.3) $$

where $f_i^\Delta(x, u) := \mathbb{E}[f(X_T^\Delta)|X_{t_i}^\Delta = x + b(t_{i-1}, x)\Delta + \sigma(t_{i-1}, x)u]$, $1 \le i \le N$, $(x, u) \in \mathbb{R}^d \times \mathbb{R}^q$. In Talay and Tubaro [122], $v^\Delta(t_i, X_{t_i}^\Delta)$ is replaced by $v(t_i, X_{t_i}^\Delta) = \mathbb{E}[f(X_T)|X_{t_i} = x]_{x=X_{t_i}^\Delta}$, that is the expectation involving the diffusion at time $T$ starting from the current value of the scheme at $t_i$. Direct computations based on B-D-G inequality and Gronwall's Lemma show that $u \mapsto f_i^\Delta(x, u)$ is Lipschitz-continuous with $[f_i^\Delta(x, \cdot)]_1 \le (\Psi(T, f, b, \sigma, q))^{1/2}$, $x \in \mathbb{R}^d$.

Since $U^{(i)}$ satisfies $\mathrm{GC}(\beta)$, from Proposition 3.1.1 it follows that

$$ \mathbb{E}[\exp(\lambda(f_i^\Delta(X_{t_{i-1}}^\Delta, \sqrt{\Delta}U^{(i)}) - \mathbb{E}[f_i^\Delta(X_{t_{i-1}}^\Delta, \sqrt{\Delta}U^{(i)})|X_{t_{i-1}}^\Delta]))|X_{t_{i-1}}^\Delta] \le \exp(\beta([f_i^\Delta(X_{t_{i-1}}^\Delta, \cdot)]_1 \sqrt{\Delta}\lambda)^2/4), \quad i = 1, \dots, N. \quad (3.4) $$

Now using the Markov exponential inequality, the decomposition (3.3) and the tower property of condition expectation one proves Theorem 3.2.1 for $M = 1$. The bound for all positive integer $M$ follows from a simple tensorization argument for independent random variables satisfying the Gaussian concentration property.

Besides our consideration for the Euler scheme, we also consider a Robbins-Monro algorithm with dynamics

$$ \theta_{n+1} = \theta_n - \gamma_{n+1} H(\theta_n, U^{(n+1)}), \quad n \ge 0, \quad \theta_0 \in \mathbb{R}^d, \quad (3.5) $$

where $(U^{(n)})_{n \ge 1}$ is an i.i.d. $\mathbb{R}^q$-valued sequence of random variables with law $\mu$ satisfying $\mathrm{GC}(\beta)$ for some $\beta > 0$ and $(\gamma_n)_{n \ge 1}$ is a sequence of non-negative deterministic steps satisfying the usual assumption (1.9). We assume that the following conditions on the function $H$ and the step sequence $(\gamma_n)_{n \ge 1}$ in (4.1) are in force:

(HL) The map $(\theta, u) \in \mathbb{R}^d \times \mathbb{R}^q \mapsto H(\theta, u)$ is uniformly Lipschitz continuous.

(HUA) The map $h: \theta \in \mathbb{R}^d \mapsto \mathbb{E}[H(\theta, U)]$ is continuously differentiable in $\theta$ and there exists $\lambda > 0$ s.t. $\forall \theta \in \mathbb{R}^d, \forall \xi \in \mathbb{R}^d, \lambda|\xi|^2 \le \langle Dh(\theta)\xi, \xi \rangle$ (Uniform Attractivity).
As we have seen in Section 1.1.3 of the introductory chapter, in order to derive a Central Limit Theorem for the sequence $(\theta_n)_{n \ge 1}$ as described in Duflo [39] or Kushner and Yin [81], it is commonly assumed that the matrix $Dh(\theta^*)$ is uniformly attractive. In our current framework, this local condition on the Jacobian matrix of $h$ at the equilibrium is replaced by the uniform assumption (HUA). In the framework of stochastic gradient algorithm, this assumption is also known as uniform convexity.

Here the global error between the stochastic approximation procedure $\theta_n$ at a given time step $n$ and its target $\theta^*$ can be decomposed as follows:

$$ 
\begin{aligned}
z_n &:= |\theta_n - \theta^*| \\
&:= (|\theta_n - \theta^*| - \mathbb{E}[|\theta_n - \theta^*|]) + \mathbb{E}[|\theta_n - \theta^*|] \\
&:= \mathcal{E}_{\text{Emp}}(\gamma, n, H, \lambda, \alpha) + \delta_n
\end{aligned} 
\quad (3.6) $$

where $\delta_n := \mathbb{E}[|\theta_n - \theta^*|]$.

The term $\mathcal{E}_{\text{Emp}}(\gamma, n, H, \lambda, \alpha)$ corresponding to the difference between the absolute value of the error at time $n$ and its mean can be viewed as an empirical error. The term $\delta_n$ is an intrinsic bias coming from the non-linear nature of the algorithm. Note that from Proposition 1.14 in Section 1.1.3 under (HL) and (HUA) one has

$$ \delta_N \le \mathbb{E}[|\theta_N - \theta^*|^2]^{\frac{1}{2}} \le C(\gamma_N)^{\frac{1}{2}}, \quad N \ge 1 $$

for the two following step sequences: $\gamma_n = c/n^\beta$, $\beta \in (1/2, 1)$, and $\gamma_n = c/n$, with $c > 1/(2\lambda)$, $n \ge 1$.

An explicit expression for the constant $C$ is provided in [7]. As for the Euler scheme, we prove that the Gaussian concentration property transfers to the empirical error under (HL) and (HUA).

**Theorem 3.2.2 (Concentration Bounds for Robbins-Monro algorithms).** Assume that the function $H$ of the recursive procedure $(\theta_n)_{n \ge 0}$ (with starting point $\theta_0 \in \mathbb{R}^d$) defined by (4.1) satisfies (HL) and (HUA), and that the step sequence $(\gamma_n)_{n \ge 1}$ satisfies the usual step assumption (1.9). Suppose that the law of the innovation satisfies $GC(\beta)$, $\beta > 0$. Then, for all $N \in \mathbb{N}^*$ and all $r \ge 0$,

$$ \mathbb{P}(|\theta_N - \theta^*| \ge r + \delta_N) \le \exp \left( - \frac{r^2}{\beta [H]_1^2 \Pi_N \sum_{k=1}^N \gamma_k^2 / \Pi_k} \right) $$

where $\Pi_N := \prod_{k=0}^{N-1} (1 - 2\lambda\gamma_{k+1} + [H]_1^2\gamma_{k+1}^2)$.

Concerning the choice of the step sequence $(\gamma_n)_{n \ge 1}$ and its impact on the concentration rate and bias, we obtain the following results:

* If we choose $\gamma_n = \frac{c}{n}$, with $c > 0$. Then, $\Gamma_N = \sum_{k=1}^N \gamma_k = c \log(N) + c_1' + r_N$, $c_1' > 0$ and $r_N \to 0$, so that $\Pi_N = (N^{-2c\lambda})$.

- If $c < \frac{1}{2\lambda}$, we obtain $\Pi_N \sum_{k=1}^N \gamma_k^2 / \Pi_k = \mathcal{O}(N^{-2c\lambda})$.

- If $c > \frac{1}{2\lambda}$, a comparison between the series and the integral yields $\Pi_N \sum_{k=1}^N \gamma_k^2 / \Pi_k = (N^{-1})$.

* If we choose $\gamma_n = \frac{c}{n^\rho}$, $c > 0$, $\frac{1}{2} < \rho < 1$, then one has $\Pi_N \sum_{k=1}^N \gamma_k^2 \Pi_k^{-1} = o(N^{-\rho+\epsilon})$, $\epsilon \in (0, 1-\rho)$.

As for the Euler like scheme, the strategy consists in introducing again a telescopic sum of conditional expectations. Denoting for all $i \in \mathbb{N}$, $\mathcal{F}_i := \sigma(U^{(j)}, j \le i)$, i.e. $(\mathcal{F}_i)_{i \in \mathbb{N}}$ is the natural filtration of the algorithm, we write for all $n \in \mathbb{N}^*$:

$$ 
\begin{aligned}
\mathcal{E}_{\text{Emp}}(\gamma, N, H, \lambda, \alpha) &= |z_N| - \mathbb{E}[|z_N|] \\
&= \sum_{i=1}^{N} \mathbb{E}[|z_N| | \mathcal{F}_i] - \mathbb{E}[|z_N| | \mathcal{F}_{i-1}] \\
&= \sum_{i=1}^{N} v_i^{\gamma}(\theta_i) - \mathbb{E}[v_i^{\gamma}(\theta_i) | \mathcal{F}_{i-1}], \\
&= \sum_{i=1}^{N} f_i^{\gamma}(\theta_{i-1}, U^{(i)}) - \mathbb{E}[f_i^{\gamma}(\theta_{i-1}, U^{(i)}) | \mathcal{F}_{i-1}],
\end{aligned} $$

where we used the Markov property for the second equality and we introduced the notations $v_i^{\gamma}(\theta) := \mathbb{E}[|\theta_n - \theta^*| | \theta_i = \theta]$, $1 \le i \le n$, $\theta \in \mathbb{R}^d$, $f_i^{\gamma}(\theta, u) = v_i^{\gamma}(\theta - \gamma_i H(\theta, u))$. The stability of the Gaussian concentration property is then derived using that the $f_i^{\gamma}$ are Lipschitz with respect to the variable $u$ and similar arguments to those employed for the Euler like scheme.
## 3.3 From Exponential to Gaussian concentration regimes

In [7], we obtained Gaussian concentration inequalities for the law of an Euler like scheme of a diffusion process and for the deviation of a Robbins-Monro algorithm from its target. It should be noted that it is the boundedness of the diffusion coefficient $\sigma$ and the uniform Lipschitz continuity of the function $u \mapsto H(\theta, u)$, $\theta \in \mathbb{R}^d$, that give rise to the Gaussian concentration regime for the deviation of the empirical error. In this section, we present the results obtained [8] in collaboration with M. Fathi from the university Paul Sabatier that extend and complete those exposed in the previous section. The key point is to quantify to properly quantify the contribution of the diffusion term to the concentration regime. We also obtain non-asymptotic bound for stochastic approximation with averaging of trajectories, in particular we prove that averaging a stochastic approximation algorithm with a slow decreasing step sequence gives rise to optimal concentration rate.

We consider the Euler-like scheme $(X_{t_i}^\Delta)_{0 \le i \le N}$ with dynamics (3.2). The innovation sequence $(U^{(i)})_{1 \le i \le N}$ satisfies $GC(\beta)$ for some $\beta > 0$. In [8], we replace assumption (A) by the following smoothness and domination assumptions

(HS) The coefficients $b, \sigma$ are uniformly Lipschitz continuous in space uniformly in time.

(HD$_\alpha$) There exists a $C^2(\mathbb{R}^d, \mathbb{R}_+^*)$ function $V$ satisfying $\exists C_V > 0, |\nabla V|^2 \le C_V V$, $\eta := \frac{1}{2} \sup_{x \in \mathbb{R}^d} \|\nabla^2 V(x)\| < +\infty$ and $\exists \alpha \in (0, 1]$, such that for all $x \in \mathbb{R}^d$,

$$ \exists C_b > 0, \quad \sup_{t \in [0, T]} |b(t, x)|^2 \le C_b V(x), \quad \exists C_\sigma > 0, \quad \sup_{t \in [0, T]} \mathrm{Tr}(a(t, x)) \le C_\sigma V^{1-\alpha}(x). $$

where $a = \sigma\sigma^*$.

The idea behind assumption (HD$_\alpha$) is to parameterize the growth of the diffusion coefficient in order to quantify its contribution to the concentration regime. Indeed, under (HS) and (HD$_\alpha$), with $\alpha \in [1/2, 1]$ we derive non-asymptotic deviation bounds for the statistical error ranging from exponential (if $\alpha = 1/2$) to Gaussian (if $\alpha = 1$) regimes.

**Theorem 3.3.1** (Transport-Entropy inequalities for Euler like schemes). Denote by $X_T^\Delta$ the value at time $T$ of the scheme (3.2) and by $\mu_T^\Delta$ its law. Denote the Lipschitz modulus of $b$ and $\sigma$ by $[b]_1$ and $[\sigma]_1$ respectively. Assume that the common law of the innovation sequence $(U^{(i)})_{1 \le i \le N}$ in (3.2) satisfies $GC(\beta)$ for some $\beta > 0$ and that the coefficients $b, \sigma$ satisfy (HS) and (HD$_\alpha$) for some $\alpha \in [1/2, 1]$.

Then, $\mu_T^\Delta$ satisfies $T_{\Phi_\alpha^*}(\lambda)$ with $\Phi_\alpha^*(\lambda) = \sup_{\rho \ge 0} \{\lambda\rho - \Phi_\alpha(\rho)\}$ and one has:

* If $\alpha \in (1/2, 1]$, for all $\rho \ge 0$

$$ \Phi_{\alpha}(\rho) = \Psi_{\alpha}(T, \Delta, b, \sigma, x)(\rho^2 \vee \rho^{\frac{2\alpha}{2\alpha-1}}), $$

* If $\alpha = 1/2$, for all $\rho \in [0, \varphi(T, b, \sigma, \Delta)^{-1/2}\lambda_{3,2})$

$$ \Phi_{1/2}(\rho) = K_{3,2} \frac{(\rho\varphi(T, b, \sigma, \Delta)^{1/2}/\lambda_{3,2})^2}{1 - (\rho\varphi(T, b, \sigma, \Delta)^{1/2}/\lambda_{3,2})}. $$

Moreover, we have $\Psi_\alpha(T, \Delta, b, \sigma, x) = K_{3,1}(\varphi(T, b, \sigma, \Delta)^2 \vee \varphi(T, b, \sigma, \Delta)^{\frac{\alpha}{2\alpha-1}})$, $\varphi(T, b, \sigma, \Delta) = C_\sigma \beta \frac{(1+C(\Delta)\Delta)}{4C(\Delta)} e^{3C(\Delta)T}$, $C(\Delta) := 2[b]_1 + [\sigma]_1^2 + \Delta[b]_1^2$ where $K_{3,1}, \lambda_{3,2}$ and $K_{3,2}$ are explicit constants.

* In the case $\alpha \in (1/2, 1]$, one easily gets the following explicit formula:

- If $\lambda \in [0, 2\Psi]$, then $\Phi_\alpha^*(\lambda) = \frac{1}{4\Psi}\lambda^2$;

- If $\lambda \in [\frac{2\alpha}{2\alpha-1}\Psi, +\infty)$, then $\Phi_\alpha^*(\lambda) = \frac{1}{2\alpha}(\frac{2\alpha-1}{2\alpha\Psi})^{2\alpha-1}\lambda^{2\alpha}$;

- If $\lambda \in (2\Psi, \frac{2\alpha}{2\alpha-1}\Psi)$, then $\Phi_\alpha^*(\lambda) = \lambda - \Psi$.

Let us note that the linear behavior of $\Phi_\alpha^*$ on a small interval is due to the fact that $\Phi_\alpha$ is not $C^1$.

* In the case $\alpha = 1/2$, tedious but simple computations show that

$$ \Phi_{1/2}^*(\lambda) = \left( \left( 1 + \frac{\lambda_{3,2}}{K_{3,2}\varphi(T, b, \sigma, \Delta)^{1/2}} \right)^{\frac{1}{2}} - 1 \right)^2. $$

We also consider the Robbins-Monro algorithm with dynamics (4.1) under the assumption that the innovation sequence $(U^{(i)})_{1 \le i \le N}$ satisfies $GC(\beta)$ for some $\beta > 0$. Instead of (HL) we will make the following assumptions:
(HL) For all $u \in \mathbb{R}^q$, the function $H(\cdot, u)$ is Lipschitz-continuous with a Lipschitz modulus having linear growth in the variable $u$, that is:

$$ \exists C_H > 0, \forall u \in \mathbb{R}^q, \sup_{(\theta, \theta') \in (\mathbb{R}^d)^2} \frac{|H(\theta, u) - H(\theta', u)|}{|\theta - \theta'|} \le C_H(1 + |u|). $$

(HLS)α (Lyapunov Stability-Domination) There exists a $C^2(\mathbb{R}^d, \mathbb{R}_+^*)$ function $L$ satisfying $\exists C_L > 0, |\nabla L|^2 \le C_L L, \eta := \frac{1}{2} \sup_{x \in \mathbb{R}^d} \|\nabla^2 L(x)\| < +\infty$ such that

$$ \forall \theta \in \mathbb{R}^d, \langle \nabla L(\theta), h(\theta) \rangle \ge 0, \quad \text{and} \quad \exists C_h > 0, \forall \theta \in \mathbb{R}^d, |h(\theta)|^2 \le C_h L(\theta). $$

and $\exists \alpha \in (0, 1]$,

$$ \exists C_{\alpha} > 0, \forall \theta \in \mathbb{R}^d, \sup_{(u, u') \in (\mathbb{R}^q)^2} \frac{|H(\theta, u) - H(\theta, u')|}{|u - u'|} \le C_{\alpha} L^{\frac{1-\alpha}{2}}(\theta). $$

Assumption (HLS)α is similar to assumption (HD)α for the Euler like scheme and allows to quantify the contribution of the diffusion term to the concentration regime of the empirical error $\mathcal{E}_{Emp}(\gamma, N, H, \lambda, \alpha)$ of the Robbins-Monro scheme.

**Theorem 3.3.2** (Transport-Entropy inequalities for stochastic approximation algorithms). Let $N \in \mathbb{N}^*$. Assume that the function $H$ of the recursive procedure $(\theta_n)_{0 \le n \le N}$ (with starting point $\theta_0 \in \mathbb{R}^d$) defined by (4.1) satisfies (HL), (HUA) and (HLS)α for $\alpha \in [\frac{1}{2}, 1]$, and that the step sequence $\gamma = (\gamma_n)_{n \ge 0}$ satisfies (1.9). Suppose that the law of the innovation satisfies $GC(\beta)$, $\beta > 0$. Denote by $\mu_N^\gamma$ the law of $\theta_N$.

Then, $\mu_N^\gamma$ satisfies $T_{\Phi_N^*}$ with $\Phi_{\alpha,N}^*(\lambda) = \sup_{\rho \ge 0} \{\lambda \rho - \Phi_{\alpha,N}(\rho)\}$ and one has:

* If $\alpha \in (\frac{1}{2}, 1]$, for all $\rho \ge 0$

$$ \Phi_{\alpha,N}(\rho) = \varphi_{\alpha}(\gamma, H, \theta_0)(C_N^{\gamma} \rho^2 \lor C_N^{\gamma, \alpha} \rho^{\frac{2\alpha}{2\alpha-1}}). $$

* If $\alpha = \frac{1}{2}$, for all $\rho \in [0, \lambda_{4.1}/\tilde{s}_N)$,

$$ \Phi_{1/2,N}(\rho) = 2\varphi_{1/2}(\gamma, H, \theta_0)C_N^{\gamma} \frac{(\rho/\lambda_{4.1})^2}{1 - (\rho\tilde{s}_N/\lambda_{4.1})}. $$

Moreover the three concentration rate sequences are defined for $N \in \mathbb{N}^*$ by

$$ C_N^\gamma := \sum_{k=0}^{N-1} \gamma_{k+1}^2 \frac{\Pi_{1,N}}{\Pi_{1,k}}, $$

$$ C_N^{\gamma, \alpha} := \sum_{k=0}^{N-1} \gamma_{k+1}^{\frac{2\alpha}{2\alpha-1}} \left( \frac{\Pi_{1,N}}{\Pi_{1,k}} \right)^{\frac{2\alpha}{2\alpha-1}} \left( (k+1) \log^2(k+4) \right)^{\frac{1-\alpha}{2\alpha-1}} $$

$$ \tilde{s}_N := \max_{0 \le k \le N-1} (k+1)^{1/2} \log(k+4) \gamma_{k+1} \left( \frac{\Pi_{1,N}}{\Pi_{1,k}} \right)^{\frac{1}{2}} \exp\left( \sum_{p=0}^{N-1} \frac{1}{(p+1) \log^2(p+4)} \right) $$

with $\Pi_{1,N} := \prod_{k=0}^{N-1} (1 - 2\lambda\gamma_{k+1} + C_{H,\mu}\gamma_{k+1}^2)$ and where $C_{H,\mu}$ and $\varphi := \varphi_\alpha(\gamma, H, \theta_0)$ are explicitly constants.

As in the case of Euler like schemes, for $\alpha \in (\frac{1}{2}, 1]$, we have:

* if $\lambda \in [0, 2\varphi(C_N^\gamma/(C_N^{\gamma,\alpha})^{2\alpha-1})^{\frac{1}{2(1-\alpha)}}]$, then $\Phi_{\alpha,N}^*(\lambda) = \lambda^2/(4\varphi C_N^\gamma)$;

* If $\lambda \in [\frac{2\alpha}{2\alpha-1}, \varphi(C_N^\gamma/(C_N^{\gamma,\alpha})^{2\alpha-1})^{\frac{1}{2(1-\alpha)}}, +\infty)$, then $\Phi_{\alpha,N}^*(\lambda) = \frac{1}{2\alpha} \left( \frac{2\alpha-1}{2\alpha\varphi} \right)^{2\alpha-1} (\lambda^{2\alpha}/(C_N^{\gamma,\alpha})^{2\alpha-1})$;

* If $\lambda \in (2\varphi(C_N^\gamma/(C_N^{\gamma,\alpha})^{2\alpha-1})^{\frac{1}{2(1-\alpha)}}, \frac{2\alpha}{2\alpha-1}, \varphi(C_N^\gamma/(C_N^{\gamma,\alpha})^{2\alpha-1})^{\frac{1}{2(1-\alpha)}})$, then $\Phi_{\alpha,N}^*(\lambda) = (\frac{C_N^\gamma}{C_N^{\gamma,\alpha}})^{\frac{2\alpha-1}{2(1-\alpha)}} \lambda - \varphi(\frac{(C_N^\gamma)^{\frac{\alpha}{1-\alpha}}}{(C_N^{\gamma,\alpha})^{\frac{2\alpha-1}{1-\alpha}}})$.

For $\alpha = \frac{1}{2}$, we obtain the following explicit bound for the Legendre transform of $\Phi_{1/2,N}$

$$ \forall \lambda \ge 0, \quad \Phi_{1/2,N}^*(\lambda) = \frac{2\varphi C_N^\gamma}{\tilde{s}_N^2} \left( \left( 1 + \frac{\tilde{s}_N \lambda_{4.1} \lambda}{2\varphi C_N^\gamma} \right)^{\frac{1}{2}} - 1 \right)^2 $$

Hence, for $N \ge 1$ being fixed, the following simple asymptotic behaviours can be easily derived:
* When $\lambda$ is small, $\Phi_{1/2,N}^*(\lambda) \sim \lambda_{4.1}^2 \lambda^2/(2\varphi C_N^\gamma)$;

* When $\lambda$ goes to infinity, $\Phi_{1/2}^*(\lambda) \sim \lambda_{4.1} \lambda / \tilde{s}_N$.

Concerning the choice of the step sequence $(\gamma_n)_{n \ge 1}$ and its impact on the concentration rate and bias, we obtain the following results:

* If we choose $\gamma_n = c/n$, with $c > 0$.

- If $c < \frac{1}{2\lambda}$, one has $C_N^\gamma = \mathcal{O}(N^{-2c\lambda})$, $C_N^{\gamma,\alpha} = \mathcal{O}(N^{-\frac{2\alpha}{2\alpha-1}c\lambda})$, $\tilde{s}_N = (N^{-c\lambda})$.

- If $c > \frac{1}{2\lambda}$, one has $C_N^\gamma = \mathcal{O}(N^{-1})$, $C_N^{\gamma,\alpha} = \mathcal{O}((\log(N))^2 \frac{1-\alpha}{2\alpha-1} N^{-\frac{\alpha}{2\alpha-1}})$, $\tilde{s}_N = \mathcal{O}(\log(N)N^{-\frac{1}{2}})$.

* If we choose $\gamma_n = \frac{c}{n\rho}$, $c > 0$, $\frac{1}{2} < \rho < 1$, then one has $C_N^\gamma = o(N^{-\rho+\epsilon})$, $\epsilon \in (0, 1-\rho)$, $C_N^{\gamma,\alpha} = o(N^{-\frac{(\rho-(1-\alpha))}{2\alpha-1}-\epsilon})$ and $\tilde{s}_N = \mathcal{O}(\log(N)N^{-\frac{\rho-1}{2}})$.

We see that the optimal convergence rate $\mathcal{O}(N^{-1})$ and $\mathcal{O}((\log(N))^2 \frac{1-\alpha}{2\alpha-1} N^{-\frac{\alpha}{2\alpha-1}})$ is obtained by choosing $\gamma_n = c/n$ for $c > 1/(2\lambda)$. Let us notice that we find the same critical level for the constant $c$ as in the CLT for stochastic algorithms, see Theorem 1.1.5. Indeed, if $c > \frac{1}{2\mathcal{R}e(\lambda_{min})}$ where $\lambda_{min}$ denotes the eigenvalue of $Dh(\theta^*)$ with the smallest real part then we know that the optimal rate $\sqrt{N}$ is achieved in the CLT.

In order to get rid of this constraint, we also proved in [8] a transport entropy inequality for the law $\tilde{\mu}_N^\gamma$ of the empirical mean $(\bar{\theta}_n)_{0 \le n \le N}$ of the sequence $(\theta_n)_{0 \le n \le N}$ defined by (4.1) when the step sequence is given by $\gamma_n = c/n^\rho$, $\frac{1}{2} < \rho < 1$. In particular, we prove that averaging the original Robbins-Monro dynamics with this slow decreasing step allows to achieve for free the optimal concentration rates, that is, the one obtained with $\gamma_n = c/n$ with $c > 1/(2\lambda)$.

We briefly explain the strategy used in [8] to prove this new transport entropy. The first step is similar to the one presented in the previous section. We consider the telescopic sum (3.3). However, under $(\mathbf{HD}_\alpha)$ the bound (3.4) becomes

$$ \mathbb{E}[\exp(\lambda(f_i^\Delta(X_{t_{i-1}}^\Delta, \sqrt{\Delta}U^{(i)}) - \mathbb{E}[f_i^\Delta(X_{t_{i-1}}^\Delta, \sqrt{\Delta}U^{(i)})|X_{t_{i-1}}^\Delta]))|X_{t_{i-1}}^\Delta] \le \exp(\beta\Psi_\alpha(T, \Delta, b, \sigma, x)\Delta \frac{\lambda^2}{4}V^{1-\alpha}(X_{t_{i-1}}^\Delta)), \quad i = 1, \dots, N. \quad (3.7) $$

Let $M_N = \sum_{i=1}^N v^\Delta(t_i, X_{t_i}^\Delta) - \mathbb{E}[v^\Delta(t_i, X_{t_i}^\Delta)|X_{t_{i-1}}^\Delta]$, $N \ge 1$ and $M_0 = 0$. From (3.7) and Hölder inequality, one gets

$$ \begin{aligned} \mathbb{E}[\exp(\lambda M_N)] & \le \mathbb{E}[\exp(\lambda M_{N-1}) \exp(\beta \Psi_\alpha(T, \Delta, b, \sigma, x) \Delta \frac{\lambda^2}{4} V^{1-\alpha}(X_{t_{N-1}}^\Delta))] \\ & \le \mathbb{E}[\exp(\lambda p M_{N-1})]^{\frac{1}{p}} \mathbb{E}[\exp(q \beta \Psi_\alpha(T, \Delta, b, \sigma, x) \Delta \frac{\lambda^2}{4} V^{1-\alpha}(X_{t_{N-1}}^\Delta))]^{\frac{1}{q}} \\ & \le \prod_{k=0}^{N-1} \left( \mathbb{E} \left[ \exp \left( \frac{\Psi_\alpha(T, \Delta, b, \sigma, x) \beta}{4} \lambda^2 q p^{2k} \Delta V^{1-\alpha} (X_{t_{N-k-1}}^\Delta) \right) \right] \right)^{\frac{1}{q}} \right)^{\frac{1}{pk}} \end{aligned} $$

where we used an induction argument for the last inequality. Hence, since $(p, q)$ are conjugate, one gets

$$ \mathbb{E}[\exp(\lambda(f(X_T^\Delta) - \mathbb{E}[f(X_T^\Delta)]))] \le \exp\left(\sup_{0 \le n \le N} \log\left(\mathbb{E}\left[e^{\frac{C\sigma\beta}{4}\lambda^2 q p^{2N}(1+C(\Delta)\Delta)^N V^{1-\alpha}(X_{t_n}^\Delta)}\right]\right)\right). $$

Finally, by stability arguments on the dynamics (3.2) with an innovation sequence satisfying $GC(\beta)$, $\beta > 0$ one gets

$$ \sup_{0 \le n \le N} \log \left( \mathbb{E} \left[ \exp \left( \lambda^2 V^{1-\alpha} (X_{t_n}^\Delta) \right) \right] \right) \le \begin{cases} K_{3.1} (\lambda \vee \lambda^{\frac{\alpha}{2\alpha-1}}) & \alpha \in (\frac{1}{2}, 1], \\ K_{3.2} \frac{(\lambda/\lambda_{3.2})^2}{1-(\lambda/\lambda_{3.2})}, & \alpha = \frac{1}{2}. \end{cases} $$

A suitable choice of $p$ and $q$ allows to conclude. The strategy of proof for the Robbins-Monro algorithm follows similar lines of arguments.

## 3.4 Perspective

An interesting problem is to investigate the case of non-regular test function $f$. In practical situations, notably in mathematical finance, one often has to deal with an irregular payoff function $f$ such as indicator functions. One possibility to address this problem is to allow $\sigma$ to be non-degenerated so that the Euler scheme admits a smooth
density. Then, one has to investigate the connection between the lack of regularity of $f$, the regularity of the density of the scheme and the resulting Lipschitz modulus of the function $y \mapsto f_i^{\Delta}(x, y)$.

Another interesting problem is to investigate concentration inequalities for the multi-level Monte Carlo estimator presented in Section 1.2.4. One expects that if the test function $f$ is Lipschitz-continuous then the Gaussian concentration regime still holds for this estimator but with a convergence rate given by the CLT established in [1]. The same problem may be considered for the multi-level stochastic approximation schemes presented in the next chapter.
# Chapter 4

# Improving the computational efficiency of stochastic approximation algorithms

In this chapter, we present [9] and [10], the latter being a joint work with L. Huang from the engineering school INSA of Toulouse. Both works aim to improve the computational efficiency for stochastic approximation schemes when one does not know how to directly sample from the exact distribution of the underlying innovation sequence. We extend the scope of the multi-level Monte Carlo and the multi-step Richardson Romberg methods briefly presented in Section 1.2.4 and Section 1.2.3.

The paradigm for both methods is described in Section 4.1. In Section 4.2 we present the results obtained in [9] related to multi-level methods for stochastic optimisation by means of stochastic approximation schemes. Contrary to the $L^2(\mathbb{P})$ analysis initiated by Giles [48], our methodology is based on establishing CLT for the proposed multi-level stochastic approximation estimators as proposed by Ben Alaya and Kebaier [1] in the standard Monte Carlo framework. In Section 4.3 we present the results obtained in [10] where we show that the principle of the multi-step Richardson-Romberg extrapolation method for Monte Carlo linear estimator as investigated in Pagès [109] can be extended to the framework of stochastic optimisation by means of stochastic approximation algorithm. For both methods we will see in particular that the new estimator outperforms the standard stochastic approximation estimator in terms of computational cost. Numerical results are presented in order to confirm the analysis.

## 4.1 Motivation

As was already mentioned in Section 1.1 of the introductory chapter, in many fields of applied mathematics such as computational finance, one often has to face some stochastic optimisation or zero search problems. We have seen some of them in Chapter 2. The stochastic approximation theory provides some recursive schemes in order find the zero of the function $h: \mathbb{R}^d \to \mathbb{R}^d$ defined by $h(\theta) := E[H(\theta, U)]$, $H: \mathbb{R}^d \times \mathbb{R}^q \to \mathbb{R}^d$, based on i.i.d. samples of the $\mathbb{R}^q$-valued random vector $U$.

Now assume, as it is the case of the Monte Carlo method, that the random vector $U$ is not directly simulatable (at a reasonable cost) but can only be approximated by another sequence of easily simulatable random vectors $((U^n)^p)_{p \ge 1}$, which strongly approximates $U$ as $n \to +\infty$ with a standard weak discretization error (or bias) $E[f(U)] - E[f(U^n)]$ of order $n^{-\alpha}$ for a specific class of functions. The computational cost required to simulate one sample of $U^n$ is assumed to be of order $n$ that is $Cost(U^n) = K \times n$ for some positive constant $K$. One standard situation corresponds to the case of a discretization of an SDE by means of an Euler-Maruyama scheme with $n$ time steps.

The important point here is that the function $h$ is generally neither known nor computable (at least at reasonable cost) and since the random variable $U$ cannot be simulated, estimating $\theta^*$ using a Robbins-Monro algorithm is not possible. Therefore, two steps are needed to compute $\theta^*$:

- the first step consists in approximating the zero $\theta^*$ of $h$ by the zero $\theta^{*,n}$ of $h^n$ defined by $h^n(\theta) := E[H(\theta, U^n)]$, $\theta \in \mathbb{R}^d$. It induces an implicit weak error which writes

$$ \mathcal{E}_D(n) := \theta^* - \theta^{*,n}. $$

Let us note that $\theta^{*,n}$ appears as a proxy of $\theta^*$ and one would naturally expect that $\theta^{*,n} \to \theta^*$ as the bias parameter $n$ tends to infinity.

- the second step consists in approximating $\theta^{*,n}$ by $M \in \mathbb{N}^*$ steps of the following SA scheme

$$ \theta_{p+1}^n = \theta_p^n - \gamma_{p+1} H(\theta_p^n, U^{n,(p+1)}), \quad 0 \le p \le M-1, \qquad (4.1) $$
where $(U^{n,(p)})_{1\le p\le M}$ is an i.i.d. sequence of random variables with the same law as $U^n$, $\theta_0^n$ is independent of the innovation of the algorithm with $\sup_{n\ge 1} \mathbb{E}[|\theta_0^n|^2] < +\infty$ and $\gamma = (\gamma_p)_{p\ge 1}$ is a sequence of non-negative deterministic and decreasing steps satisfying the usual step assumption (1.9). This induces a statistical error which writes

$$ \mathcal{E}_S(n, M, \gamma) := \theta^{*,n} - \theta_M^n. $$

The global error between $\theta^*$, the quantity to estimate, and its implementable approximation $\theta_M^n$ can be decomposed as follows:

$$ \begin{aligned} \mathcal{E}_{\text{glob}}(n, M, \gamma) &= \theta^* - \theta^{*,n} + \theta^{*,n} - \theta_M^n \\ &:= \mathcal{E}_D(n) + \mathcal{E}_S(n, M, \gamma). \end{aligned} $$

Concerning the asymptotic behaviour of the implicit weak error $\mathcal{E}_D(n)$, the following result is established in [9].

**Proposition 4.1.1 (Implicit discretisation error).** For all $n \in \mathbb{N}^*$, assume that $h$ and $h^n$ satisfy the mean reverting assumption (1.11) of Corollary 1.1.2. Moreover, suppose that $(h^n)_{n\ge 1}$ converges locally uniformly towards $h$. Then, one has

$$ \theta^{*,n} \to \theta^* \quad \text{as} \quad n \to +\infty. $$

Moreover, suppose that $h$ and $h^n$, $n \ge 1$, are continuously differentiable and that $Dh(\theta^*)$ is non-singular. Assume that $(Dh^n)_{n\ge 1}$ converges locally uniformly to $Dh$. If there exists $\alpha \in \mathbb{R}^*$ such that

$$ \forall \theta \in \mathbb{R}^d, \lim_{n \to +\infty} n^\alpha (h^n(\theta) - h(\theta)) = \mathcal{E}(h, \alpha, \theta), $$

then, one has

$$ \lim_{n \to +\infty} n^\alpha (\theta^{*,n} - \theta^*) = -Dh^{-1}(\theta^*) \mathcal{E}(h, \alpha, \theta^*). $$

Hence we see that the standard weak rate of convergence transfers to the implicit weak error under mild assumptions. The above proposition is the first step toward an expansion of $\theta^{*,n} - \theta^*$ in powers of $n^{-\alpha}$. This will be discussed in further details in Section 4.3 for the development of a multi-step Richardson-Romberg method.

Regarding the statistical error $\mathcal{E}_S(n, M, \gamma) := \theta^{*,n} - \theta_M^n$, it is well-known that under standard assumptions, the Robbins-Monro theorem, guarantees that $\lim_M \mathcal{E}_S(n, M, \gamma) = 0$ a.s. for each fixed $n \in \mathbb{N}^*$, see Corollary 1.1.2. Moreover, under mild technical conditions, a CLT holds at rate $\gamma^{-1/2}(M)$, that is, for each fixed $n \in \mathbb{N}^*$, $\gamma^{-1/2}(M)\mathcal{E}_S(n, M, \gamma)$ converges in distribution to a normally distributed random variable with mean zero and finite covariance matrix, see Theorem 1.1.5. In particular if we set $\gamma(p) = \gamma_0/p$, $\gamma_0 > 0$, $p \ge 1$, the weak convergence rate is $\sqrt{M}$ provided that $2\mathbb{R}\mathcal{E}(\lambda_{min})\gamma_0 > 1$ where $\lambda_{min}$ denotes the eigenvalue of $Dh(\theta^*)$ with the smallest real part. Moreover, the Ruppert and Polyak averaging principle allows one to achieve the optimal rate of convergence for free by devising a standard Robbins-Monro algorithm with step sequence $\gamma(p) = \gamma_0/p^\rho$, $\frac{1}{2} < \rho < 1$ and computing the empirical mean of the trajectory $(\theta_p^n)_{1\le p\le M}$.

Given the order of the implicit weak error and a step sequence $\gamma$ satisfying (1.9) a natural question is to find the optimal balance between the value of $n$ and the number $M$ of steps in the recursive scheme (4.1) in order to achieve a given global error. This problem was originally investigated by Duffie and Glynn [38] for the standard Monte Carlo method. The error between $\theta^*$ and the approximation $\theta_M^n$ writes $\theta_M^n - \theta^* = \theta_M^n - \theta^{*,n} + \theta^{*,n} - \theta^*$ suggesting to select $M = \gamma^{-1}(1/n^{2\alpha})$, where $\gamma^{-1}$ is the inverse function of $\gamma$, when the weak error is of order $n^{-\alpha}$.

However, due to the non-linearity of the SA algorithm (4.1), the methodology developed in [38] does not apply in our context. The key ingredients to tackle this question consists in linearising the dynamic of $(\theta_p^n)_{1\le p\le M}$ around its target $\theta^{*,n}$, quantifying the contribution of the non linearities in the space variable $\theta_p^n$ and the innovations and finally exploiting stability arguments from SA schemes. The next result is established in [9].

**Theorem 4.1.2 (Optimal tradeoff).** Suppose that the assumptions of Proposition 4.1.1 are satisfied and that $h$ satisfies the assumptions of Theorem 1.1.5. Assume that $h^n$ is twice continuously differentiable with $Dh^n$ Lipschitz continuous uniformly in $n$ and that it satisfies the strong mean-reverting assumption (1.14) with a constant $\lambda$ (instead of $\alpha$) that does not depend on $n$. Then, under some additional assumptions that we omit for sake of simplicity, one has

$$ n^\alpha (\theta_{\gamma^{-1}(1/n^{2\alpha})}^n - \theta^*) \implies -Dh^{-1}(\theta^*) \mathcal{E}(h, \alpha, \theta^*) + \mathcal{N}(0, \Sigma^*), $$

$$ \Sigma^* := \int_0^\infty \exp(-s(Dh(\theta^*) - \zeta I_d))^T \mathbb{E}[H(\theta^*, U)H(\theta^*, U)^T] \exp(-s(Dh(\theta^*) - \zeta I_d)) ds \quad (4.2) $$

$$ \text{with } \zeta = 0 \text{ if } \gamma(n) = \gamma_0/n^\rho, \quad \frac{1}{2} < \rho < 1 \text{ and } \zeta = 1/(2\gamma_0) \text{ if } \gamma(n) = \gamma_0/n, \quad \gamma_0 > 1/(2\lambda). $$
Define the empirical mean sequence $(\bar{\theta}_p^n)_{p \ge 1}$ of the sequence $(\theta_p^n)_{p \ge 1}$ by setting

$$ \bar{\theta}_p^n = \frac{\theta_0 + \theta_1^n + \dots + \theta_p^n}{p+1} = \bar{\theta}_{p-1}^n - \frac{1}{p+1} (\bar{\theta}_{p-1}^n - \theta_p^n), $$

with a the step sequence $\gamma_p = \gamma_0/p^\rho$ with $\rho \in (1/2, 1)$. Then, one has

$$ n^{\alpha} (\bar{\theta}_{n^{2\alpha}}^n - \theta^*) \implies -Dh^{-1}(\theta^*) \mathcal{E}(h, \alpha, \theta^*) + \mathcal{N}(0, Dh(\theta^*)^{-1}\mathbb{E}[H(\theta^*, U)H(\theta^*, U)^T](Dh(\theta^*)^{-1})^T), \quad n \to \infty. $$

Optimising with respect to the usual choice of the step sequence, the minimal computational cost (to achieve an error of order $n^{-\alpha}$) is given by $C_{SA} = K \times n \times \gamma^{-1}(1/n^{2\alpha})$ and is optimal for $\gamma(p) = \gamma_0/p$, $p \ge 1$, provided that $\gamma_0 > \frac{1}{2\alpha}$ or for the Ruppert and Polyak averaged algorithm, leading to a complexity of order $n^{2\alpha+1}$. In this case, we see that the computational cost is similar to the one achieved by the classical Monte Carlo algorithm for the computation of $\mathbb{E}[f(X_T)]$

In Section 4.2, we present a multi-level estimator which allows to significantly increase the efficiency the standard stochastic algorithm designed according to the optimal tradeoff.

## 4.2 Multi-level stochastic approximation algorithm

The multi-level stochastic approximation method uses the following decomposition $\theta^{*,n} = \theta^{*,1} + \sum_{\ell=1}^L (\theta^{*,m^\ell} - \theta^{*,m^{\ell-1}})$ with $m^L = n$, $m \ge 2$, and implements $L+1$ stochastic schemes with a sequence of bias parameter $(m^\ell)_{0 \le \ell \le L}$. More precisely, the target $\theta^*$ is estimated by the quantity

$$ \Theta_n^{ml} = \theta_{M_0}^1 + \sum_{\ell=1}^L \left( \theta_{M_\ell}^{m^\ell} - \theta_{M_\ell}^{m^{\ell-1}} \right). $$

It is important to point out here that for each level $\ell$ the couple $(\theta_{M_\ell}^{m^\ell}, \theta_{M_\ell}^{m^{\ell-1}})$ is computed using $M_\ell$ i.i.d. copies of $(U^{m^{\ell-1}}, U^{m^\ell})$. Moreover the random variables $U^{m^{\ell-1}}$ and $U^{m^\ell}$ use two different bias parameter but are perfectly correlated. Finally, for two different levels, the SA schemes are based on independent samples.

**Theorem 4.2.1.** Suppose that $h$ and $(h^n)_{n \in \mathbb{N}}$ satisfy the assumptions of Theorem 4.1.1. Assume that $h^n$ is twice continuously differentiable in a neighbourhood of $\theta^*$, with $Dh^n$ Lipschitz-continuous uniformly in $n$. Assume that the following CLT holds:

$$ m^{\ell\rho}(U^{m^\ell} - U^{m^{\ell-1}}) \stackrel{\text{stably}}{\implies} V^m \quad \text{as} \quad \ell \to +\infty $$

where $V^m$ is an $\mathbb{R}^q$-valued random variable eventually defined on an extension $(\tilde{\Omega}, \tilde{\mathcal{F}}, \tilde{\mathbb{P}})$ of $(\Omega, \mathcal{F}, \mathbb{P})$.

Suppose that $\tilde{\mathbb{E}}[(D_x H(\theta^*, U)V^m - \tilde{\mathbb{E}}[D_x H(\theta^*, U)V^m])(D_x H(\theta^*, U)V^m - \tilde{\mathbb{E}}[D_x H(\theta^*, U)V^m])^T]$ is a positive definite matrix. Assume that the step sequence is given by $\gamma_p = \gamma(p)$, $p \ge 1$, where $\gamma$ is a positive function defined on $[0, +\infty[$ decreasing to zero, satisfying one of the following assumptions:

* $\gamma$ varies regularly with exponent $(-a)$, $a \in (1/2, 1)$, that is, for any $x > 0$, $\lim_{t \to +\infty} \gamma(tx)/\gamma(t) = x^{-a}$.

* for $t \ge 1$, $\gamma(t) = \gamma_0/t$ and $\gamma_0$ satisfies $\lambda\gamma_0 > 1$.

Suppose that $\rho$ satisfies one of the following assumptions:

* if $\rho \in (0, 1/2)$, then assume that $\alpha > 2\rho$, $\lambda\gamma_0 > \alpha/(\alpha - 2\rho)$ (if $\gamma(t) = \gamma_0/t$) and

$$ \exists \beta > \rho, \quad \forall \theta \in \mathbb{R}^d, \quad \sup_{n \ge 1} n^\beta \|Dh^n(\theta) - Dh(\theta)\| < +\infty. $$

In this case we set:

$$ M_0 = \gamma^{-1}(1/(n^{2\alpha} \log(n))), \quad M_l = \gamma^{-1}(m^{\ell \frac{(1+2\rho)}{2}}(m^{\frac{1-2\rho}{2}} - 1)/(n^{2\alpha}(n^{\frac{(1-2\rho)}{2}} - 1))), \ell = 1, \dots, L. $$

* if $\rho = 1/2$, then assume that $\alpha = 1$, $\theta_0^{m^\ell} = \theta_0$, $\ell = 1, \dots, L$, with $\mathbb{E}[|\theta_0|^2] < +\infty$ and

$$ \exists \beta > 1/2, \quad \forall \theta \in \mathbb{R}^d, \quad \sup_{n \ge 1} n^\beta \|Dh^n(\theta) - Dh(\theta)\| < +\infty. $$

In this case we set:

$$ M_0 = \gamma^{-1}(1/(n^2 \log(n))), \quad M_l = \gamma^{-1}(m^l \log(m)/(n^2 \log(n)(m-1))), \ell = 1, \dots, L. $$
Then one has

$$n^{\alpha}(\Theta_{n}^{ml} - \theta^{*}) \implies -Dh^{-1}(\theta^{*})\mathcal{E}(h, 1, \theta^{*}) + \mathcal{N}(0, \Sigma^{*}), \quad n \to +\infty$$

with

$$\Sigma^* := \int_0^\infty \left( e^{-s(Dh(\theta^*) - \zeta I_d)} \right)^T \tilde{\mathbb{E}} \left[ \left( D_x H(\theta^*, U) V - \tilde{\mathbb{E}}[D_x H(\theta^*, U) V] \right) \left( D_x H(\theta^*, U) V - \tilde{\mathbb{E}}[D_x H(\theta^*, U) V] \right)^T \right] e^{-s(Dh(\theta^*) - \zeta I_d)} ds.$$

**Remark 4.2.1.** The previous CLT for the multi-level stochastic approximation estimator of $\theta^*$ holds if the standard weak error (and thus the implicit weak error), is of order $1/n^\alpha$ and the strong rate error is of order $1/n^\rho$ with $\alpha > 2\rho$ or $\alpha = 1$ and $\rho = 1/2$. Due to the non-linearity of the scheme, which leads to annoying remainder terms in the Taylor's expansions used in the proof, those results do not seem to easily extend to a weak discretisation error of order $1/n^\alpha$ with $\alpha < 1$ and $\rho = 1/2$ or a faster strong convergence rate $\rho > 1/2$. Moreover, for the same reason this result does not seem to extend to the empirical sequence associated to the multi-level estimator according to the Ruppert & Polyak averaging principle.

The assumption related to the CLTs satisfied by the two sequences $(U^{(n)} - U)_{n \ge 1}$ and $(U^{m^\ell} - U^{m^{\ell-1}})_{\ell \ge 1}$ includes the case of the value at time $T$ of a SDE, namely $U = X_T$ approximated by its continuous Euler-Maruyama scheme $U^n = X_T^n$ with $n$ steps. Under mild assumption on the coefficients of the SDE, one has $\rho = 1/2$. Moreover, $U$ may depend on the whole path of the process. For instance, one may have $U = L_T^0(X)$, the symmetric local time at level 0 of the one-dimensional continuous diffusion process $X$ and the approximations may be given by

$$U^n = \sum_{i=1}^{[nt]} f \left( u_n X_{\frac{i-1}{n}}, \sqrt{n} \left( X_{\frac{i}{n}} - X_{\frac{i-1}{n}} \right) \right).$$

Then under some assumptions on the function $f$ and the coefficients $b$, $\sigma$, the weak and strong rate of convergence are $\rho = 1/4$, see Jacod [65] for more details. Let us note that we do not know what happens in our results when $\rho > 1/2$ which includes the case of higher order schemes for discretisation schemes of SDE, like the Milstein scheme.

The CLT presented in Theorem 4.2.1 shows that in order to obtain an error of order $1/n^\alpha$, one has to set $M_0 = \gamma^{-1}(1/n^{2\alpha})$ and $M_l = \gamma^{-1}(m^{\ell(1+2\rho)/2}(m^{\frac{1-2\rho}{2}} - 1)/(n^{2\alpha}(n^{\frac{1-2\rho}{2}} - 1)))$, if $\rho \in (0, 1/2)$ or $M_0 = \gamma^{-1}(1/n^2)$ and $M_l = \gamma^{-1}(m^\ell \log(m)/(n^2 \log(n)(m-1)))$ if $\rho = 1/2$, $\ell = 1, \dots, L$ with $L = \log(n)/\log(m)$. In both cases the complexity of the multi-level stochastic approximation method is given by

$$C_{\text{ML-SA}}(\gamma) = C \times \left( \gamma^{-1}(1/(n^{2\alpha} \log(n))) + \sum_{\ell=1}^{L} M_{\ell}(m^{\ell} + m^{\ell-1}) \right). \quad (4.3)$$

In order to compute this complexity, we distinguish the two following cases:

* If $\gamma(p) = \gamma_0/p$ then the optimal complexity is given by

$$C_{\text{ML-SA}}(\gamma) = C \left( n^{2\alpha} \log(n) + \frac{n^2(n^{\frac{1-2\rho}{2}} - 1)}{m^{\frac{1-2\rho}{2}} - 1} \sum_{\ell=1}^{L} m^{-\frac{(1+2\rho)}{2}\ell} (m^{\ell} + m^{\ell-1}) \right) = \mathcal{O}(n^{2\alpha} n^{1-2\rho}),$$

if $\rho \in (0, 1/2)$ under the constraint $\lambda\gamma_0 > \alpha(\alpha - 2\rho)$ and

$$C_{\text{ML-SA}}(\gamma) = C \left( n^2 \log(n) + n^2 (\log n)^2 \frac{m^2 - 1}{m(\log m)^2} \right) = \mathcal{O}(n^2 (\log(n))^2),$$

if $\rho = 1/2$ under the constraint $\lambda\gamma_0 > 1$. These computational costs are similar to those achieved by the multi-level Monte Carlo method for the computation of $\mathbb{E}[f(X_T)]$, see Giles [48] and Ben Alaya and Kebaier [1]. As discussed in Giles [48], this complexity attains a minimum near $m = 7$.

* If we choose $\gamma(p) = \gamma_0/p^a$, $\frac{1}{2} < a < 1$ then simple computations show that the computational cost is given by

$$C_{\text{ML-SA}}(\gamma) = C \left( n^{\frac{2\alpha}{a}} \log^{\frac{1}{a}}(n) + n^{\frac{2}{a}}(n^{1-2\rho} - 1)^{\frac{1}{a}} \sum_{\ell=1}^{L} m^{-\frac{(1+2\rho)}{a}\ell} (m^{\ell} + m^{\ell-1}) \right) = \mathcal{O}(n^{\frac{2\alpha}{a}} n^{\frac{1-2\rho}{a}}),$$

if $\rho \in (0, 1/2)$ and

$$C_{\text{ML-SA}}(\gamma) = C \times \left( n^{\frac{2}{a}} \log^{\frac{1}{a}}(n) + n^{\frac{2}{a}} (\log n)^{\frac{1}{a}} \frac{(m-1)^{\frac{1}{a}}(m+1)}{m(\log m)^{\frac{1}{a}}} \sum_{\ell=1}^{L} m^{-\ell(\frac{1}{a}-1)} \right) = \mathcal{O}(n^{\frac{2}{a}} (\log n)^{\frac{1}{a}})$$

if $\rho = 1/2$.
Hence we see that the optimal complexity is achieved for $\gamma(p) = \gamma_0/p$ under a constraint on $\gamma_0$ involving $\lambda$. We also obtained a CLT for the statistical Romberg stochastic approximation estimator that is defined by $\Theta_n^{sr} = \theta_{M_1}^{n\beta} + \theta_{M_2}^{n\beta} - \theta_{M_2}^{n\beta}$. The couple $(\theta_{M_2}^n, \theta_{M_2}^{n\beta})$ is computed using i.i.d. copies of $(U^{n\beta}, U^n)$, the random variables $U^{n\beta}$ and $U^n$ being perfectly correlated. Moreover, the random variables used to obtain $\theta_{M_1}^{n\beta}$ are independent to those used for the computation of $(\theta_{M_2}^n, \theta_{M_2}^{n\beta})$. The optimal parameters are $M_1 = \gamma^{-1}(1/n^{2\alpha})$ and $M_2 = \gamma^{-1}(1/(n^{2\alpha-2\rho\beta}))$. When $\gamma(p) = \gamma_0/p$, a constraint on $\gamma_0$ involving $\lambda$ appears. Here we can circumvent this issue by considering the averaged version of the last estimator namely $\Theta_n^{sr} = \bar{\theta}_{M_3}^{n\beta} + \bar{\theta}_{M_4}^{n\beta} - \bar{\theta}_{M_4}^{n\beta}$, with $M_3 = n^{2\alpha}$ and $M_4 = n^{2\alpha-2\rho\beta}$, for any $\beta \in (0, 1)$ (note that $\rho \le \alpha$ so that $\rho\beta < \alpha$). The optimal choice is $\beta^* = 1/(1+2\rho)$ leading to an optimal complexity $C_{SR-SA}(\gamma) = C_{SR-RP}(\gamma) = C'n^{2\alpha+\frac{1}{1+2\rho}}$.

As a numerical illustration of these results we consider the following toy example. We turn our attention to the computation of the level of the function $\theta \mapsto e^{-rT}\mathbb{E}(X_T - \theta)_+$ (European call option), $X_T = x \exp((r - \frac{\sigma^2}{2})T + \sigma W_T)$, which is given by the Black-Scholes formula.

Therefore, we first fix a value $\theta^*$ (the target of our procedure) and compute the corresponding level $\ell = e^{-rT}\mathbb{E}(X_T - \theta^*)_+$ by the Black-Scholes formula. Accordingly, we set $h(\theta) = \mathbb{E}[H(\theta, X_T)]$ with $H(\theta, x) = \ell - (x - \theta)_+$. We use the following values for the parameters: $x = 100$, $r = 0.05$, $\sigma = 0.4$, $T = 1$. We plot in Figure 4.1 the behaviors of $nh^n(\theta^*)$ and $n(\theta^{*,n} - \theta^*)$ for $n = 100, \dots, 500$. The value of $h^n(\theta^*)$ is approximated by its Monte Carlo estimator and $\theta^{*,n}$ is estimated by $\theta_M^n$, both estimators being computed with $M = 10^8$ samples. The variance of the Monte Carlo estimator ranges from $9.73 \times 10^6$ for $n = 100$ to $9.39 \times 10^7$ for $n = 500$.

To compare the three methods in terms of computational costs, we compute the different estimators, namely $\theta_{\gamma^{-1}(1/n^2)}^n$ where $(\theta_p^n)_{p \ge 1}$ is given by (4.1), $\Theta_n^{sr}$ and $\Theta_n^{ml}$ for a set of $N = 200$ values of the target $\theta^*$ equidistributed on the interval $[90, 110]$ and for different values of $n$. For each value $n$ and for each method we compute its complexity as the number of operations and the root-mean-squared error which is given by

$$ \text{RMSE} = \left( \frac{1}{N} \sum_{k=1}^{N} (\Theta_k^n - \theta_k^{*n})^2 \right)^{1/2} $$

where $\Theta_k^n = \theta_{\gamma^{-1}(1/n^2)}^n$, $\Theta_n^{sr}$ or $\Theta_n^{ml}$ is the considered estimator. For each given $n$, we provide a couple (RMSE, Complexity) which is plotted on Figure 4.2. Let us note that the multi-level SA estimator has been computed for different values of $m$ (ranging from $m = 2$ to $m = 7$) and different values of $L$. We set $\gamma(p) = \gamma_0/p$, with $\gamma_0 = 2$, $p \ge 1$, so that $\beta^* = 1/2$.

We briefly comment Figure 4.2. The curves of the statical Romberg SA and the multi-level SA methods are displaced below the curve of the SA method. Therefore, for a given error, the complexity of both methods are much lower than the one of the crude SA. The difference in terms of computational cost becomes more significant as the RMSE is small, which corresponds to large values of $n$. The difference between the statistical Romberg and the multi-level SA method is not significant for small values of $n$, i.e. for a RMSE between 1 and 0.1. For a RMSE lower than $5.10^{-2}$, which corresponds to a number of steps $n$ greater than about 600-700, we observe that the multi-level SA procedure becomes much more effective than both methods. For a RMSE fixed around 1 (which corresponds to $n = 100$ for the SA algorithm and Statical Romberg SA), one divides the complexity by a factor of approximately 5 by using the statistical romberg SA. For a RMSE fixed at $10^{-1}$, the computational cost gain is approximately equal to 10 by using either the statistical romberg SA algorithm or the multi-level SA one. Finally, for a RMSE fixed at $5.5 \times 10^{-2}$, the complexity gain achieved by using the multi-level SA procedure instead of the statistical Romberg one is approximately equal to 5.

## 4.3 Multi-step Richardson-Romberg for stochastic approximation algorithm

This section is devoted to the work [10] related to a Multi-step Richardson-Romberg method for the stochastic approximation algorithms introduced in [9]. Our first result is an expansion of the implicit discretisation error $\theta^{*,n} - \theta^*$ presented in the previous expansion. We will denote by $[\mathbf{H}-\mathbf{k}]$ the following set of assumptions:

1. For all $\theta \in \mathbb{R}^d$,

$$ h(\theta) - h^n(\theta) = \frac{\Lambda_1^n(\theta)}{n^\alpha} + \dots + \frac{\Lambda_k^n(\theta)}{n^{\alpha k}} + o\left(\frac{1}{n^{\alpha k}}\right). \quad (4.4) $$

2. $h, h^n \in C^k(\mathbb{R}^d, \mathbb{R}^d)$ and for all $l \le k - 1$, for all $\theta \in \mathbb{R}^d$,

$$ D^l h^n(\theta) - D^l h(\theta) = \frac{\Lambda_1^l(\theta)}{n^\alpha} + \dots + \frac{\Lambda_{k-l}^l(\theta)}{n^{\alpha(k-l)}} + o\left(\frac{1}{n^{\alpha(k-l)}}\right). \quad (4.5) $$