# Particle Method for the McKean-Vlasov equation with common noise

Théophile Le Gall

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

To cite this version:

Théophile Le Gall. Particle Method for the McKean-Vlasov equation with common noise. 2024. [hal-04851158]

HAL Id: hal-04851158

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

Preprint submitted on 20 Dec 2024

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Particle method for the McKean-Vlasov equation with common noise

Théophile Le Gall*

December 20, 2024

## Abstract

This paper studies the numerical simulation of the solution to the McKean-Vlasov equation with common noise. We begin by discretizing the solution in time using the Euler scheme, followed by spatial discretization through the particle method, inspired by the propagation of chaos property. Assuming Hölder continuity in time, as well as Lipschitz continuity in the state and measure arguments of the coefficient functions $b, \sigma$ and $\sigma^0$, we establish the convergence rate of the Euler scheme and the particle method. These results extend those in [24] for the standard McKean-Vlasov equation without common noise. Finally, we present two simulation examples: a modified conditional Ornstein Uhlenbeck process with common noise and an interbank market model presented in [29].

**Keywords.** Euler scheme, McKean-Vlasov equation with common noise, Mean-field limits, Numerical analysis of the particle method.

## 1 Introduction

We consider the $\mathbb{R}^d$-valued McKean-Vlasov stochastic differential equation (SDE) with common noise defined for $t \in [0, T]$ by

$$dX_t = b(t, X_t, \mathcal{L}^1(X_t))dt + \sigma(t, X_t, \mathcal{L}^1(X_t))dW_t + \sigma^0(t, X_t, \mathcal{L}^1(X_t))dW_t^0, \quad (1.1)$$

where, for some $T > 0$, $W = (W_t)_{t \in [0,T]}$ and $W^0 = (W_t^0)_{t \in [0,T]}$ are two independent Brownian motions respectively called idiosyncratic noise and common noise. The coefficient $b$ is a mapping from $[0, T] \times \mathbb{R}^d \times \mathcal{P}_p(\mathbb{R}^d)$ to $\mathbb{R}^d$ where $\mathcal{P}_p(\mathbb{R}^d)$ is the space of probability measure having a $p$-th finite moment for $p \ge 1$. We endow $\mathcal{P}_p(\mathbb{R}^d)$ with the $p$-Wasserstein metric defined later in (1.3). The coefficients $\sigma, \sigma^0$ are mappings from $[0, T] \times \mathbb{R}^d \times \mathcal{P}_p(\mathbb{R}^d)$ to $\mathbb{R}^d \times \mathbb{R}^q$ which represent respectively the intensity of the idiosyncratic noise $W$ and the common noise $W^0$. The notation $\mathcal{L}^1(\cdot)$ represents for the conditional law given the trajectory of the common noise (see further (1.2), Proposition 1.1 and 3.3). The initial condition of Equation (1.1) is a random variable $X_0$ independent of $W$ and $W^0$.

The McKean-Vlasov equation, initially introduced by H. McKean [26] is a nonlinear partial differential equation (PDE) associated with a class of stochastic differential equations (SDEs) where the drift and diffusion coefficients depend not only on the time and the state of the process, but also on its marginal laws (see e.g. [30]). The distribution-dependent structure of the McKean-Vlasov equation is extensively applied for modeling phenomena in statistical physics (see e.g. [5], [25]), mathematical biology (see e.g. [3]), social sciences and quantitative finance, both frequently driven by advancements in mean field

<footnote>*CEREMADE, CNRS, UMR 7534, Université Paris-Dauphine, PSL University, 75016 Paris, France, legall@ceremade.dauphine.fr.</footnote>
games and interacting diffusion models (see e.g. [7], [23]). The McKean-Vlasov equation equipped with a common noise (see further (1.1) for a precise definition) was first introduced in [1], [10], [16] and [22], where the term common noise served to model a type of shared risk in a particle system. Several papers such as [13] or [14] explore how the introduction of the common noise can restore uniqueness in mean-field games, which are derived from deterministic differential games involving a large number of players.

This paper aims to develop a numerical method for the McKean-Vlasov equation with common noise, accompanied by an analysis of the associated convergence rate. For the standard McKean-Vlasov equation without common noise, under Lipschitz continuity assumptions on the coefficient functions, we refer to [6] and [24] for the simulation of the solution to the SDE, and to [2] and [17] for the estimation of the density solution to the McKean-Vlasov PDE. Additionally, recent advancements in handling super-linear growth coefficient functions can be found in [28], [12] among others. For the simulation of the invariant measure, we refer to [11]. For the McKean-Vlasov equation with common noise, a recent study [4] provides the convergence rate of a numerical scheme in a different setting from the one considered in this paper. For a detailed comparison, see further Remark 1.4.

Following the construction as presented in [6] and [24] for standard McKean-Vlasov equation without common noise, our approach employs the Euler scheme, defined further in (1.4), as a temporal discretization, and the particle method as a spatial discretization, defined further in (1.6), that was introduced in [24] and we extend it to account for the case with common noise. Notice that the addition of common noise needs further careful consideration, particularly in terms of the conditional distributions given the common noise in the measure argument of the coefficient functions. In this context, the empirical measure serves as an estimator for the law of the solution process, conditioned on the common noise. This method relies on the propagation of chaos property, initially introduced by Kac [18] and further studied in [21], [19], and [30].

## 1.1 Probabilistic settings

In this paper, we consider two filtered probability spaces $(\Omega^0, \mathcal{F}^0, \mathbb{F}^0, \mathbb{P}^0)$ and $(\Omega^1, \mathcal{F}^1, \mathbb{F}^1, \mathbb{P}^1)$ satisfying the usual conditions, where $\mathbb{F}^0 = (\mathcal{F}_t^0)_{t \in [0,T]}$ and $\mathbb{F}^1 = (\mathcal{F}_t^1)_{t \in [0,T]}$. In addition, we provide $W^0 = (W_t^0)_{t \in [0,T]}$ and $W = (W_t)_{t \in [0,T]}$ two $q$-dimensional $\mathbb{F}^0$-adapted and $\mathbb{F}^1$-adapted Wiener processes respectively supported on $(\Omega^0, \mathcal{F}^0, \mathbb{P}^0)$ and $(\Omega^1, \mathcal{F}^1, \mathbb{P}^1)$ which respectively represent the common noise and the idiosyncratic noise. We introduce naturally the product space $(\Omega, \mathcal{F}, \mathbb{F}, \mathbb{P})$ where $\Omega = \Omega^0 \times \Omega^1$, $(\mathcal{F}, \mathbb{P})$ is the completion of $(\mathcal{F}^0 \otimes \mathcal{F}^1, \mathbb{P}^0 \otimes \mathbb{P}^1)$ and $\mathbb{F} = (\mathcal{F}_t)_{t \in [0,T]}$ is the complete and right-continuous augmentation of $(\mathcal{F}_t^0 \otimes \mathcal{F}_t^1)_{t \in [0,T]}$. Consider a random variable $X$ defined on the filtered probability space $(\Omega, \mathcal{F}, \mathbb{F}, \mathbb{P})$. Then, for $\mathbb{P}^0$-a.e., $\omega^0 \in \Omega^0$, $X(\omega^0, \cdot)$ is a random variable on $(\Omega^1, \mathcal{F}^1, \mathbb{P}^1)$ (see e.g. [9, Section 2.1.3]). In particular, we may define

$$ \mathcal{L}^1(X) : \omega^0 \in \Omega^0 \mapsto \mathcal{L}(X(\omega^0, \cdot)) \in \mathcal{P}(\mathbb{R}^d), \quad (1.2) $$

for almost every $\omega^0 \in \Omega^0$. On the exceptional event where $\mathcal{L}(X(\omega^0, \cdot))$ cannot be computed, we may assign it arbitrary values in $\mathcal{P}(\mathbb{R}^d)$.

**Proposition 1.1** ([9, Lemma 2.4]). Given a random variable $X : (\Omega, \mathcal{F}, \mathbb{P}) \to (\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d))$, the mapping $\mathcal{L}^1(X)$ defined by (1.2) is almost surely well defined under $\mathbb{P}^0$, and forms a random variable from $(\Omega^0, \mathcal{F}^0, \mathbb{P}^0)$ into $\mathcal{P}(\mathbb{R}^d)$ endowed with its Borel $\sigma$-field generated by the Lévy-Prokhorov metric (see [8, Section 5.1.1 and Proposition 5.7]). Moreover, the random variable $\mathcal{L}^1(X)$ provides a conditional law of $X$ given $\mathcal{F}^0$.

We endow $\mathcal{P}_p(\mathbb{R}^d)$ with the $p$-Wasserstein metric

$$ W_p(\mu, \nu) = \inf_{\pi \in \Pi(\mu, \nu)} \left[ \int_{\mathbb{R}^d \times \mathbb{R}^d} |x - y|^p \pi(dx, dy) \right]^{\frac{1}{p}}, \quad (1.3) $$
where $\Pi(\mu, \nu)$ denotes the set of probability measure on $\mathbb{R}^d \times \mathbb{R}^d$ with respective marginals $\mu$ and $\nu$.

## 1.2 Construction of the Euler scheme and particle method

Let $M$ be the number of time discretization and $h = \frac{T}{M}$ be the time step. For every $m \in \{0, \dots, M\}$, we define $t_m = hm$. We consider $W^1, \dots, W^N$ i.i.d. copies of the Brownian motion $W$, and define the re-normalized increments $Z_m^i, Z_m, Z_m^0, 1 \le m \le M, 1 \le i \le N$ as follows

$$Z_{m+1}^i = \frac{1}{\sqrt{h}}(W_{t_{m+1}}^i - W_{t_m}^i), Z_{m+1}^0 = \frac{1}{\sqrt{h}}(W_{t_{m+1}}^0 - W_{t_m}^0), Z_{m+1} = \frac{1}{\sqrt{h}}(W_{t_{m+1}} - W_{t_m}).$$

The theoretical Euler scheme of the McKean-Vlasov equation with common noise (1.1) is defined, for $0 \le m \le M-1$, by $\bar{X}_0 = X_0$ and

$$\bar{X}_{t_{m+1}}^M = \bar{X}_{t_m}^M + h \cdot b(t_m, \bar{X}_{t_m}^M, \mathcal{L}^1(\bar{X}_{t_m}^M)) + \sqrt{h} \cdot \sigma(t_m, \bar{X}_{t_m}^M, \mathcal{L}^1(\bar{X}_{t_m}^M)) Z_{m+1} + \sqrt{h} \cdot \sigma^0(t_m, \bar{X}_{t_m}^M, \mathcal{L}^1(\bar{X}_{t_m}^M)) Z_{m+1}^0, \quad (1.4)$$

equipped with its natural continuous extension defined, for every $t \in [t_m, t_{m+1})$, by

$$\bar{X}_t^M = \bar{X}_{t_m}^M + b(t_m, \bar{X}_{t_m}^M, \mathcal{L}^1(\bar{X}_{t_m}^M))(t-t_m) + \sigma(t_m, \bar{X}_{t_m}^M, \mathcal{L}^1(\bar{X}_{t_m}^M))(W_t - W_{t_m}) + \sigma^0(t_m, \bar{X}_{t_m}^M, \mathcal{L}^1(\bar{X}_{t_m}^M))(W_t^0 - W_{t_m}^0). \quad (1.5)$$

At each time step $t_m$, we build an $N$-particle system $(\bar{X}_{t_m}^{1,N}, \dots, \bar{X}_{t_m}^{N,N})$ such that for every $i \in \{1, \dots, N\}$, we have

$$\left\{ \begin{array}{l} \bar{X}_{t_{m+1}}^{i,N} = \bar{X}_{t_m}^{i,N} + h \cdot b(t_m, \bar{X}_{t_m}^{i,N}, \bar{\mu}_{t_m}^N) + \sqrt{h} \cdot \sigma(t_m, \bar{X}_{t_m}^{i,N}, \bar{\mu}_{t_m}^N) Z_{m+1}^i \\ \qquad + \sqrt{h} \cdot \sigma^0(t_m, \bar{X}_{t_m}^{i,N}, \bar{\mu}_{t_m}^N) Z_{m+1}^0 \\ \bar{\mu}_{t_m}^N = \frac{1}{N} \sum_{i=1}^N \delta_{\bar{X}_{t_m}^{i,N}}. \end{array} \right. \quad (1.6)$$

also equipped with its natural continuous extension defined, for every $t \in [t_m, t_m + 1)$, by

$$\bar{X}_t^{i,N} = \bar{X}_{t_m}^{i,N} + \int_{t_m}^{t} b(t_m, \bar{X}_{t_m}^{i,N}, \bar{\mu}_{t_m}^N) ds + \int_{t_m}^{t} \sigma(t_m, \bar{X}_{t_m}^{i,N}, \bar{\mu}_{t_m}^N) dW_s^i + \int_{t_m}^{t} \sigma^0(t_m, \bar{X}_{t_m}^{i,N}, \bar{\mu}_{t_m}^N) dW_s^0. \quad (1.7)$$

In the system (1.6), at each time step $t_m$, the particles $\bar{X}_{t_m}^{i,N}, 1 \le i \le N$, have interaction through the empirical measure $\bar{\mu}_{t_m}^N$. The idea of the particle method is to use $\bar{\mu}_{t_m}^N$ as an estimator of $\mathcal{L}^1(\bar{X}_{t_m}^M)$ in definition (1.4), at each time step $t_m, 0 \le m \le M$.

## 1.3 Assumptions and main results

The main results of this paper will be established under the following assumptions, which are assumed to be held for a fixed $p \in [2, \infty)$.

**Assumption 1.** The random variable $X_0$ is defined on $(\Omega^1, \mathcal{F}^1, \mathbb{P}^1)$ such that $\mathbb{E}^1(|X_0|^p) < \infty$.

**Assumption 2.** The coefficient mappings $b, \sigma$ and $\sigma^0$ are continuous in time and Lipschitz continuous in the state and measure arguments, that is, there exists a constant $L > 0$ such that for every $t \in [0, T]$, $x, y \in \mathbb{R}^d$ and $\mu, \nu \in \mathcal{P}_p(\mathbb{R}^d)$, we have

$$\max(|b(t, x, \mu) - b(t, y, \nu)|, |\sigma(t, x, \mu) - \sigma(t, y, \nu)|, |\sigma^0(t, x, \mu) - \sigma^0(t, y, \nu)|) \le L(|x - y| + W_p(\mu, \nu)).$$
Assumptions 1 and 2 guarantee the existence and strong uniqueness of a solution $X = (X_t)_{t \in [0,T]}$ to the McKean-Vlasov equation with common noise (1.1) satisfying the following estimate

$$ \left\| \sup_{t \in [0,T]} |X_t| \right\|_p \le C (1 + \|X_0\|_p), \quad (1.8) $$

where $C$ is a positive constant depending on $p, T, L, b, \sigma$ and $\sigma^0$. For the proof in the case $p = 2$, we refer to Proposition 2.8 in [9]. The proof for $p > 2$ follows a similar approach, with only minor differences.

**Assumption 3.** The coefficient mappings $b, \sigma$ and $\sigma^0$ are $\rho$-Hölder continuous in time, for some $\rho \in (0, 1]$ uniformly in space and measure, in the sense that there exists a constant $L > 0$ such that for every $s, t \in [0, T]$, $x \in \mathbb{R}^d$, $\mu \in \mathcal{P}_p(\mathbb{R}^d)$, we have

$$ \max(|b(t, x, \mu) - b(s, x, \mu)|, |\sigma(t, x, \mu) - \sigma(s, x, \mu)|, |\sigma^0(t, x, \mu) - \sigma^0(s, x, \mu)|) \le L(1 + |x| + W_p(\mu, \delta_0)) |t - s|^\rho. $$

The main results of this paper are the following two theorems whose proofs are presented in Section 4.

**Theorem 1.2.** Assume that Assumptions 1, 2 hold for some $p \in [2, \infty)$. Fix $M \in \mathbb{N}$ and set $h = \frac{T}{M}$. For every $m = 0, \dots, M$, let $\bar{\mu}_{t_m} = \mathcal{L}^1(\bar{X}_{t_m})$, where $\bar{X}_{t_m}$ is defined by (1.4) and let $(\bar{\mu}_{t_m})_{1 \le m \le M}$ denote the empirical measure of the particles $(\bar{X}_{t_m}^{i,N})_{1 \le i \le N}$ defined by the particle method (1.6).

(i) We have

$$ \sup_{m \in \{1, \dots, M\}} \left\| W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m}) \right\|_p \xrightarrow[N \to +\infty]{} 0. $$

(ii) Moreover, if we assume that Assumption 1 holds for $p + \varepsilon \in (2, \infty)$ for some $\varepsilon > 0$, we have the following rates of convergence

$$ \sup_{m \in \{1, \dots, M\}} \left\| W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m}) \right\|_p \le C \begin{cases} N^{-\frac{1}{2p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > \frac{d}{2} \text{ and } \varepsilon \ne p, \\ N^{-\frac{1}{2p}} (\log(1+N))^{\frac{1}{p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = \frac{d}{2} \text{ and } \varepsilon \ne p, \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, \frac{d}{2}) \text{ and } p + \varepsilon \ne \frac{d}{d-p}, \end{cases} $$

where $C$ is a positive constant which depends on $d, p, L, T, \|X_0\|_{p+\varepsilon}, b, \sigma$ and $\sigma^0$.

**Theorem 1.3.** Consider $(X_t)_{t \in [0,T]}$ the solution to the McKean-Vlasov equation with common noise (1.1) and let $(\bar{X}_t^{1,N})_{t \in [0,T]}$ be the process defined by the particle method (1.7) which depends on the same Brownian motions as $(X_t)_{t \in [0,T]}$. Then there exists a positive constant $C$ which depends on $d, p, L, T, \|X_0\|_{p+\varepsilon}, b, \sigma, \sigma^0$, such that

$$ \left\| \sup_{t \in [0,T]} |\bar{X}_t^{1,N} - X_t| \right\|_p \le C(h^{\frac{1}{2}} \wedge \rho + \mathcal{E}_N), \quad (1.9) $$

where

$$ \mathcal{E}_N = \begin{cases} N^{-\frac{1}{2p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > \frac{d}{2} \text{ and } \varepsilon \ne p, \\ N^{-\frac{1}{2p}} (\log(1+N))^{\frac{1}{p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = \frac{d}{2} \text{ and } \varepsilon \ne p, \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, \frac{d}{2}) \text{ and } p + \varepsilon \ne \frac{d}{d-p}. \end{cases} $$
**Remark 1.4.** We highlight the differences between this paper and [4]. In [4], the authors also analyzed the convergence rate of the particle method for the McKean-Vlasov equation with common noise. The key distinctions can be summarized in two aspects: (i) the construction of the numerical approaches, and (ii) the assumptions imposed on the initial random variable and the coefficient functions, along with the resulting convergence rate for the time discretization.

Regarding the first aspect, the approach in [4] begins with the particle system used in the conditional propagation of chaos property for the McKean-Vlasov equation (see e.g. [9, Section 2.1.4]), which is defined by

$$dX_t^{n,N} = b(t, X_t^{n,N}, \frac{1}{N}\sum_{i=1}^N \delta_{X_t^{i,N}})dt + \sigma(t, X_t^{n,N}, \frac{1}{N}\sum_{i=1}^N \delta_{X_t^{i,N}})dW_t^n + \sigma^0(t, X_t^{n,N}, \frac{1}{N}\sum_{i=1}^N \delta_{X_t^{i,N}})dW_t^0 \\ \text{with } X_0^{1,N}, \dots, X_0^{1,N} \stackrel{\text{i.i.d.}}{\sim} X_0, W^1, \dots, W^N \stackrel{\text{i.i.d.}}{\sim} W, \qquad (1.10)$$

and subsequently applies a time discretization by using a Milstein-type scheme to this particle system. It is worth noting that (1.10) can be regarded as a high-dimensional equation in terms of state arguments, without involving the measure argument. This feature enables the application of classical numerical analysis methods for diffusion processes to study the system (1.10). In our paper, we first apply time discretization using the Euler scheme, retaining the measure argument within the scheme. This approach allows for the potential integration of other spatial discretization methods in future work, such as the optimal quantization method, as discussed in [24] for the standard McKean-Vlasov equation without common noise.

As for the second aspect, the approach in [4] imposes additional regularity conditions on the coefficient functions $\sigma$ and $\sigma^0$ with respect to both the state and measure arguments. Specifically, it requires:

$$|\partial_x \sigma_{\ell}^u(t, x, \mu) \sigma_{\ell}^v(t, x, \mu) - \partial_x \sigma_{\ell}^u(t, x', \mu') \sigma_{\ell}^v(t, x', \mu')| \le L\{|x - x'| + W_2(\mu, \mu')\}, \\ |\partial_{\mu} \sigma_{\ell}^u(t, x, \mu, y) \sigma_{\ell}^v(t, y, \mu) - \partial_{\mu} \sigma_{\ell}^u(t, x', \mu', y') \sigma_{\ell}^v(t, y', \mu')| \le L\{|x - x'| + |y - y'| + W_2(\mu, \mu')\},$$

for all $u, v \in \{0, 1\}, \ell \in \{1, \dots, m_u\}, \ell_1 \in \{1, \dots, m_v\}, t \in [0, T], x, x', y, y' \in \mathbb{R}^d$ and $\mu, \mu' \in \mathcal{P}_2(\mathbb{R}^d)$, where $\sigma^1$ in their paper corresponds to $\sigma$ here. Additionally, they assume that $X_0$ has a finite p-th moment with $p \ge 4$. These conditions enable a faster convergence rate with respect to the time step h. However, they exclude certain coefficient functions, such as $x \mapsto \sigma(x) = |x|$ or $x \in \mathbb{R} \mapsto \sigma(x) = \sqrt{x} \mathbf{1}_{[0,1]}(x) + \mathbf{1}_{(1,\infty]}(x)$, which can be handled within the framework proposed in this paper.

## 1.4 Organization of the paper

The paper is organized as follows. Some preliminary results are gathered in Section 2 along with some notations. Sections 3 and 4 respectively present the proofs for the convergence rates of the Euler scheme (see further Proposition 3.1) and the particle method (Theorem 1.2, Theorem 1.3). Section 5 provides numerical examples to illustrate the methods discussed in this paper. The first example is a modified Ornstein-Uhlenbeck process to which Brownian common noise has been added. For the second example, we simulate the Interbank market model presented in [29][Section 5] which is an application of a risk-sensitive mean field games with common noise. Appendix A is dedicated to presenting the detailed proofs of the lemmas referenced throughout the paper, that are essential to supporting the proofs of the main results.

## 2 Preliminary results

In this paper, we fix a terminal time $T > 0$, and denote the space of continuous function from $[0, T]$ to a Polish space $S$ by $C([0, T], S)$. We also use the notation $L^p([0, T] \times \Omega)$ for the set of $(\mathcal{F}_t)_{t \in [0, T]}$-progressively
measurable continuous processes $X = (X_t)_{t \in [0,T]}$ such that

$$ \|X\|_{L^p([0,T] \times \Omega)} = \left\| \sup_{t \in [0,T]} |X_t| \right\|_p < \infty. $$

We now list some key lemmas that will support the subsequent proofs.

**Lemma 2.1** (Lemma 2.5 in [9]). Given an $\mathbb{R}^d$-valued process $(X_t)_{t \in [0,T]}$, adapted to the filtration $\mathbb{F}$, consider for any $t \in [0,T]$, a version of $\mathcal{L}^1(X_t)$ as defined in (1.2). Then, the $\mathcal{P}(\mathbb{R}^d)$-valued process $(\mathcal{L}^1(X_t))_{t \in [0,T]}$ is adapted to $\mathbb{F}^0$. If, moreover, $(X_t)_{t \in [0,T]}$ has continuous paths and satisfies $\mathbb{E}[\sup_{0 \le t \le T} |X_t|^p] < \infty$, then we can find a version of each $\mathcal{L}^1(X_t)$, $t \in [0,T]$, such that the process $(\mathcal{L}^1(X_t))_{t \in [0,T]}$ has continuous paths in $\mathcal{P}_p(\mathbb{R}^d)$ and is $\mathbb{F}^0$-adapted.

Consider now the unique strong solution $X = (X_t)_{t \in [0,T]}$ of (1.1). Lemma 2.1 above, and Proposition 2.9, Remark 2.10 in [9] imply that there exists a version of each $\mathcal{L}^1(X_t)$, $t \in [0,T]$, such that the process $(\mathcal{L}^1(X_t))_{t \in [0,T]}$ has continuous paths in $\mathcal{P}_p(\mathbb{R}^d)$ and that it provides a version of the conditional law of $X_t$ given $\mathbb{W}^0$.

**Remark 2.2** (Remark 2.3 in [9]). With a slight abuse of notation, we shall not distinguish a random variable $X$ constructed on $(\Omega^0, \mathcal{F}^0, \mathbb{P}^0)$ (resp. $(\Omega^1, \mathcal{F}^1, \mathbb{P}^1)$) with its natural extension $\tilde{X} : (\omega^0, \omega^1) \mapsto X(\omega^0)$ (resp. $\tilde{X} : (\omega^0, \omega^1) \mapsto X(\omega^1)$) on $(\Omega, \mathcal{F}, \mathbb{P})$. Similarly, for a sub-$\sigma$-algebra $\mathcal{G}^0$ of $\mathcal{F}^0$ (resp. $\mathcal{G}^1$ of $\mathcal{F}^1$), we shall often just write $\mathcal{G}^0$ (resp. $\mathcal{G}^1$) for the sub-$\sigma$-algebra $\mathcal{G}^0 \otimes \{\emptyset, \Omega^1\}$ (resp. $\{\emptyset, \Omega^0\} \otimes \mathcal{G}^1$).

**Lemma 2.3** (General Minkowski inequality). For every $p \in [1, \infty)$, for every process $(X_t)_{t \in [0,T]}$ and for every $T > 0$,

$$ \left\| \int_{0}^{T} X_t dt \right\|_{p} \le \int_{0}^{T} \| X_t \|_{p} dt. $$

**Lemma 2.4** (Burkölder-Davis-Gundy inequality). For every $p \in (0, \infty)$, there exist two positive constants $c_p, C_p$ such that, for every continuous local martingale $X = (X_t)_{t \in [0,T]}$ which vanishes at 0,

$$ c_p \|(\langle X \rangle_T)^{\frac{1}{2}}\|_p \le \left\| \sup_{0 \le t \le T} |X_t| \right\|_p \le C_p \|(\langle X \rangle_T)^{\frac{1}{2}}\|_p. $$

**Lemma 2.5** ('A la Gronwall' Lemma). Let $f : [0, T] \to \mathbb{R}_+$ be a Borel, locally bounded, non-negative and non-decreasing mapping, let $g : [0, T] \to \mathbb{R}_+$ be a non-negative and non-decreasing mapping such that:

$$ f(t) \le C_1 \int_0^t f(s) ds + C_2 \left( \int_0^t f^2(s) s ds \right)^{\frac{1}{2}} + g(t), \quad \forall t \in [0, T], $$

where $C_1$ and $C_2$ are two positive constants. Then for any $t \in [0, T]$,

$$ f(t) \le 2 \exp\{(2C_1 + C_2^2)t\}g(t). $$

We refer to Section 7.8 in [27] for the proofs of Lemma 2.3, Lemma 2.4 and Lemma 2.5. Moreover, we have the following result on the p-Wasserstein distance between conditional laws, whose proof is postponed to Appendix A.

**Lemma 2.6.** For $Y_1, Y_2$ two random variables on $(\Omega, \mathcal{F}, \mathbb{P})$ with a finite p-th moment, $p \in [1, \infty)$, we have

$$ \left\| \mathcal{W}_p(\mathcal{L}^1(Y_1), \mathcal{L}^1(Y_2)) \right\|_p \le \|Y_1 - Y_2\|_p. $$

**Remark 2.7.** From Lemma 2.6, we deduce that for every random variable $X \in L^p(\Omega)$, we have

$$ \| \mathcal{W}_p(\mathcal{L}^1(X), \delta_0) \|_p \le \| X \|_p. $$
### 3 Convergence rate of the Euler scheme

In this section, we prove the convergence rate of the Euler scheme as described in the following proposition.

**Proposition 3.1.** Let $X = (X_t)_{t \in [0,T]}$ denote the unique solution of the McKean-Vlasov equation (1.1), and let $\bar{X}^M = (\bar{X}_t^M)_{t \in [0,T]}$ denote the process defined by the continuous Euler scheme (1.5). Under Assumptions 1, 2 and 3, there exists a constant $C$ depending on $d, p, L, T, \rho$ and $\|X_0\|_p$ such that

$$ \|X - \bar{X}^M\|_{L^p} = \left\| \sup_{0 \le t \le T} |X_t - \bar{X}_t^M| \right\|_p \le C h^{\frac{1}{2} \wedge \rho}. $$

The proof of Proposition 3.1 needs the two following lemmas whose proofs are postponed to Appendix A. In this paper, a constant denoted by $C_{p_1, ..., p_n}$ is a constant depending on parameters $p_1, ..., p_n$, whose value can change line ton line.

**Lemma 3.2.** Let $(\bar{X}_t^M)_{t \in [0,T]}$ be the process defined by the continuous extension of the Euler scheme (1.5). Assume Assumptions 1 and 2 hold. Then, for every $M \ge 1$, $(\bar{X}_t^M)_{t \in [0,T]} \in L^p([0, T] \times \Omega)$ and there exists a non-negative constant $C$ such that

$$ \left\| \sup_{0 \le t \le T} |\bar{X}_t^M| \right\|_p \le C(1 + \|X_0\|_p), $$

where $C$ depends on $p, T, L$ and the coefficients $b, \sigma$ and $\sigma^0$.

The following lemma establishes that $\mathcal{L}^1(\bar{X}_t)$ is a version of the conditional law of $\bar{X}_t$ given the common noise $W^0$, whose proof is postponed to Appendix A.

**Lemma 3.3.** Assume that Assumption 1 holds with some $p \in [2, +\infty)$. Let $\bar{X}^M = (\bar{X}_t^M)_{t \in [0,T]}$ be the process defined by the continuous Euler scheme (1.5). Then for any $t \in [0, T]$, $\mathcal{L}^1(\bar{X}_t)$ provides a version of the conditional distribution of $\bar{X}_t$ given $W^0$. Moreover, we can find a version of $(\mathcal{L}^1(\bar{X}_t))_{t \in [0,T]}$ such that it is $\mathbb{F}^0$-adapted and has continuous paths in $\mathcal{P}_p(\mathbb{R}^d)$.

In the subsequent discussion, to simplify notation, we directly denote by $(\mathcal{L}^1(X_t))_{t \in [0,T]}$ and $(\mathcal{L}^1(\bar{X}_t))_{t \in [0,T]}$ the versions with continuous paths in $\mathcal{P}_p(\mathbb{R}^d)$.

**Lemma 3.4.** Consider $X = (X_t)_{t \in [0,T]}$ the solution to the McKean-Vlasov equation with common noise (1.1). Assume Assumption 1 and 2 hold. For every $0 \le s \le t \le T$, we have

$$ \|X_t - X_s\|_p \le C(t-s)^{\frac{1}{2}}, $$

where the constant $C$ depends on $p, L, T$ and the data $(\|X_0\|_p, b, \sigma, \sigma^0)$.

For the sake of clarity in the proof of Proposition 3.1, we introduce the following notation, for every $m \in \{0, ..., M-1\}$ and for every $t \in [t_m, t_{m+1})$, we define

$$ \underline{t} := t_m. \qquad (3.1) $$

**Proof of Proposition 3.1.** Denote $(\mu_t)_{t \in [0,T]} = (\mathcal{L}^1(X_t))_{t \in [0,T]}$ and $(\bar{\mu}_t^M)_{t \in [0,T]} = (\mathcal{L}^1(\bar{X}_t^M))_{t \in [0,T]}$ that is well-defined by Lemma 3.3. We write $\bar{X}_s$ and $\bar{\mu}_s$ instead of $\bar{X}_s^M$ and $\bar{\mu}_s^M$ when there is no ambiguity. Using Inequality (1.8) and Lemma 3.2, we deduce that $X = (X_t)_{t \in [0,T]}$, $(\bar{X}_t^M)_{t \in [0,T]}$ belongs to $L^p([0, T] \times \Omega)$.
Consequently $(\mu_t)_{t \in [0,T]}$ and $(\bar{\mu}_t^M)_{t \in [0,T]}$ take values in $C([0,T], \mathcal{P}_p(\mathbb{R}^d))$. Fix $t \in [0,T]$, by Minkowski's inequality, we get that

$$ \left\| \sup_{0 \le u \le t} |X_u - \bar{X}_u^M| \right\|_p \le \left\| \sup_{0 \le u \le t} \int_0^u \left( b(s, X_s, \mu_s) - b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) ds \right\|_p \\ + \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) dW_s \right\|_p \\ + \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma^0(s, X_s, \mu_s) - \sigma^0(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) dW_s^0 \right\|_p. \quad (3.2) $$

In the following proof, we provide an upper bound for each term of the right-hand side of (3.2). For the first term, by general Minkowski's inequality (Lemma 2.3), we have

$$ \begin{aligned} \left\| \sup_{0 \le u \le t} \int_0^u \left( b(s, X_s, \mu_s) - b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) ds \right\|_p &\le \left\| \int_0^t \left| b(s, X_s, \mu_s) - b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right| ds \right\|_p \\ &\le \int_0^t \left\| b(s, X_s, \mu_s) - b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right\|_p ds \\ &\le \int_0^t \left\| b(s, X_s, \mu_s) - b(\underline{s}, X_s, \mu_s) \right\|_p ds + \int_0^t \left\| b(\underline{s}, X_s, \mu_s) - b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right\|_p ds \\ &\le L \int_0^t \left\| (s - \underline{s})^\rho (1 + |X_s| + W_p(\mu_s, \delta_0)) \right\|_p ds + L \int_0^t \left\| |X_s - \bar{X}_s| + W_p(\mu_s, \bar{\mu}_s) \right\|_p ds, \end{aligned} $$

where we use the $\rho$-Hölder and $L$-Lipschitz continuity of $b$. From the estimate (1.8) of the solution process $(X_t)_{t \in [0,T]}$, we can deduce that

$$ \begin{aligned} & \left\| \sup_{0 \le u \le t} \int_0^u \left( b(s, X_s, \mu_s) - b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) ds \right\|_p \\ &\le L \sup_{0 \le u \le t} |u - \underline{u}|^\rho \int_0^t (1 + \|X_s\|_p + \|W_p(\mu_s, \delta_0)\|_p) ds + L \int_0^t (\|X_s - \bar{X}_s\|_p + \|W_p(\mu_s, \bar{\mu}_s)\|_p) ds \\ &\le L h^\rho \int_0^t (1 + 2\|X_s\|_p) ds + 2L \int_0^t \|X_s - \bar{X}_s\|_p ds \\ &\le L T h^\rho \left( 1 + \left\| \sup_{0 \le s \le t} |X_s| \right\|_p \right) + 2L \int_0^t (\|X_s - X_s\|_p + \|X_s - \bar{X}_s\|_p) ds \\ &\le C_{p,T,L,b,\sigma,\sigma^0} (1 + \|X_0\|_p) h^\rho + 2L \int_0^t \left( \|X_s - X_s\|_p + \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p \right) ds. \end{aligned} \quad (3.3) $$

We apply Lemma 3.4 to Inequality (3.3) to derive that

$$ \begin{aligned} & \left\| \sup_{0 \le u \le t} \int_0^u \left( b(s, X_s, \mu_s) - b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) ds \right\|_p \\ &\le C_{p,T,L,b,\sigma,\sigma^0} (1 + \|X_0\|_p) h^\rho + L \int_0^t \left( \|X_s - X_s\|_p + \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p \right) ds \\ &\le C_{p,T,L,b,\sigma,\sigma^0} (1 + \|X_0\|_p) h^\rho + L T \kappa \sup_{0 \le s \le T} |s - \underline{s}|^{\frac{1}{2}} + \int_0^t \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p ds \end{aligned} $$
$$ \leq C_{p,T,L,b,\sigma,\sigma^0,\|X_0\|_p} h^{\frac{1}{2}\wedge\rho} + \int_0^t \left\| \sup_{0\leq u\leq s} |X_u - \bar{X}_u| \right\|_p ds. \quad (3.4) $$

For the second and the last terms of Inequality (3.2), the computations are very similar since $\sigma$ and $\sigma^0$ have the same regularity as that of $b$. By the Burkölder-Davis-Gundy inequality (Lemma 2.4), we have

$$ \begin{aligned} & \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) dW_s \right\|_p \\ & \le C_{d,p} \left\| \left( \int_0^t \left| \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right|^2 ds \right)^{\frac{1}{2}} \right\|_p = C_{d,p} \left\| \int_0^t \left| \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right|^2 ds \right\|_{\underline{p}}^{\frac{1}{2}} \\ & \le C_{d,p} \left( \int_0^t \left\| \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right\|^2 \right)_{\underline{p}}^{\frac{1}{2}} ds \right\|_{\underline{p}}^{\frac{1}{2}} \\ & \le C_{d,p} \left( \int_0^t \left\| \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, X_s, \mu_s) \right\|_p^2 ds + \int_0^t \left\| \sigma(\underline{s}, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right\|_p^2 ds \right)^{\frac{1}{2}}, \end{aligned} $$

where we used the Minkowski's and Young's inequalities. Due to the $\rho$-Hölder and $L$-Lipschitz continuity of the mapping $\sigma$, we obtain that

$$ \begin{aligned} & \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) dW_s \right\|_p \\ & \le C_{d,p} \left( \int_0^t \left\| (s-\underline{s})^{2\rho} (1+|X_s| + W_p(\mu_s, \delta_0)) \right\|_p^2 ds + L^2 \int_0^t \left\| |X_s - \bar{X}_s| + W_p(\mu_s, \bar{\mu}_s) \right\|_p^2 ds \right)^{\frac{1}{2}} \\ & \le C_{d,p} \left[ \sup_{0 \le u \le T} |s-\underline{s}|^{2\rho} \int_0^t \left( 1 + \|X_s\|_p + \|W_p(\mu_s, \delta_0)\|_p \right) ds + 2L^2 \int_0^t \left( \|X_s - \bar{X}_s\|_p^2 + \|W_p(\mu_s, \bar{\mu}_s)\|_p^2 \right) ds \right]^{\frac{1}{2}} \\ & \le C_{d,p} \left( h^{2\rho} \int_0^t \left( 1 + 2 \sup_{0 \le u \le s} |X_u| \right) ds + 4L^2 \int_0^t \|X_s - \bar{X}_s\|_p^2 ds \right)^{\frac{1}{2}}. \end{aligned} $$

Applying the estimate (1.8) yields to

$$ \begin{aligned} & \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) dW_s \right\|_p \\ & \le C_{d,p} \left( C_{p,L,T,b,\sigma,\sigma^0} (1 + \|X_0\|_p)^2 h^{2\rho} + 8L^2 \int_0^t \left( \|X_s - X_\underline{s}\|_p^2 + \|X_\underline{s} - \bar{X}_s\|_p^2 \right) ds \right)^{\frac{1}{2}} \\ & \le C_{d,p,L,T,b,\sigma,\sigma^0} \left[ (1 + \|X_0\|_p) h^\rho + \left( \int_0^t \left( \|X_s - X_\underline{s}\|_p^2 + \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right)^2 ds \right)^{\frac{1}{2}} \right]. \end{aligned} \quad (3.5) $$

Applying again Lemma 3.4 to Inequality (3.5) yields to

$$ \begin{aligned} & \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma(s, X_s, \mu_s) - \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) dW_s \right\|_p \\ & \le C_{d,p,L,T,b,\sigma,\sigma^0} \left[ (1 + \|X_0\|_p) h^\rho + \left( \int_0^t \left( \|X_s - X_\underline{s}\|_p^2 + \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right)^2 ds \right)^{\frac{1}{2}} \right] \end{aligned} $$
$$ 
\begin{aligned}
&\le C_{d,p,L,T,b,\sigma,\sigma^0} \left[ (1 + \|X_0\|_p) h^\rho + \sup_{0 \le s \le T} |s - \bar{s}|^{\frac{1}{2}} + \left( \int_0^t \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p^2 ds \right)^{\frac{1}{2}} \right] \\
&\le C_{d,p,L,T,b,\sigma,\sigma^0, \|X_0\|_p} h^{\frac{1}{2}\wedge\rho} + \left( \int_0^t \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p^2 ds \right)^{\frac{1}{2}}.
\end{aligned} 
\quad (3.6) $$

We can repeat the same reasoning for $\sigma^0$ and $W^0$ in order to deduce a similar upper bound

$$ 
\begin{aligned}
&\left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma^0(s, X_s, \mu_s) - \sigma^0(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right) dW_s^0 \right\|_p \\
&\le C_{d,p,L,T,b,\sigma,\sigma^0, \|X_0\|_p} h^{\frac{1}{2}\wedge\rho} + \left( \int_0^t \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p^2 ds \right)^{\frac{1}{2}}.
\end{aligned} 
\quad (3.7) $$

We plug the Inequalities (3.4), (3.6) and (3.7) into Inequality (3.2) to get that

$$ 
\begin{aligned}
\left\| \sup_{0 \le s \le t} |X_s - \bar{X}_s| \right\|_p &\le \int_0^t \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p ds + \left( \int_0^t \left\| \sup_{0 \le u \le s} |X_u - \bar{X}_u| \right\|_p^2 ds \right)^{\frac{1}{2}} \\
&\quad + C_{d,p,L,T,b,\sigma,\sigma^0, \|X_0\|_p} h^{\frac{1}{2}\wedge\rho}.
\end{aligned} 
\quad (3.8) $$

Since $X = (X_t)_{t \in [0,T]}$ and $(\bar{X}_t^M)_{t \in [0,T]}$ belongs to $\mathbb{L}^p$, the application

$$ t \mapsto \left\| \sup_{0 \le s \le t} |X_s - \bar{X}_s| \right\|_p $$

is continuous, non-decreasing and non-negative on $[0, T]$. We conclude this proof by applying Lemma 2.5 to Inequality (3.8) and deduce the existence of a constant $C$ depending on the parameters $d, p, L, T$, and the data $(\|X_0\|_p, b, \sigma, \sigma^0)$ such that we have

$$ \left\| \sup_{0 \le s \le t} |X_s - \bar{X}_s| \right\|_p \le C h^{\frac{1}{2} \wedge \rho}. \quad \square $$

# 4 Convergence rate of the particle system

This section is devoted to proving the convergence rate of the particle method, as described in Theorem 1.2 and Theorem 1.3. To do this, we need the following $N$-particle system $(\bar{Y}^1, \dots, \bar{Y}^N)$ without interaction. Recall the definition of $\underline{s}$ in (3.1). Let $W^0, W^1, \dots, W^N$ be the same Wiener processes as defined in (1.7).

$$ \forall 1 \le i \le N, \quad \bar{Y}_t^i = X_0^i + \int_0^t b(\underline{s}, \bar{Y}_s^i, \mathcal{L}^1(\bar{Y}_s^i)) ds + \int_0^t \sigma(\underline{s}, \bar{Y}_s^i, \mathcal{L}^1(\bar{Y}_s^i)) dW_s^i + \int_0^t \sigma^0(\underline{s}, \bar{Y}_s^i, \mathcal{L}^1(\bar{Y}_s^i)) dW_s^0. \quad (4.1) $$

We have the following property of the particles in the system (4.1), whose proof is postponed to Appendix A.

**Lemma 4.1.** The particles $\bar{Y}^1, \dots, \bar{Y}^N$ are identically distributed having the same distribution as $\bar{X}$ defined by the continuous Euler scheme (1.5), and independent conditionally to $W^0$.

Hence, we may still use the same notation $\bar{\mu}_t$ for $\mathcal{L}^1(\bar{Y}_t^i)$ in the proof. In order to prove Theorem 1.2, we will need the following results in addition of Lemma 4.1 (see [20, Corollary 2.14] and [15, Theorem 1] for the proof).
**Lemma 4.2** (Corollary 2.14 in [20]). Suppose $(X_i)_{1\le i\le N}$ are i.i.d. $\mathbb{R}^d$-valued random variables with law $\mu$ and let $\mu_N = \frac{1}{N}\sum_{i=1}^N \delta_{X_i}$ denote the empirical measure. If $\mu \in \mathcal{P}(\mathbb{R}^d)$, $p \ge 1$, then $\mathcal{W}_p(\mu_N, \mu) \to 0$ almost surely, and also $\mathbb{E}[\mathcal{W}_p^p(\mu_N, \mu)] \to 0$ when $N \to +\infty$.

**Lemma 4.3** (Theorem 1 in [15]). Let $\mu \in \mathcal{P}(\mathbb{R}^d)$. Assume that $M_q(\mu) := \int_{\mathbb{R}^d} |x|^q \mu(dx) < \infty$ for some $q > p$. We consider an i.i.d. sequence $(X_k)_{k \ge 1}$ of $\mu$-distributed random variables and, for $N \ge 1$, the empirical measure $\mu_N := \frac{1}{N} \sum_{k=1}^N \delta_{X_k}$. There exists a positive constant $C_{d,p,q}$ such that, for every $N \ge 1$, we have

$$ \mathbb{E}[\mathcal{W}_p^p(\mu_N, \mu)] \le C_{d,p,q} M_q^{p/q}(\mu) \begin{cases} N^{-\frac{1}{2}} + N^{-(q-p)/q} & \text{if } p > \frac{d}{2} \text{ and } q \neq 2p, \\ N^{-\frac{1}{2}} \log(1+N) + N^{-(q-p)/q} & \text{if } p = \frac{d}{2} \text{ and } q \neq 2p, \\ N^{-\frac{p}{d}} + N^{-(q-p)/q} & \text{if } p \in (0, \frac{d}{2}) \text{ and } q \neq \frac{d}{d-p}. \end{cases} $$

**Remark 4.4.** We deduce from Lemma 4.3 that

$$ \| \mathcal{W}_p(\mu_N, \mu) \|_p \le C_{d,p,q} M_q^{1/q}(\mu) \begin{cases} N^{-\frac{1}{2p}} + N^{-(q-p)/qp} & \text{if } p > \frac{d}{2} \text{ and } q \neq 2p, \\ N^{-\frac{1}{2p}} \log(1+N)^{\frac{1}{p}} + N^{-(q-p)/qp} & \text{if } p = \frac{d}{2} \text{ and } q \neq 2p, \\ N^{-\frac{1}{d}} + N^{-(q-p)/qp} & \text{if } p \in (0, \frac{d}{2}) \text{ and } q \neq \frac{d}{d-p}. \end{cases} $$

*Proof of Theorem 1.2.* (i) Recall that, for every $t \in [0, T]$, $\bar{\mu}_t = \mathcal{L}^1(\bar{X}_t)$, $\bar{\mu}_t^N$ is the empirical distribution of the particles $(\bar{X}_t^{i,N})_{1\le i\le N}$ defined in (1.7). Using the Minkowski inequality, we can bound $\sup_{m \in \{1, \dots, M\}} \|\mathcal{W}_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m})\|_p$ by two terms, which we later demonstrate that they converge to 0 at the desired convergence rate. Consider $(\bar{Y}^1, \dots, \bar{Y}^N)$ identically distributed copies of the continuous Euler scheme $\bar{X}$ defined in (4.1). For every $t \in [0, T]$, we define $\nu_t^N = \frac{1}{N} \sum_{i=1}^N \delta_{\bar{Y}_t^i}$. It follows that, for every $m \in \{1, \dots, M\}$,

$$ \begin{aligned} \sup_{k \in \{1, \dots, m\}} \| \mathcal{W}_p(\bar{\mu}_{t_k}^N, \bar{\mu}_{t_k}) \|_p & \le \sup_{k \in \{1, \dots, m\}} \| \mathcal{W}_p(\bar{\mu}_{t_k}^N, \nu_{t_k}^N) \|_p + \sup_{k \in \{1, \dots, m\}} \| \mathcal{W}_p(\bar{\mu}_{t_k}, \nu_{t_k}^N) \|_p \\ & \le \frac{1}{N} \sum_{i=1}^N \sup_{k \in \{1, \dots, m\}} \| \bar{X}_{t_k}^{i,N} - \bar{Y}_{t_k}^i \|_p + \sup_{0 \le t \le t_m} \| \mathcal{W}_p(\bar{\mu}_t, \nu_t^N) \|_p \\ & \le \left\| \sup_{k \in \{1, \dots, m\}} \left| \bar{X}_{t_k}^{1,N} - \bar{Y}_{t_k}^1 \right| \right\|_p + \sup_{0 \le t \le t_m} \| \mathcal{W}_p(\bar{\mu}_t, \nu_t^N) \|_p \\ & \le \left\| \sup_{0 \le t \le t_m} \left| \bar{X}_t^{1,N} - \bar{Y}_t^1 \right| \right\|_p + \sup_{0 \le t \le t_m} \| \mathcal{W}_p(\bar{\mu}_t, \nu_t^N) \|_p, \end{aligned} \quad (4.2) $$

where the first inequality follows from the triangle inequality of the Wasserstein distance and the Minkowski inequality, the second inequality comes from the fact that the measure $\frac{1}{N} \sum_{i=1}^N \delta_{(x_i, y_i)}$ is a coupling of the measure $\frac{1}{N} \sum_{i=1}^N \delta_{x_i}$ and $\frac{1}{N} \sum_{i=1}^N \delta_{y_i}$, and the third inequality is due to the fact that $(\bar{X}^{i,N}, \bar{Y}^i)_{i \in \{1, \dots, N\}}$ are identically distributed random variables. Remark that $(\bar{Y}_t^i)_{t \in [0, T]}$ are conditionally independent copies of the continuous extension $(\bar{X}_t)_{t \in [0, T]}$ of the Euler scheme defined by (1.5). For the first term of (4.2), we apply the Lipschitz continuity of the coefficients $b$, $\sigma$, $\sigma^0$ and the Burkölder-Davis-Gundy inequality (Lemma 2.4) to get

$$ \left\| \sup_{0 \le s \le t} \left| \bar{X}_s^{1,N} - \bar{Y}_s^1 \right| \right\|_p \le \left\| \sup_{0 \le u \le t} \int_0^u \left( b(\underline{s}, \bar{Y}_s^1, \bar{\mu}_s) - b(\underline{s}, \bar{X}_s^{1,N}, \bar{\mu}_s^N) \right) ds \right\|_p \\ + \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma(\underline{s}, \bar{Y}_s^1, \bar{\mu}_s) - \sigma(\underline{s}, \bar{X}_s^{1,N}, \bar{\mu}_s^N) \right) dW_s \right\|_p $$
$$ + \left\| \sup_{0 \le u \le t} \int_0^u \left( \sigma^0(\underline{s}, \bar{Y}_s^1, \bar{\mu}_s) - \sigma^0(\underline{s}, \bar{X}_s^{1,N}, \bar{\mu}_s^N) \right) dW_s^0 \right\|_p \\ \le L \int_0^t \left( \|\bar{X}_u^{1,N} - \bar{Y}_u^1\|_p + \|W_p(\bar{\mu}_u^N, \bar{\mu}_u)\|_p \right) du + \left[ C_{d,p} L \int_0^t \left( \|\bar{X}_u^{1,N} - \bar{Y}_u^1\|_p^2 + \|W_p(\bar{\mu}_u^N, \bar{\mu}_u)\|_p^2 \right) du \right]^{\frac{1}{2}}, $$

where the positive constant $C_{d,p}$ comes from the Burkölder-Davis-Gundy inequality, see Lemma 2.4. It follows from Inequality (4.2)

$$ \left\| \sup_{0 \le s \le t} |\bar{X}_s^{1,N} - \bar{Y}_s^1| \right\|_p \le 2L \int_0^t \left\| \sup_{0 \le r \le u} |\bar{X}_r^{1,N} - \bar{Y}_r^1| \right\|_p du + C_{d,p,L} \left( \int_0^t \left\| \sup_{0 \le r \le u} |\bar{X}_r^{1,N} - \bar{Y}_r^1| \right\|_p^2 du \right)^{\frac{1}{2}} + g(t), \quad (4.3) $$

where

$$ t \in [0, T] \mapsto g(t) = L \int_0^t \sup_{0 \le r \le u} \|W_p(\bar{\mu}_r, \nu_r^N)\|_p du + C_{d,p,L} \left( \int_0^t \sup_{0 \le r \le u} \|W_p(\bar{\mu}_r, \nu_r^N)\|_p^2 du \right)^{\frac{1}{2}}. $$

We also define on $[0, T]$ the mapping

$$ t \mapsto f(t) = \left\| \sup_{0 \le r \le t} |\bar{X}_r^{1,N} - \bar{Y}_r^1| \right\|_p. $$

It follows from Lemma 3.2 that $(\bar{Y}_t^1)_{t \in [0,T]} \in L^p([0,T] \times \Omega)$. By a similar reasoning, we can deduce that $(\bar{X}_t^{1,N})_{t \in [0,T]} \in L^p([0,T] \times \Omega)$. Moreover $(\bar{\mu}_t)_{t \in [0,T]}$ and $(\nu_t^N)_{t \in [0,T]}$ are random variables on $\Omega^0$ taking values in $C([0,T], P_p(\mathbb{R}^d))$. Hence the mappings $f$ and $g$ are continuous, non-decreasing and non-negative on $[0,T]$. We deduce from Inequality (4.3) and Lemma 2.5 that

$$ \left\| \sup_{0 \le s \le t} |\bar{X}_s^{1,N} - \bar{Y}_s^1| \right\|_p \le 2e^{(2L+C_{d,p,L}^2)t} g(t). \quad (4.4) $$

Merging Inequalities (4.2) and (4.4), we apply Lemma 2.5 once more to derive the following essential inequality

$$ \sup_{m \in \{1, \dots, M\}} \|W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m})\|_p \le C_{d,p,L,T} \sup_{0 \le t \le T} \|W_p(\bar{\mu}_t, \nu_t^N)\|_p, \quad (4.5) $$

where $C_{d,p,L,T}$ is a positive constant. Due to Lemma 4.1, the $\bar{Y}^i$'s are identically distributed and independent conditionally to $W^0$ and in $L^p$, by Lemma 4.2, we deduce that for every $t \in [0, T]$, for every $\omega_0 \in \Omega^0$,

$$ \mathbb{E}^1 \left[ W_p^p (\bar{\mu}_t(\omega_0), \nu_t^N(\omega_0)) \right] \xrightarrow[N \to +\infty]{} 0. $$

Moreover we have the following control by the Young's inequality

$$ \begin{aligned} \sup_{0 \le t \le T} \mathbb{E}^1 [W_p^p (\bar{\mu}_t, \nu_t^N)] &\le 2^{p-1} \sup_{0 \le t \le T} \mathbb{E}^1 [W_p^p (\bar{\mu}_t, \delta_0) + W_p^p (\delta_0, \nu_t^N)] \\ &\le 2^{p-1} \frac{1}{N} \sum_{i=1}^N \sup_{0 \le t \le T} \mathbb{E}^1 [|\bar{Y}_t^1|^p] + 2^{p-1} \sup_{0 \le t \le T} \mathbb{E}^1 [|\bar{Y}_t^1|^p] \\ &= 2^p \sup_{0 \le t \le T} \mathbb{E}^1 [|\bar{Y}_t^1|^p]. \end{aligned} \quad (4.6) $$

We conclude by using the $L^p$ version of the dominated convergence theorem and the Inequality (4.5) to deduce that

$$ \sup_{m \in \{1, \dots, M\}} \|W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m})\|_p \le C_{d,p,L,T} \sup_{0 \le t \le T} \|W_p(\bar{\mu}_t, \nu_t^N)\|_p \xrightarrow[N \to +\infty]{} 0. $$
(ii) Additionally, we assume that $\|X_0\|_{p+\varepsilon} < \infty$ for $\varepsilon > 0$. Hence, it follows from Lemma 3.2 that $\bar{Y}^i \in \mathbb{L}^{p+\varepsilon}$ for every $i \in \{1, \dots, N\}$. Moreover the $\bar{Y}^i$'s are i.i.d. given the path of $W^0$ according to Lemma 4.1. Applying Lemma 4.3, we establish the following rate of convergence for every $\omega_0 \in \Omega^0$,

$$ \sup_{0 \le t \le T} \mathbb{E}^1 \left[ W_p^p \left( \bar{\mu}_t(\omega_0), \nu_t^N(\omega_0) \right) \right] \le \tilde{C} \left\{ \begin{array}{ll} N^{-\frac{1}{2}} + N^{-\frac{\varepsilon}{p+\varepsilon}} & \text{if } p > \frac{d}{2} \text{ and } \varepsilon \neq p, \\ N^{-\frac{1}{2}} \left( \log(1+N) \right) + N^{-\frac{\varepsilon}{p+\varepsilon}} & \text{if } p = \frac{d}{2} \text{ and } \varepsilon \neq p, \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p+\varepsilon}} & \text{if } p \in (0, \frac{d}{2}) \text{ and } p + \varepsilon \neq \frac{d}{d-p}, \end{array} \right. $$

where $\tilde{C}$ is a constant independent of $N$. By Inequality (4.6), we can apply the $L^p$ version of the dominated convergence theorem to deduce a similar convergence rate for $\sup_{0 \le t \le T} \|W_p^p(\bar{\mu}_t, \nu_t^N)\|_p$. From the Inequality (4.5) and Remark 4.4, we can conclude that

$$ \begin{aligned} \max_{m \in \{1, \dots, M\}} \|W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m})\|_p &\le \sup_{t \in [0, T]} \|W_p(\bar{\mu}_t^N, \bar{\mu}_t)\|_p \\ &\le C \left\{ \begin{array}{ll} N^{-\frac{1}{2p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > \frac{d}{2} \text{ and } \varepsilon \neq p, \\ N^{-\frac{1}{2p}} \left(\log(1+N)\right)^{\frac{1}{p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = \frac{d}{2} \text{ and } \varepsilon \neq p, \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, \frac{d}{2}) \text{ and } p + \varepsilon \neq \frac{d}{d-p}, \end{array} \right. \end{aligned} $$

where $C$ is a positive constant that depends on the parameters $d, p, L, T$ and the data ($\|X_0\|_{p+\varepsilon}, b, \sigma, \sigma^0$). $\square$

*Proof of Theorem 1.3.* Recall that the processes $(\bar{Y}_t^i)_{t \in [0, T]}$, $1 \le i \le N$ are defined by the continuous Euler scheme (4.1). Then we use the Minkowski inequality to get that

$$ \left\| \sup_{t \in [0, T]} |\bar{X}_t^{1,N} - X_t| \right\|_p \le \left\| \sup_{t \in [0, T]} |\bar{X}_t^{1,N} - \bar{Y}_t^1| \right\|_p + \left\| \sup_{t \in [0, T]} |\bar{Y}_t^1 - X_t| \right\|_p. \quad (4.7) $$

By the inequality (4.4), we have

$$ \left\| \sup_{t \in [0, T]} |\bar{X}_t^{1,N} - \bar{Y}_t^1| \right\|_p \le C_{d,p,L,T} \sup_{t \in [0, T]} \|W_p(\bar{\mu}_t^N, \bar{\mu}_t)\|_p. $$

From Theorem 1.2, we derive the following rate of convergence

$$ \sup_{t \in [0, T]} \|W_p(\bar{\mu}_t^N, \bar{\mu}_t)\|_p \le C_{d,p,L,T} \mathcal{E}_N, \quad (4.8) $$

where

$$ \mathcal{E}_N = \left\{ \begin{array}{ll} N^{-\frac{1}{2p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > \frac{d}{2} \text{ and } \varepsilon \neq p, \\ N^{-\frac{1}{2p}} \left( \log(1+N) \right)^{\frac{1}{p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = \frac{d}{2} \text{ and } \varepsilon \neq p, \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, \frac{d}{2}) \text{ and } p + \varepsilon \neq \frac{d}{d-p}. \end{array} \right. $$

Moreover, it follows from Proposition 3.1 that

$$ \left\| \sup_{t \in [0, T]} |X_t - \bar{Y}_t^1| \right\|_p \le C_{d,p,L,T} h^{\frac{1}{2} \wedge \rho}, \quad (4.9) $$

since $\bar{X}_t^{1,N}$ shares the same Brownian motions as $X_t$, so does $\bar{Y}_t^1$. We plug Inequalities (4.8) and (4.9) into Inequality (4.7) to deduce the rate of convergence in Inequality (1.9). $\square$
# 5 Numerical simulations

This section presents two simulation examples, the first one is a conditional Ornstein-Uhlenbeck process with common noise and the second one is an application from [29] representing an Interbank market model. The simulation codes are available via Google Colab (https://bit.ly/49xgGHG, https://bit.ly/3D4zwK3).

## 5.1 Conditional Ornstein-Uhlenbeck process with common noise

In this section, we present the following simulation example in $\mathbb{R}^d$

$$dX_t = -(X_t - \mathbb{E}^1(X_t))dt + \sigma dW_t + \sigma^0 dW_t^0, \quad X_0 = x_0 \in \mathbb{R}^d \quad (5.1)$$

where $\sigma$ and $\sigma^0$ are two positive scalar. This example generalizes, incorporating common noise, the one analyzed in detail in [20, Section 3.1]. It is obvious that Equation (5.1) admits a unique solution $X = (X_t)_{t \in [0,T]}$, as the drift coefficient $b$ is linear in $x$ and in $\mu$, and the diffusion coefficients $\sigma$ and $\sigma^0$ are constants. A straightforward computation yields the following closed-form expression for the unique solution

$$X_t = e^{-t}x_0 + \sigma \int_0^t e^{-(t-s)}dW_s + \sigma^0 W_t^0. \quad (5.2)$$

We set $x_0 = 0$, $T = 1$ and $\sigma = \sigma^0 = \sqrt{0.2}$. In a first time, we try to compute the simulation error by computing the $L^2$-error between the solution process $(X_t)_{t \in [0,T]}$ and its approximation.

The numerical study of this equation is divided into two parts. In the first part, we consider the case where the dimension $d = 2$ and analyze the $L^2$-error $\{\mathbb{E}[\sup_{t \in [0,T]} |X_t - \bar{X}_t^{1,N}|^2]\}^{1/2}$, where $(X_t)_{t \in [0,T]}$ is the true solution defined by (5.2) and $(\bar{X}_t^{1,N})_{t \in [0,T]}$ is the first particle path defined by (1.7) sharing the same Brownian motion of $(X_t)_{t \in [0,T]}$. This $L^2$-error is approximated by

$$\varepsilon_N^2 = \frac{1}{30} \sum_{j=1}^{30} \frac{1}{N} \sum_{i=1}^{N} \max_{1 \le m \le M} |\bar{X}_{t_m}^{i,N,j} - X_{t_m}^{i,j}|^2, \quad (5.3)$$

using the time discretization number $M = 100$. Here, in (5.3), for every fixed $i \in \{1, ..., N\}$, $j \in \{1, ..., 30\}$, the random variables $\bar{X}_{t_m}^{i,N,j}$ is simulated using the particle method (1.6), and $X_{t_m}^{i,j}$ is computed using the explicit solution (5.2) both sharing the same Brownian motion. Note that the mean error over $j \in \{1, ..., 30\}$ is used for the Monte-Carlo estimation of the expectation. In Figure 1, we display the log-log error of $\varepsilon_N$ as a mapping of the number of particles $N \in \{2^6, ..., 2^{16}\}$.

The second part focuses on the one-dimensional case. Remark that this example has an explicit formula for the density of $(\mu_t)_{t \in [0,T]} = (\mathcal{L}^1(X_t))_{t \in [0,T]}$, given by

$$\frac{\mu_t(x)}{dx} = \frac{1}{\sqrt{2\pi\sigma^2(1-e^{-t})}} \exp\left(-\frac{|x - X_0 - \sigma_0 W_t^0|^2}{2\sigma^2(1-e^{-t})}\right).$$

Hence, we can compute an approximated density $(\hat{\mu}_t^{N,h,\eta})_{t \in [0,T]}$ by using the kernel method as presented in [17] for a standard McKean-Vlasov equation. The estimator $\hat{\mu}_T^{N,h,\eta}$ of the density $\mu_T$ at the final time $T = 1$ is defined, for every $x \in \mathbb{R}$, by

$$\hat{\mu}_T^{N,h,\eta}(x) := N^{-1} \sum_{i=1}^{N} \eta^{-1} K \left( \eta^{-1} \left( x - \bar{X}_T^{i,N} \right) \right),$$
where $K(\cdot)$ is the kernel function and $\eta$ is the bandwidth. We will apply a Gaussian-based kernel of order $l = 5$, $K(x) = \frac{1}{8}(15 - 10x^2 + x^4)\phi(x)$ where $\phi(x) = (2\pi)^{-1/2} \exp(-\frac{x^2}{2})$, with the bandwidth $\eta = N^{-1/(2(l+1)+1)}$ chosen according to Corollary 2.11 of [17], since $\mu_T \in C_b^\infty(\mathbb{R})$. The estimation error is computed by

$$ \mathcal{E}_N^2 = \frac{1}{30} \sum_{j=1}^{30} \max_{x \in \mathcal{D}} \left| (\hat{\mu}_T^{N,h,\eta})_j(x) - \mu_T(x) \right|^2, \quad (5.4) $$

where the domain $\mathcal{D}$ is a uniform grid on $[-3, 3]$ chosen according to the trajectory of the simulations $(\bar{X}^{i,N})_{1 \le i \le N}$. Figure 2 shows the log-log error of this density estimation.

<figure><img src="image_3.png" /><figcaption>Figure 1: Log-log error (5.3) between $(X_t)_{t \in [0,T]}$ and Figure 2: Log-log error (5.4) between $\hat{\mu}_T^{N,h,\eta}$ and $\mu_T$ $(\bar{X}_{t_m})_{0 \le m \le M}$ (slope = -0.35) (slope = -0.4)</figcaption></figure>

## 5.2 Interbank market model

We consider an application from [29] where they study risk-sensitive mean field games (MFG) with common noise. It is an infinite population model where the log-reserve $(X_t^{i,N})_{t \in [0,T]}$ of the bank $i \in I$ at time $t \in [0, T]$ satisfies the following dynamics

$$ dX_t^{i,N} = (a(\bar{X}_t - X_t^{i,N}) + u_t^i + b(t))dt + \sigma\sqrt{1-\rho^2}dW_t^i + \sigma\rho dW_t^0, \quad (5.5) $$

with $\bar{X}_t = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^N X_t^{i,N}$ represents the limiting market state which is following the dynamic

$$ d\bar{X}_t = (\bar{u}_t + b(t))dt + \sigma\rho dW_t^0, \quad (5.6) $$

where $\bar{u}_t = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^N u_t^i$. The transaction $u_t^i \in \mathbb{R}$ represents the money that the bank lends to or borrows from the central bank during the market activity at each time $t$, the market shock is simulated by $W_t^0 \in \mathbb{R}$ which is independent of the shock received by the bank $W_t^i \in \mathbb{R}$.

The parameter $a \in \mathbb{R}$ represents the mean reversion rate of the bank's reserve towards the market state. The bank's liquidity before market activity at each time $t$ is denoted by $b(t)$. The volatility of the log-reserve of the bank with respect to its own local shock (underlying uncertainty source) is given by $\sigma\rho \in \mathbb{R}$. Meanwhile, the volatility of the log-reserve with respect to the global shock that affects the market (i.e. the macroeconomic factors), is characterized by $\sigma\sqrt{1-\rho^2} \in \mathbb{R}$. As can be seen from above equations,
an instantaneous coefficient $0 \le \rho \le 1$ is a common multiplier factor for the shock delivered by the bank itself and by the environment.

The error analysis of this Interbank market model will proceed as follows. First, we fix $b(t) = 1$ for every $t \in [0, T]$, set $a = 10$ and choose $\rho \in (0, 1)$ to ensure the presence of both non-zero idiosyncratic noises and common noise. Next, we compute the optimal transaction rate $(u_t^{i,*})_{t \in [0,T]}$ with respect to the cost function defined by [29, Equation (134)]. It is important to note that the explicit formula for $(u_t^{i,*})_{t \in [0,T]}$ depends on parameters $(\Pi_t, \Lambda_t, \Upsilon_t, \Delta_t, \Gamma_t, \Psi_t)_{t \in [0,T]}$, as specified in [29, Equations (139-144)] and these parameters are determined by a system of ODEs, which can be numerically solved using the Python library `solve_ivp` from `scipy.integrate`.

After computing the optimal transaction rate $(u_t^{i,*})_{t \in [0,T]}$, the value of $\bar{u}_t^*$ can be directly deduced from [29, Equation (145)] and plugged into the dynamics of $(\bar{X}_t)_{t \in [0,T]}$ defined by (5.6). Consequently, it can be numerically solved using a standard Euler scheme for diffusion processes (see e.g. [27, Section 7.1]). On the other hand, the dynamics described by (5.5), incorporating the optimal transaction rates $(u_t^{i,*})_{t \in [0,T]}$, can be computed using the particle method presented in this paper.

The simulation error is therefore evaluated as the difference between $\frac{1}{N} \sum_{i=1}^N X_t^{i,N}$, obtained from the following dynamics

$$dX_t^{i,N} = \left\{ 10 \left( \frac{1}{N} \sum_{i=1}^{N} X_t^{i,N} - X_t^{i,N} \right) + u_t^{i,*} + 1 \right\} dt + \sigma \sqrt{1 - \rho^2} dW_t^i + \sigma \rho dW_t^0, \quad (5.7)$$

and $\bar{X}_t$ computed by (5.6), and we perform 30 independent experiments for Monte-Carlo approximation of the expectation. More specifically, we set the time discretization number $M = 100$ for both (5.6) and (5.7), and compute the simulation error using the following formula:

$$\mathcal{E}_N^2 = \frac{1}{30} \sum_{j=1}^{30} \max_{1 \le m \le M} \left\| \frac{1}{N} \sum_{i=1}^{N} X_{t_m}^{i,N,j} - \bar{X}_{t_m}^j \right\|^2. \quad (5.8)$$

where $j$ in the superscript of $\bar{X}_{t_m}^{i,N,j}$ and $\bar{X}_{t_m}^j$ denotes the index of the independent experiment, and for each fixed $j \in \{1, ..., 30\}$, $\bar{X}_{t_m}^{i,N,j}$, $1 \le i \le N$ is computed using the particle method (1.6) based on the dynamics in (5.7), while $\bar{X}_{t_m}^j$ is computed using a standard Euler scheme for the diffusion process (5.6). Figure 3 illustrates the log-log error of $\mathcal{E}_N^2$ defined by (5.8) and Figure 4 displays 10 simulated paths of $(X_t^{i,N})_{t \in [0,T]}$ under the above setting.

<figure><img src="image_4.png" /><figcaption>Figure 3: Log-log error of $\mathcal{E}_N$ defined by (5.8) (slope = -0.4)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: 10 paths of particles $(X_t^{i,N})_{t \in [0,T]}$</figcaption></figure>
# A Appendix

*Proof of Lemma 2.6.* For a fixed $\omega^0 \in \Omega^0$, we have

$$ \begin{aligned} \mathcal{W}_p^p(\mathcal{L}^1(Y^1)(\omega^0), \mathcal{L}^1(Y^2)(\omega^0)) &= \mathcal{W}_p^p(\mathcal{L}(Y^1(\omega^0, \cdot)), \mathcal{L}(Y(\omega^0, \cdot))) \\ &\le \mathbb{E}^1 [|Y_1(\omega^0, \cdot) - Y_2(\omega^0, \cdot)|^p] \end{aligned} $$

This inequality is true for every $\omega^0 \in \Omega^0$, hence we have

$$ \mathbb{E}[\mathcal{W}_p^p(\mathcal{L}^1(Y_1), \mathcal{L}^1(Y_2))] \le \mathbb{E}[|Y_1 - Y_2|^p] $$

Then we have

$$ \|\mathcal{W}_p(\mathcal{L}^1(Y_1), \mathcal{L}^1(Y_2))\|_p \le \|Y_1 - Y_2\|_p. \quad \square $$

We present the proofs of Lemma 3.2 and Lemma 3.4 required by the proof of Proposition 3.1.

*Proof of Lemma 3.2.* We reason by forward induction. By Minkowski's inequality, we have

$$ \begin{aligned} \|\bar{X}_{t_1}\|_p &= \|\bar{X}_0 + h \cdot b(0, \bar{X}_0, \bar{\mu}_0)t_1 + \sigma(0, \bar{X}_0, \bar{\mu}_0)W_{t_1} + \sigma^0(0, \bar{X}_0, \bar{\mu}_0)W_{t_1}^0\|_p \\ &\le \|\bar{X}_0\|_p + h\|b(0, \bar{X}_0, \bar{\mu}_0)\|_p + \sqrt{h}\|\sigma(0, \bar{X}_0, \bar{\mu}_0)Z_1\|_p + \sqrt{h}\|\sigma^0(0, \bar{X}_0, \bar{\mu}_0)Z_1^0\|_p \end{aligned} $$

Since the solution $(X_t)_{t \in [0,T]}$ and the continuous scheme $(\bar{X}_t)_{t \in [0,T]}$ have the same initial condition $X_0$, we have, $\mathbb{P}$-almost surely

$$ \bar{X}_0 = X_0 \in L^p(\Omega) \text{ and } \bar{\mu}_0 = \mathcal{L}(X_0). $$

For every $0 \le m \le M - 1$, $Z_{m+1}$ and $Z_{m+1}^0$ are independent of $(\bar{X}_{t_m}, \dots, \bar{X}_{t_1}, X_0)$ and identically distributed. Moreover, since $b$, $\sigma$ and $\sigma^0$ have a linear growth, we have by Minkowski's inequality

$$ \begin{aligned} \|\bar{X}_{t_1}\|_p &\le \|X_0\|_p + C_{p,L,T,b,\sigma,\sigma^0}(1 + \|X_0\|_p + \|\mathcal{W}_p(\mu_0, \delta_0)\|_p)(h + \sqrt{h}\|Z_1\|_p + \sqrt{h}\|Z_1^0\|_p) \\ &\le \|X_0\|_p + C_{p,L,T,b,\sigma,\sigma^0}(1 + 2\|X_0\|_p)(h + \sqrt{h}). \end{aligned} $$

For $1 \le m \le M - 1$, we repeat the same reasoning to show that

$$ \begin{aligned} \|\bar{X}_{t_{m+1}}\|_p &\le \|\bar{X}_{t_m}\|_p + C_{p,L,T,b,\sigma,\sigma^0}(1 + 2\|\bar{X}_{t_m}\|_p)(h + C_p\sqrt{h}) \\ &\le C_{p,L,T,b,\sigma,\sigma^0,h}^0 \|\bar{X}_{t_m}\|_p + C_{p,L,T,b,\sigma,\sigma^0,h}^1. \end{aligned} $$

By forward induction, we have for every $0 \le m \le M$,

$$ \|\bar{X}_{t_m}\|_p \le C(1 + \|X_0\|_p), \quad (A.1) $$

where $C$ is a positive constant depending on the parameters $p$, $L$, $T$, $h$, $M$ and the coefficients $b$, $\sigma$ and $\sigma^0$. For every $t \in [0, T]$, we recall the notation $\underline{t} = t_m$ where $m \in \{0, \dots, M\}$ such that $t_m \le t < t_{m+1}$. By definition (1.5) of the continuous Euler scheme $(\bar{X}_t)_{t \in [0,T]}$, it is $\mathbb{F}$-adapted and right continuous so is progressively measurable. Moreover, we have

$$ \begin{aligned} \left\| \sup_{0 \le t \le T} |\bar{X}_t| \right\|_p &= \left\| \sup_{0 \le t \le T} \left| X_0 + \int_0^t b(\underline{s}, \bar{X}_s, \bar{\mu}_s) ds + \int_0^t \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) dW_s + \int_0^t \sigma^0(\underline{s}, \bar{X}_s, \bar{\mu}_s) dW_s^0 \right| \right\|_p \\ &\le \|X_0\|_p + \left\| \sup_{0 \le t \le T} \left| \int_0^t b(\underline{s}, \bar{X}_s, \bar{\mu}_s) ds \right| \right\|_p + \left\| \sup_{0 \le t \le T} \left| \int_0^t \sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s) dW_s \right| \right\|_p \end{aligned} $$
$$
\begin{align}
&+ \left\| \sup_{0 \le t \le T} \left| \int_0^t \sigma^0(\underline{s}, \bar{X}_s, \bar{\mu}_s) dW_s^0 \right| \right\|_p \nonumber \\
&\le \|X_0\|_p + \left\| \int_0^T \left| b(\underline{s}, \bar{X}_s, \bar{\mu}_s) \right| ds \right\|_p + C_{d,p} \left\| \left( \int_0^T |\sigma(\underline{s}, \bar{X}_s, \bar{\mu}_s)|^2 ds \right)^{\frac{1}{2}} \right\|_p \tag{A.2} \\
&\quad + C_{d,p} \left\| \left( \int_0^T |\sigma^0(\underline{s}, \bar{X}_s, \bar{\mu}_s)|^2 ds \right)^{\frac{1}{2}} \right\|_p \nonumber \\
&\le \|X_0\|_p + C_{L,T,b,\sigma,\sigma^0} \left[ \int_0^T \left\| (1 + |\bar{X}_s| + W_p(\bar{\mu}_s, \delta_0)) \right\|_p ds + \left\| \left( \int_0^T |1 + |\bar{X}_s| + W_p(\bar{\mu}_s, \delta_0)|^2 ds \right)^{\frac{1}{2}} \right\|_p \right] \tag{A.3}
\end{align}
$$

where we use Minkowski's inequality for the first inequality, the second inequality follows from Lemma 2.4 and the linear growth of the coefficients. By applying Young's inequality on (A.3), we can infer that

$$
\begin{align*}
& \left\| \sup_{0 \le t \le T} |\bar{X}_t| \right\|_p \\
&\le \|X_0\|_p + C_{L,T,d,p,b,\sigma,\sigma^0} \left[ \int_0^T (1+2\|\bar{X}_s\|_p) ds + \sqrt{T} + \left( \int_0^T \|\bar{X}_s\|_p^2 ds \right)^{\frac{1}{2}} + \left( \int_0^T W_p^2(\bar{\mu}_s, \delta_0) ds \right)^{\frac{1}{2}} \right] \\
&\le \|X_0\|_p + C_{L,T,d,p,b,\sigma,\sigma^0} \left[ \int_0^T (1+2\|\bar{X}_s\|_p) ds + \left( \int_0^T \|\bar{X}_s\|_p^2 ds \right)^{\frac{1}{2}} \right],
\end{align*}
$$

from Lemma 2.6. Then we deduce the following inequality to which we apply Lemma 2.5

$$
\left\| \sup_{0 \le t \le T} |\bar{X}_t| \right\|_p \le \|X_0\|_p + C_{L,T,d,p,b,\sigma,\sigma^0} \left[ \int_0^T \left\| \sup_{0 \le u \le s} |\bar{X}_u^M| \right\|_p ds + \left( \int_0^T \left\| \sup_{0 \le u \le s} |\bar{X}_u^M| \right\|_p^2 ds \right)^{\frac{1}{2}} \right]. \quad (A.4)
$$

Since $\bar{X}^M \in \mathbb{L}^p$, the application $t \mapsto \left\| \sup_{0 \le u \le t} \bar{X}_u^M \right\|_p$ is continuous, non-negative and non-decreasing on $[0, T]$, the final estimate follows from Lemma 2.5

$$
\left\| \sup_{0 \le t \le T} |\bar{X}_t| \right\|_p \le \tilde{C} (1 + \|X_0\|_p),
$$

where $\tilde{C}$ is a positive constant depending on $p$, $L$, $T$, $h$, $M$ and the coefficients $b$, $\sigma$ and $\sigma^0$. Hence
$\bar{X} \in \mathbb{L}^p([0, T] \times \Omega)$. $\square$

Proof of Lemma 3.3. The first step of the proof is inspired by the proof of [9, Proposition 2.9]. The first additional step is to combine the results of the Lemmas 2.1 and 3.2, that is $(\bar{X}_t)_{t \in [0,T]} \in \mathbb{L}^p([0,T] \times \Omega)$ and $\mathcal{L}^1(\bar{X}_t)_{t \in [0,T]}$ has continuous paths in $\mathcal{P}_p(\mathbb{R}^d)$ and is $\mathbb{F}^0$-adapted. The second one follows directly from the construction of the Euler scheme $(\bar{X}_t)_{t \in [0,T]}$. There exists a unique process that satisfies Equation (1.5) for the Brownian motions $(W_t)_{t \in [0,T]}$ and $(W_t^0)_{t \in [0,T]}$. The second step is a consequence of Lemma 2.1 and Lemma 3.2. $\square$

*Proof of Lemma 3.4.* Set $(\mu_t)_{t \in [0,T]} = (\mathcal{L}^1(X_t))_{t \in [0,T]}$. Fix $s, t \in [0,T]$ such that $s \le t$. Since $(X_t)_{t \in [0,T]}$ solves the McKean-Vlasov equation with common noise (1.1), we have

$$
\|X_t - X_s\|_p = \left\| \int_s^t b(u, X_u, \mu_u) du + \int_s^t \sigma(u, X_u, \mu_u) dW_u + \int_s^t \sigma^0(u, X_u, \mu_u) dW_u^0 \right\|_p
$$
$$
\begin{align*}
&\le \left\| \int_s^t b(u, X_u, \mu_u) du \right\|_p + \left\| \int_s^t \sigma(u, X_u, \mu_u) dW_u \right\|_p + \left\| \int_s^t \sigma^0(u, X_u, \mu_u) dW_u^0 \right\|_p \\
&\le \int_s^t \|b(u, X_u, \mu_u)\|_p du + C_{d,p} \left\| \left( \int_s^t \sigma(u, X_u, \mu_u)^2 du \right)^{\frac{1}{2}} \right\|_p + C_{d,p} \left\| \left( \int_s^t \sigma^0(u, X_u, \mu_u)^2 du \right)^{\frac{1}{2}} \right\|_p \\
&= \int_s^t \|b(u, X_u, \mu_u)\|_p du + C_{d,p} \left\| \int_s^t \sigma(u, X_u, \mu_u)^2 du \right\|_{\frac{p}{2}}^{\frac{1}{2}} + C_{d,p} \left\| \int_s^t \sigma^0(u, X_u, \mu_u)^2 du \right\|_{\frac{p}{2}}^{\frac{1}{2}} \\
&\le \int_s^t \|b(u, X_u, \mu_u)\|_p du + C_{d,p} \left( \int_s^t \|\sigma(u, X_u, \mu_u)^2\|_{\frac{p}{2}} du \right)^{\frac{1}{2}} + C_{d,p} \left( \int_s^t \|\sigma^0(u, X_u, \mu_u)^2\|_{\frac{p}{2}} du \right)^{\frac{1}{2}} \\
&= \int_s^t \|b(u, X_u, \mu_u)\|_p du + C_{d,p} \left( \int_s^t \|\sigma(u, X_u, \mu_u)\|_p^2 du \right)^{\frac{1}{2}} + C_{d,p} \left( \int_s^t \|\sigma^0(u, X_u, \mu_u)\|_p^2 du \right)^{\frac{1}{2}}.
\end{align*}
$$

where we used the Minkowski's inequality at the first inequality, the second one follows from the gen-
eral Minkowski's (Lemma 2.3) and Burkölder-Davis-Gundy inequalities (Lemma 2.4) and the general
Minkowski's inequality provides the last one. We use the linear growth of the coefficients and Minkowski's
and Young's inequalities to get that

$$
\begin{align*}
\|X_t - X_s\|_p &\le C_{L,T,b,\sigma,\sigma^0} \left( \int_s^t \|1 + |X_u| + W_p(\mu_u, \delta_0)\|_p du + 2 \left[ \int_s^t \|1 + |X_u| + W_p(\mu_u, \delta_0)\|_p^2 du \right]^{\frac{1}{2}} \right) \\
&\le C_{L,T,b,\sigma,\sigma^0} \left( \int_s^t (1 + \|X_u\|_p + \|W_p(\mu_u, \delta_0)\|_p) du + 2 \left[ \int_s^t (1 + \|X_u\|_p + \|W_p(\mu_u, \delta_0)\|_p)^2 du \right]^{\frac{1}{2}} \right) \\
&\le C_{L,T,b,\sigma,\sigma^0} \left( \int_s^t (1 + \|X_u\|_p) du + \left[ \int_s^t (1 + \|X_u\|_p)^2 du \right]^{\frac{1}{2}} \right) \\
&\le C_{L,T,b,\sigma,\sigma^0} \left( \int_s^t (1 + \|X_u\|_p) du + \left[ \int_s^t (1 + \|X_u\|_p)^2 du \right]^{\frac{1}{2}} \right) \\
&\le C_{L,T,b,\sigma,\sigma^0} \left[ \left( 1 + \left\| \sup_{0 \le t \le T} |X_t| \right\|_p \right) (t-s) + \left( 1 + \left\| \sup_{0 \le t \le T} |X_t| \right\|_p \right)^2 \sqrt{t-s} \right] \\
&\le \kappa \sqrt{t-s} \quad (\text{by Inequality (1.8)}).
\end{align*}
$$

where $\kappa$ is a positive constant depending on $p, L, T, b, \sigma, \sigma^0$, and $\|X_0\|_p$. Then we have

$$
\|X_t - X_s\|_p \le \kappa \sqrt{t-s}. \quad \square
$$

Proof of Lemma 4.1. The particles $(\bar{Y}_t^i)_{1\le i\le N}$ are copies of the continuous of the continuous extension $\bar{X}$ of the Euler scheme (1.5). From (4.1), for every $1 \le i \le N$, there exists a measurable function $\Psi^i$ on $\mathbb{R}^d \times \mathcal{C}([0, T], \mathbb{R}^q) \times \mathcal{C}([0, T], \mathbb{R}^q)$ such that:

$$
\bar{Y}^i = \Psi^i \left( X_0^i, (W_t^i)_{0 \le t \le T}, (W_t^0)_{0 \le t \le T} \right).
$$

As the idiosyncratic noises and the initial random variables $(W^i, X_0^i)_{1\le i\le N}$ are i.i.d. by definition, the
particles $(\bar{Y}^i)_{1\le i\le N}$ are identically distributed and independent conditionally to $W^0$.
$\square$

References

[1] S. Ahuja. “Wellposedness of mean field games with common noise under a weak monotonicity condition”. In: *SIAM J. Control Optim.* 54.1 (2016), pp. 30–48. ISSN: 0363-0129.
[2] F. Antonelli and A. Kohatsu-Higa. "Rate of convergence of a particle method to the solution of the McKean-Vlasov equation". In: *Ann. Appl. Probab.* 12.2 (2002), pp. 423–476. ISSN: 1050-5164.

[3] J. Baladron et al. "Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons". In: *J. Math. Neurosci.* 2 (2012), Art. 10, 50.

[4] S. Biswas et al. "An explicit Milstein-type scheme for interacting particle systems and McKean-Vlasov SDEs with common noise and non-differentiable drift coefficients". In: *Ann. Appl. Probab.* 34.2 (2024), pp. 2326–2363. ISSN: 1050-5164.

[5] M. Bossy, J.-F. Jabir, and K. Martínez Rodríguez. "Instantaneous turbulent kinetic energy mod-
elling based on Lagrangian stochastic approach in CFD and application to wind energy". In: *J. Comput. Phys.* 464 (2022), Paper No. 110929, 29. ISSN: 0021-9991.

[6] M. Bossy and D. Talay. "A stochastic particle method for the McKean-Vlasov and the Burgers equation". In: *Math. Comp.* 66.217 (1997), pp. 157–192. ISSN: 0025-5718.

[7] P. Cardaliaguet and C.-A. Lehalle. "Mean field game of controls and an application to trade crowd-ing". In: *Math. Financ. Econ.* 12.3 (2018), pp. 335–363. ISSN: 1862-9679.

[8] R. Carmona and F. Delarue. *Probabilistic theory of mean field games with applications.* I. Vol. 83. Mean field FBSDEs, control, and games. Springer, Cham, 2018, pp. xxv+713.

[9] R. Carmona and F. Delarue. *Probabilistic theory of mean field games with applications.* II. Vol. 84. Mean field games with common noise and master equations. Springer, Cham, 2018, pp. xxiv+697.

[10] R. Carmona, F. Delarue, and D. Lacker. "Errata: Mean field games with common noise". In: *Ann. Probab.* 48.5 (2020), pp. 2644–2646. ISSN: 0091-1798.

[11] J.-F. Chassagneux and G. Pagès. "Computing the invariant distribution of McKean-Vlasov SDEs by ergodic simulation". In: *arXiv preprint arXiv:2406.13370* (2024).

[12] X. Chen and G. dos Reis. "A flexible split-step scheme for solving McKean-Vlasov stochastic dif-
ferential equations". In: *Appl. Math. Comput.* 427 (2022), Paper No. 127180, 23. ISSN: 0096-3003.

[13] F. Delarue. "Restoring uniqueness to mean-field games by randomizing the equilibria". In: *Stoch. Partial Differ. Equ. Anal. Comput.* 7.4 (2019), pp. 598–678. ISSN: 2194-0401.

[14] F. Delarue and R. Foguen Tchuendom. "Selection of equilibria in a linear quadratic mean-field game". In: *Stochastic Process. Appl.* 130.2 (2020), pp. 1000–1040. ISSN: 0304-4149.

[15] N. Fournier and A. Guillin. "On the rate of convergence in Wasserstein distance of the empirical measure". In: *Probab. Theory Related Fields* 162.3-4 (2015), pp. 707–738. ISSN: 0178-8051.

[16] P. Jameson Graber. "Linear quadratic mean field type control and mean field games with common noise, with application to production of an exhaustible resource". In: *Appl. Math. Optim.* 74.3 (2016), pp. 459–486. ISSN: 0095-4616.

[17] M. Hoffmann and Y. Liu. "A statistical approach for simulating the density solution of a McKean-Vlasov equation". working paper or preprint. 2023.

[18] M. Kac. "Foundations of kinetic theory". In: *Proceedings of The third Berkeley symposium on mathematical statistics and probability*. Vol. 3. 600. 1956, pp. 171–197.

[19] D. Lacker. "Hierarchies, entropy, and quantitative propagation of chaos for mean field diffusions". In: *Probab. Math. Phys.* 4.2 (2023), pp. 377–432. ISSN: 2690-0998.

[20] D. Lacker. *Mean Field Games and interacting particle systems*.

[21] D. Lacker. "On a strong form of propagation of chaos for McKean-Vlasov equations". In: *Electron. Commun. Probab.* 23 (2018), Paper No. 45, 11.
[22] D. Lacker and K. Webster. "Translation invariant mean field games with common noise". In: *Electron. Commun. Probab.* 20 (2015), no. 42, 13.

[23] J.-M. Lasry and P.-L. Lions. "Mean-field games with a major player". In: *Comptes Rendus Mathématique* 356.8 (2018), pp. 886-890. ISSN: 1631-073X.

[24] Y. Liu. "Particle method and quantization-based schemes for the simulation of the McKean-Vlasov equation". In: *ESAIM Math. Model. Numer. Anal.* 58.2 (2024), pp. 571-612. ISSN: 2822-7840.

[25] N. Martzel and C. Aslangul. "Mean-field treatment of the many-body Fokker-Planck equation". In: *J. Phys. A* 34.50 (2001), pp. 11225-11240. ISSN: 0305-4470.

[26] H. P. McKean Jr. "A class of Markov processes associated with nonlinear parabolic equations". In: *Proc. Nat. Acad. Sci. U.S.A.* 56 (1966), pp. 1907-1911. ISSN: 0027-8424.

[27] G. Pagès. *Numerical Probability*. Springer, 2018. ISBN: 978-3-319-90276-0.

[28] G. dos Reis, S. Engelhardt, and G. Smith. "Simulation of McKean-Vlasov SDEs with super-linear growth". In: *IMA J. Numer. Anal.* 42.1 (2022), pp. 874-922. ISSN: 0272-4979.

[29] X. Ren and D. Firoozi. "Risk-Sensitive Mean Field Games with Common Noise: A Theoretical Study with Applications to Interbank Markets". In: *arXiv preprint arXiv:2403.03915* (2024).

[30] A.-S. Sznitman. "Topics in propagation of chaos". In: *École d'Été de Probabilités de Saint-Flour XIX—1989*. Vol. 1464. Springer, Berlin, 1991, pp. 165-251.