# CONVERGENCE OF PARTICLES AND TREE BASED SCHEME FOR SINGULAR FBSDES

JEAN-FRANÇOIS CHASSAGNEUX<fnref n="1" />,* AND MOHAN YANG<fnref n="2" />

**Abstract.** We study an implementation of the theoretical splitting scheme introduced in Chassagneux and Yang (J. Comput. Phys. 468 (2022) 111459) for singular FBSDEs (Carmona and Delarue, Probab. Theory Relat. Fields 157 (2013) 333–388) and their associated quasi-linear degenerate PDEs. The fully implementable algorithm is based on particles approximation of the transport operator and tree like approximation of the diffusion operator appearing in the theoretical splitting. We prove the convergence with a rate of our numerical method under some reasonable conditions on the coefficients functions. This validates a posteriori some numerical results obtained in Chassagneux and Yang (J. Comput. Phys. 468 (2022) 111459). We conclude the paper with a numerical section presenting various implementations of the algorithm and discussing their efficiency in practice.

**Mathematics Subject Classification.** 65C30, 65C20, 65M12, 60H35.

Received March 6, 2023. Accepted December 30, 2024.

## 1. INTRODUCTION

In this work, we study the theoretical and numerical convergence of a particles and tree based scheme for singular FBSDEs. The class of singular Forward Backward Stochastic Differential Equations (FBSDE) is motivated by application in the modeling of carbon markets [4, 8, 10]. In these models, the FBSDE solution represents the equilibrium price of an allowance (to emit one ton of CO₂) in a perfect market. Let (Ω, F, P) be a stochastic basis supporting a d-dimensional Brownian motion W and T > 0 a terminal time. We denote by F := (F_t)_{t≥0} the filtration generated by the Brownian motion (augmented). The singular FBSDE system, with solution (P_t, E_t, Y_t, Z_t)_{0≤t≤T}, has the following form:

$$ \left\{ \begin{aligned} dP_t &= b(P_t) dt + \sigma(P_t) dW_t \\ dE_t &= \mu(P_t, Y_t) dt \\ dY_t &= Z_t dW_t \end{aligned} \right. . \qquad (1.1) $$

The terminal condition is typically given by

$$ Y_T = \phi(E_T) \text{ where } \phi(e) = \mathbf{1}_{e \ge \Lambda}, \qquad (1.2) $$

*Keywords and phrases.* Singular FBSDEs, probabilistic methods, splitting scheme, particle method.

<footnote n="1">¹ ENSAE-CREST and Institut Polytechnique de Paris, Palaiseau, France.</footnote>

<footnote n="2">² UFR de Mathématiques & LPSM, Université Paris Cité, Bâtiment Sophie Germain, 8 place Aurélie Nemours, 75013 Paris, France.</footnote>

<footnote>*Corresponding author: chassagneux@ensae.fr</footnote>
for some $\Lambda > 0$. Existence and uniqueness to the above FBSDE is not straightforward has it is fully coupled (the backward process $Y$ appears in the coefficient of the forward process $E$), it is degenerate in the forward direction $E$ and the terminal condition is discontinuous. However, assuming some structural conditions, see Theorem 2.1 below for a precise statement, Carmona and Delarue [3] managed to prove the wellposedness of such singular FBSDEs. They also show that the following Markovian representation holds true for the $Y$-process, namely:

$$ Y_t := \mathcal{V}(t, P_t, E_t), \quad 0 \le t < T, \qquad (1.3) $$

where the measurable function $\mathcal{V}: [0, T] \times \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$ is classically named the decoupling field.

The decoupling field $\mathcal{V}$ associated to (1.1) is a weak solution to

$$ \partial_t u + \partial_e(\mathfrak{M}(p, u)) + \partial_p u b(p) + \frac{1}{2} \mathrm{Tr}[\sigma(p)\sigma^\top(p)\partial_{pp}^2 u] = 0 \quad \text{and} \quad u(T, p, e) = \phi(e), \qquad (1.4) $$

with

$$ \mathfrak{M}(p, y) := \int_{0}^{y} \mu(p, v) \, dv, \quad 0 \le y \le 1, \qquad (1.5) $$

where $\partial_p$ denotes the Jacobian with respect to $p$, $\partial_t$ the time derivative, $\partial_e$ the derivative with respect to the $e$ variable, $\top$ is the transpose and $\tilde{c}_{pp}^2$ is the matrix of second derivatives with respect to the $p$ variable.

The algorithm we study is based on a splitting scheme designed for this kind of equation and introduced in [6]. Though it could be defined at a general level for the PDE, we should note that the convergence proof relies on the FBSDE setting. Indeed, under appropriate conditions, the FBSDE appears as well-posed random characteristics for the PDE (1.4). The splitting scheme involves the iteration on a discrete time grid of the composition of two operators: a transport operator $\mathcal{T}$ (in this case the $P$ variable is fixed) and a diffusion operator $\mathcal{D}$ (in this case the $E$ variable is fixed). Precise definitions are given in the next section. From [6], we already know that the theoretical splitting scheme is convergent with a rate one half. This result is obtained under the same structural conditions which guarantee the well-posedness of the FBSDE. In [6], the splitting scheme is then implemented using various finite difference approximations of the transport operator and a non-linear regression using Deep Neural Networks for the diffusion operator. An alternative scheme ([6], Sect. 3.3.1) is also suggested to verify the convergence result of the main numerical procedure. The goal of this paper is to analyse precisely this alternative scheme by proving its convergence with a rate, under some regularity conditions and assuming a constant volatility coefficient in the $P$ dynamics.

The rest of the paper is organized as follows. In the next section, we recall the theoretical setting for the study of singular FBSDEs and recall the theoretical splitting scheme. We also introduce the framework used to study the convergence of the numerical methods and give some properties of the solution $\mathcal{V}$ in this framework. In Section 3, we present the numerical algorithm and some of its variants that will be used for numerical tests. We also state our main convergence result in Theorem 3.1. Section 4 is dedicated to the proof of the convergence by studying precisely all sources of errors. The last section presents numerical results for the various schemes introduced, illustrating the theoretical convergence.

**Notations.** In the following we will use the following spaces:

Let $\mathbb{H} = (\mathcal{H}_t)_{t \ge 0}$ denotes a generic right-continuous filtration,

- For fixed $0 \le a < b < +\infty$ and $I = [a, b]$ or $I = [a, b) \cup (b, +\infty)$, $S^{2,k}(I, \mathbb{H})$ is the set of $\mathbb{R}^k$-valued càdlàg<fnref n="1" /> $\mathcal{H}_t$-adapted processes $Y$, s.t.

$$ \|Y\|_{S^2}^2 := \mathbb{E} \left[ \sup_{t \in I} |Y_t|^2 \right] < \infty. $$

Note that we may omit the dimension and the terminal date in the norm notation as this will be clear from the context. $S_c^{2,k}(I, \mathbb{H})$ is the subspace of processes with continuous sample paths.

<footnote n="1">¹French acronym for right continuous with left limits.</footnote>
- For fixed $0 \le a < b < +\infty$, and $I = [a, b]$, we denote by $\mathcal{H}^{2,k}(I, \mathbb{H})$ the set of $\mathbb{R}^k$-valued progressively measurable processes $Z$, such that

$$ \|Z\|_{\mathcal{H}^2}^2 := \mathbb{E}\left[\int_I |Z_t|^2 dt\right] < \infty. $$

We denote by $\mathcal{P}(\mathbb{R})$ the space of probability measures on $\mathbb{R}$ and $\mathcal{P}_q(\mathbb{R})$, $q \ge 1$ the subset of probability measure which have a $q$-moment finite.

We denote by $\mathcal{I}$ the space of cumulative distribution function (CDF) namely functions $\theta$ given by

$$ \mathbb{R} \ni x \mapsto \theta(x) = \mu((-\infty, x]) \in [0, 1] $$

for some $\mu \in \mathcal{P}(\mathbb{R})$.

The set $\mathcal{C}^{1,2,1}(I \times \mathbb{R}^d \times \mathbb{R})$ (resp. $\mathcal{C}^{1,2,2}(I \times \mathbb{R}^d \times \mathbb{R})$) with $I = [0, T)$ or $I = [0, T]$, is the set of functions $f: I \times \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$ continuously differentiable: once in their first variable, twice in their second variable and once (resp. twice) in their third variable.

## 2. REVIEW OF THEORETICAL RESULT FOR SINGULAR FBSDES

In this section, we first recall the main properties of solution to (1.1). We then define the theoretical splitting scheme studied in [6] and used to design our implemented algorithms introduced in the next section. Finally, we consider a smooth framework for the solution of (1.4) and obtain some properties useful to prove the convergence of the numerical schemes.

### 2.1. Well-posedness of Singular FBSDES

We first define a class of functions to which the solution belongs.

**Definition 2.1.** Let $\mathcal{K}$ be the class of functions $\Phi: \mathbb{R}^d \times \mathbb{R} \to [0, 1]$ such that $\Phi$ is $L_\phi$-Lipschitz in the first variable for some $L_\Phi > 0$, namely

$$ |\Phi(p, e) - \Phi(p', e)| \le L_\Phi |p - p'| \quad \text{for all} \quad (p, p', e) \in \mathbb{R}^d \times \mathbb{R}^d \times \mathbb{R}, \qquad (2.1) $$

and moreover satisfying, for all $p \in \mathbb{R}^d$,

$$ \Phi(p, e) = \nu(p, (-\infty, e]) \text{ where } \nu(p, \cdot) \in \mathcal{P}(\mathbb{R}). \qquad (2.2) $$

From the above definition, we have that $\Phi(p, \cdot) \in \mathcal{I}$ and is thus a non-decreasing, right continuous function and satisfies, for all $p \in \mathbb{R}^d$,

$$ \lim_{e \to -\infty} \Phi(p, e) = 0 \quad \text{and} \quad \lim_{e \to -\infty} \Phi(p, e) = 1. \qquad (2.3) $$

We now recall the existence and uniqueness result for the singular FBSDE. It is also key to define the theoretical splitting scheme. This result is obtained for the following class of admissible coefficient functions.

**Definition 2.2.** Let $\mathcal{A}$ be the class of functions $B: \mathbb{R}^d \to \mathbb{R}^d$, $\Sigma: \mathbb{R}^d \to \mathcal{M}_d$, $F: \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$ which are $L$-Lipschitz continuous functions. Moreover, $F$ is strictly decreasing in $y$ and satisfies, for all $p \in \mathbb{R}^d$,

$$ \ell_1|y - y'|^2 \le (y - y')(F(p, y') - F(p, y)) \le \ell_2|y - y'|^2, \qquad (2.4) $$

where $L$, $\ell_1$ and $\ell_2$ are positive constants.

The well-posedness result in this setting reads as follows.
**Theorem 2.1** ([3], Prop. 2.10, [8], Prop. 3.2). Let $\tau > 0$, $(B, \Sigma, F) \in \mathcal{A}$ and $\Phi \in \mathcal{K}$.

Given any initial condition $(t_0, p, e) \in [0, \tau) \times \mathbb{R}^d \times \mathbb{R}$, there exists a unique progressively measurable 4-tuple of processes $(P_t^{t_0,p,e}, E_t^{t_0,p,e}, Y_t^{t_0,p,e}, Z_t^{t_0,p,e})_{t_0 \le t \le \tau} \in S_c^{2,d}([t_0, \tau], \mathbb{F}) \times S_c^{2,1}([t_0, \tau], \mathbb{F}) \times S_c^{2,1}([t_0, \tau), \mathbb{F}) \times \mathcal{H}^{2,d}([t_0, \tau], \mathbb{F})$ satisfying the dynamics

$$ 
\begin{aligned}
dP_t^{t_0,p,e} &= B(P_t^{t_0,p,e}) dt + \Sigma(P_t^{t_0,p,e}) dW_t, & P_t^{t_0,p,e} &= p \in \mathbb{R}^d, \\
dE_t^{t_0,p,e} &= F(P_t^{t_0,p,e}, Y_t^{t_0,p,e}) dt, & E_t^{t_0,p,e} &= e \in \mathbb{R}, \\
dY_t^{t_0,p,e} &= Z_t^{t_0,p,e} dW_t,
\end{aligned} 
\quad (2.5) $$

and such that

$$ \mathbb{P} \left[ \Phi_{-}(P_{\tau}^{t_0,p,e}, E_{\tau}^{t_0,p,e}) \le \lim_{t \uparrow \tau} Y_{t}^{t_0,p,e} \le \Phi_{+}(P_{\tau}^{t_0,p,e}, E_{\tau}^{t_0,p,e}) \right] = 1. \quad (2.6) $$

The unique decoupling field defined by

$$ [0, \tau) \times \mathbb{R}^d \times \mathbb{R} \ni (t_0, p, e) \to w(t_0, p, e) = Y_{t_0}^{t_0, p, e} \in \mathbb{R} $$

is continuous and satisfies

(1) For any $t \in [0, \tau)$, the function $w(t, \cdot, \cdot)$ is $1/(l_1(\tau - t))$-Lipschitz continuous with respect to $e$.

(2) For any $t \in [0, \tau)$, the function $w(t, \cdot, \cdot)$ is C-Lipschitz continuous with respect to $p$, where $C$ is a constant depending on $L$, $\tau$ and $L_\phi$ only.

(3) Given $(p, e) \in \mathbb{R}^d \times \mathbb{R}$, for any family $(p_t, e_t)_{0 \le t < \tau}$ converging to $(p, e)$ as $t \uparrow \tau$, we have

$$ \Phi_{-}(p, e) \le \liminf_{t \to \tau} w(t, p_t, e_t) \le \limsup_{t \to \tau} w(t, p_t, e_t) \le \Phi_{+}(p, e). \quad (2.7) $$

(4) For any $t \in [0, \tau)$, the function $w(t, \cdot, \cdot) \in \mathcal{K}$.

The previous result leads us to define a non-linear operator associated to (1.1) under the following assumption, that will hold true for the rest of the paper:

**Standing assumptions.** The coefficients $(b, \sigma, \mu)$ of the singular FBSDE (1.1) belong to $\mathcal{A}$.

**Definition 2.3.** We define the semi-group of operator $(\Theta_h)_{h>0}$ by

$$ (0, \infty) \times \mathcal{K} \ni (h, \psi) \mapsto \Theta_h(\psi) = v(0, \cdot) \in \mathcal{K} \quad (2.8) $$

where $v$ is the decoupling field given in Theorem 2.1 with parameters $\tau = h$, $B = b$, $\Sigma = \sigma$, $F = \mu$ and $\Phi = \psi$.

Let us now be given a discrete time grid of $[0, T]$:

$$ \pi := \{t_0 := 0 \le \dots \le t_n \le \dots t_N := T\}, \quad (2.9) $$

for $N \ge 1$. For ease of presentation, we assume that the time grid $\pi$ is equidistant and thus, for $0 \le n \le N - 1$,

$$ t_{n+1} - t_n = \frac{T}{N} =: \mathfrak{h}. \quad (2.10) $$

Using Definition 2.3, we observe that

$$ \mathcal{V}(0, \cdot) := \Theta_T(\phi) = \Theta_{t_1-t_0} \circ \dots \circ \Theta_{t_N-t_{N-1}}(\phi). \quad (2.11) $$

We recall now a result that arises from the proof of the previous Theorem 2.1, see [8].
**Corollary 2.1** (Approximation result). Let $\tau > 0$, $(B, \Sigma, F) \in \mathcal{A}$ and $\Phi \in \mathcal{K}$. Let $(\Phi^k)_{k \ge 0}$ be a sequence of smooth functions belonging to $\mathcal{K}$ and converging pointwise towards $\Phi$ as $k$ goes to $+\infty$. For $\eta > 0$, consider then $w^{\eta,k}$ the solution to:

$$ \partial_t u + F(p, u)\partial_e u + \mathcal{L}_p u + \frac{1}{2}\eta^2(\partial_{ee}^2 u + \Delta_{pp} u) = 0 \quad \text{and} \quad u(\tau, \cdot) = \Phi^k \qquad (2.12) $$

where $\Delta_{pp}$ is the Laplacian with respect to $p$, and $\mathcal{L}_p$ is the operator

$$ \mathcal{L}_p(\varphi)(t, p, e) = \partial_p \varphi(t, p, e)B(p) + \frac{1}{2}\mathrm{Tr}[A(p)\partial_{pp}^2](\varphi)(t, p, e), \qquad (2.13) $$

with $A = \Sigma\Sigma^\top$. For later use, we also define $\mathcal{L}^\eta := \mathcal{L}_p + \frac{1}{2}\eta^2(\partial_{ee}^2 + \Delta_{pp})$. Then the functions $w^{\eta,k}$ are $\mathcal{C}^{1,2,2}([0, \tau] \times \mathbb{R}^d \times \mathbb{R})$ and $\lim_{k \to \infty} \lim_{\eta \to 0} w^{\eta,k} = w$ where the convergence is locally uniform in $[0, \tau) \times \mathbb{R}^d \times \mathbb{R}$. Moreover, for all $k, \eta$, $w^{\eta,k}(t, \cdot) \in \mathcal{K}$.

For later use, we introduce the system of FBSDEs associated to (2.12). We thus consider $(W, \tilde{B}, \tilde{\tilde{B}})$ to be independent Brownian motions with dimension respectively $d, d$ and one, and denote $\mathcal{G} = (\mathcal{G}_t)_{t \ge 0}$ the filtration they generate (augmented). For $t_0 \in [0, T)$ and $(p, e) \in \mathbb{R}^d \times \mathbb{R}$, let

$$ (P^\eta, E^{\eta,k}) \in \mathcal{S}_c^{2,d}([t_0, T], \mathcal{G}) \times \mathcal{S}_c^{2,1}([t_0, T], \mathcal{G}) $$

and

$$ (Y^{\eta,k}, Z^{\eta,k}, \tilde{Z}^{\eta,k}, \tilde{\tilde{Z}}^{\eta,k}) \in \mathcal{S}_c^{2,d}([t_0, T], \mathcal{G}) \times \mathcal{H}^{2,d}([t_0, T], \mathcal{G}) \times \mathcal{H}^{2,d}([t_0, T], \mathcal{G}) \times \mathcal{H}^{2,1}([t_0, T], \mathcal{G}) $$

be the solution to

$$ \left\{ \begin{array}{l} dP_t^\eta = B(P_t^\eta) dt + \Sigma(P_t^\eta) dW_t + \eta d\tilde{B}_t, \quad P_{t_0}^\eta = p, \\ dE_t^{\eta,k} = F(P_t^\eta, Y_t^{\eta,k}) dt + \eta d\tilde{\tilde{B}}_t, \quad E_{t_0}^\eta = e, \\ dY^{\eta,k} = Z^{\eta,k} dW_t + \tilde{Z}^{\eta,k} d\tilde{B}_t + \tilde{\tilde{Z}}^{\eta,k} d\tilde{\tilde{B}}_t. \end{array} \right. \qquad (2.14) $$

Observe that $Y_t^{\eta,k} = w^{\eta,k}(t, P_t^\eta, E_t^{\eta,k})$ for $t \in [t_0, T]$, in particular $Y_T^{\eta,k} = \Phi^k(P_T, E_T)$. Let us note that the above processes depend upon the initial condition $(t_0, p, e)$ but we slightly abuse the notation in this regard by omitting to indicate the initial condition.

## 2.2. A theoretical splitting scheme

The numerical algorithm is based on a splitting method for the quasilinear PDE (1.4), which has been introduced in [6]. We now recall the splitting approach at a theoretical level.

We first introduce the transport step where the diffusion part is frozen.

**Definition 2.4** (Transport step). We set

$$ (0, \infty) \times \mathcal{K} \ni (h, \psi) \mapsto \mathcal{T}_h(\psi) = \tilde{v}(0, \cdot) \in \mathcal{K}, $$

where $\tilde{v}$ is the decoupling field defined in Theorem 2.1 with parameters $\tau = h$, $B = 0$, $\Sigma = 0$, $F = \mu$ and terminal condition $\Phi = \psi$.

In the definition above, $\tilde{v}(\cdot)$, defined on $[0, h] \times \mathbb{R}^d \times \mathbb{R}$, is the unique entropy solution to

$$ \partial_t w + \partial_e(\mathfrak{M}(p, w)) = 0 \quad \text{and} \quad \tilde{v}(h, \cdot) = \psi, \qquad (2.15) $$

recall (1.4).
**Remark 2.1.** For later use in the convergence analysis, we introduce a small modification of the previous operator that acts on function in $\mathcal{J}$ instead of $\mathcal{K}$. Namely, we set

$$ (0, \infty) \times \mathbb{R}^d \times \mathcal{J} \ni (h, p, \theta) \mapsto \tilde{T}_h(p, \theta) = \hat{v}_p(0, \cdot) \in \mathcal{J} $$

where $\hat{v}_p(0, \cdot)$ is the unique entropy solution to

$$ \partial_t w + \partial_e(\mathfrak{M}(p, w)) = 0 \quad \text{and} \quad \hat{v}_p(h, \cdot) = \theta. \qquad (2.16) $$

One observes that $\mathcal{T}_h(\psi)(p, \cdot) = \tilde{T}_h(p, \psi(p, \cdot))$.

Let us mention that, from e.g. Theorem 1.1 of [11], we have the following stability result:

$$ \int_{\mathbb{R}} |\tilde{T}_h(p, \psi_1(p, \cdot))(e) - \tilde{T}_h(p, \psi_2(p, \cdot))(e)| \, de \le \int_{\mathbb{R}} |\psi_1(p, e) - \psi_2(p, e)| \, de, \qquad (2.17) $$

for $\psi_1, \psi_2 \in \mathcal{K}$.

We now introduce the diffusion step, where conversely, the E-process is frozen to its initial value.

**Definition 2.5 (Diffusion step).** We set

$$ (0, \infty) \times \mathcal{K} \ni (h, \psi) \mapsto \mathcal{D}_h(\psi) = \bar{v}(0, \cdot) \in \mathcal{K} $$

where $\bar{v}(0, \cdot)$ is the decoupling field in Theorem 2.1 with parameters $\tau = h$, $B = b$, $\Sigma = \sigma$, $F = 0$ and terminal condition $\Phi = \psi$.

Observe that, for $t \in [0, h)$,

$$ \bar{v}(t, p, e) = \mathbb{E}[\psi(P_h^{t,p}, e)], \quad (p, e) \in \mathbb{R}^d \times \mathbb{R} \quad \text{and} \quad \bar{v}(t, \cdot) \in \mathcal{K}. \qquad (2.18) $$

We can now define the theoretical scheme on $\pi$ by a backward induction. Since

$$ \mathcal{V}(0, \cdot) = \Theta_T(\phi) = \prod_{0 \le n < N} \Theta_{t_{n+1}-t_n}(\phi), \qquad (2.19) $$

the main point of the scheme is to replace one step of $\Theta$ by one step of $\mathcal{T} \circ \mathcal{D}$. This leads to the following.

**Definition 2.6 (Theoretical splitting scheme).** We set

$$ (0, \infty) \times \mathcal{K} \ni (h, \psi) \mapsto \mathcal{S}_h(\psi) := \mathcal{T}_h \circ \mathcal{D}_h(\psi) \in \mathcal{K}. $$

For $n \le N$, we denote by $u_n^\pi$ the solution of the following backward induction on $\pi$:

- for $n = N$, set $u_N^\pi := \phi$,

- for $n < N$, $u_n^\pi = S_{t_{n+1}-t_n}(u_{n+1}^\pi)$.

The $(u_n^\pi)_{0 \le n \le N}$ stands for the approximation of the decoupling field $\mathcal{V}(t, \cdot)$ for $t \in \pi$.

It is proven in Proposition 2.1 of [6] a key result concerning the scheme's truncation error.

**Proposition 2.1 (truncation error).** Under our standing assumptions on the coefficients $(\mu, b, \sigma)$, the following holds, for $\psi \in \mathcal{K}$:

$$ \int_{\mathbb{R}} |\mathcal{S}_h(\psi)(p, e) - \Theta_h(\psi)(p, e)| \, de \le C_{L_\psi} (1 + |p|^2)^{\frac{3}{2}}, \qquad (2.20) $$

for $p \in \mathbb{R}^d$, $h > 0$.

This allows in particular to obtain the convergence with a rate $\frac{1}{2}$ of the theoretical splitting scheme to the decoupling field $\mathcal{V}$, namely

$$ \int_{\mathbb{R}} |\mathcal{V}(0, p, e) - u_0^\pi(p, e)| \, de \le C(1 + |p|^2) \sqrt{\bar{h}}, \quad p \in \mathbb{R}^d, $$

for a positive constant $C$, see Theorem 2.2 of [6].
## 2.3. Smooth setting for convergence results

We now consider a more restricted framework. In this setting, we give further properties of the value function $V$ that will be used to prove our convergence results. Indeed, even if the scheme can be defined in a quite general setting, the convergence with a rate is obtained in a more regular framework.

Therefore, the main assumption we shall use is the following, see Remark 2.2 below.

**Assumption 2.1.** The decoupling field $V$ is a $C^{1,2,1}([0,T) \times \mathbb{R}^d \times \mathbb{R})$ solution to

$$ \partial_t V + \mu(p, V)\partial_e V + b(p)\partial_p V + \frac{1}{2}\bar{\sigma}^2 \Delta_p V = 0. \quad (2.21) $$

In particular, the function $\sigma(\cdot)$ is constant $\sigma(\cdot) := \bar{\sigma}I_d$ with $\bar{\sigma} > 0$ and $I_d$ the $d \times d$ identity matrix, recall (1.4). Moreover, the function $b$ is twice differentiable with bounded and Lipschitz derivatives and the function $\mu$ is differentiable with Lipschitz and bounded derivatives.

We first list some further properties of the function $V$ in the regular setting of the previous assumption.

**Proposition 2.2.** Under Assumption 2.1, the function $V$ satisfies, for $(t, p, e) \in [0, T) \times \mathbb{R}^d \times \mathbb{R}$,

$$ 0 \le \partial_e V(t, p, e) \le \frac{C}{T-t} \quad \text{and} \quad |\partial_p V(t, p, e)| \le C. \quad (2.22) $$

Moreover, for $(t, p) \in [0, T) \times \mathbb{R}^d$,

$$ \int_{\mathbb{R}} |\partial_p V(t, p, e)| \, de \le C(T-t), \quad (2.23) $$

and, thus, for $p' \in \mathbb{R}^d$,

$$ \int_{\mathbb{R}} |V(t, p, e) - V(t, p', e)| \, de \le C(T-t)|p - p'|. \quad (2.24) $$

*Proof.* The estimates (2.22) are obtained directly from the property of $V$ given in Theorem 2.1. We now study the $L^1$-control.

Let $V^{\eta,k}$ be the smooth approximation given in Corollary 2.1 to $V$. We consider also the associated FBSDE $(P^{\eta}, E^{\eta,k}, Y^{\eta,k})$ in (2.14). We adapt the computations in the proof of Lemma 3.8 in [8], to obtain the equation corresponding to equation (3.35) there. In our simpler framework ($r = 0, \phi$ does not depend on $p$), it reads

$$ \partial_p V^{\eta,k}(t_0, p, e) = \mathbb{E} \left[ \int_{t_0}^T \partial_p \mu(P_t^\eta, Y_t^{\eta,k}) \partial_p P_t^\eta \partial_e E_t^{\eta,k} \partial_e V^{\eta,k}(t, P_t^\eta, E_t^{\eta,k}) dt \right] \quad (2.25) $$

where $(\partial_p P^{\eta}, \partial_e E^{\eta,k})$ are the tangent processes associated to $(P^{\eta}, E^{\eta,k})$. They are given by

$$ \partial_e E_t^{\eta,k} = 1 + \int_{t_0}^t \partial_y \mu(P_s^\eta, V^{\eta,k}(s, P_s^\eta, E_s^{\eta,k})) \partial_e V^{\eta,k}(s, P_s^\eta, E_s^{\eta,k}) \partial_e E_s^{\eta,k} \, ds, \quad (2.26) $$

$$ (\partial_{p^{\prime\ell}} P_t^{\eta})^{\ell} = \mathbf{1}_{\{\ell=\ell'\}} + \int_{t_0}^{t} \partial_p b^{\ell}(P_s^{\eta}) \partial_{p^{\prime\ell}} P_s^{\eta} \, ds \quad \text{for } \ell, \ell' \in \{1, \dots, d\}. \quad (2.27) $$

Observe, in particular, that

$$ \partial_e E_t^{\eta,k} = \exp\left(\int_{t_0}^t \partial_y \mu(P_s^\eta, V^{\eta,k}(s, P_s^\eta, E_s^{\eta,k})) \partial_e V^{\eta,k}(s, P_s^\eta, E_s^{\eta,k}) \, ds\right), \quad t \in [t_0, T]. $$
Moreover, it is well known in our Lipschitz setting that

$$ \mathbb{E} \left[ \sup_{t \in [t_0, T]} |\hat{\partial}_p P_t^{\eta}|^{\kappa} \right] \le C_{\kappa}, \qquad (2.28) $$

for $\kappa \ge 1$.

For a positive real number $R$, we obtain from (2.25) that

$$ \int_{-R}^{R} |\hat{\partial}_p \mathcal{V}^{\eta,k}(t_0, p, e)| \, de \le \mathbb{E} \left[ \int_{t_0}^{T} \int_{-R}^{R} |\hat{\partial}_p \mu(P_t^{\eta}, Y_t^{\eta,k}) \hat{\partial}_p P_t^{\eta}| |\hat{\partial}_e[\mathcal{V}^{\eta,k}(t, P_t^{\eta}, E_t^{\eta,k})]| \, dt \, de \right] $$

recalling that $\hat{\partial}_e E^{\eta,k}$ and $\hat{\partial}_e \mathcal{V}^{\eta,k}(\cdot)$ are non-negative. Moreover, since $\mathcal{V}^{\eta,k}$ is bounded and $\mu$ Lipschitz continuous, we compute

$$ \int_{-R}^{R} |\hat{\partial}_p \mathcal{V}^{\eta,k}(t_0, p, e)| \, de \le 2 \|\mathcal{V}^{\eta,k}\|_{\infty} \|\hat{\partial}_p \mu\|_{\infty} (T - t_0) \mathbb{E} \left[ \sup_{t \in [t_0, T]} |\hat{\partial}_p P_t^{\eta}| \right]. $$

Using (2.28), we eventually obtain

$$ \int_{-R}^{R} |\hat{\partial}_p \mathcal{V}^{\eta,k}(t_0, p, e)| \, de \le C(T - t_0). \qquad (2.29) $$

Since $\mathcal{V}^{\eta,k}(t_0, p + \delta, e) - \mathcal{V}^{\eta,k}(t_0, p, e) = \delta \int_0^1 \hat{\partial}_p \mathcal{V}^{\eta,k}(t_0, p + \lambda \delta, e) \, d\lambda$, we have

$$ \int_{-R}^{R} |\mathcal{V}^{\eta,k}(t_0, p + \delta, e) - \mathcal{V}^{\eta,k}(t_0, p, e)| \, de \le \delta \int_{0}^{1} \int_{-R}^{R} |\hat{\partial}_p \mathcal{V}^{\eta,k}(t_0, p + \lambda \delta, e)| \, de \, d\lambda, \qquad (2.30) $$

$$ \le C(T - t_0)\delta \qquad (2.31) $$

where we used (2.29) for the last inequality. The dominated convergence theorem leads to

$$ \int_{-R}^{R} \frac{1}{\delta} |\mathcal{V}(t_0, p + \delta, e) - \mathcal{V}(t_0, p, e)| \, de \le C(T - t_0), \qquad (2.32) $$

recalling Corollary 2.1. Taking limit as $\delta \to 0$, we obtain

$$ \int_{-R}^{R} |\hat{\partial}_p \mathcal{V}(t_0, p, e)| \, de \le C(T - t_0). \qquad (2.33) $$

The estimate (2.23) is then obtained by letting $R \to +\infty$ and invoking monotone convergence. $\square$

The study of the regularity of the decoupling field, in particular how to prove Assumption 2.1, is outside the scope of this paper and is left for further research. We, however, give now an example of model for which Assumption 2.1 is known to hold.

**Remark 2.2.** Let us consider the following *toy model*, which has been studied in [4]

$$ \begin{cases} dP_t = \bar{\sigma} \, dW_t, \quad P_0 = p, \\ dE_t = (P_t - Y_t) \, dt, \quad E_0 = e, \\ dY_t = Z_t \, dW_t. \end{cases} \qquad (2.34) $$

Here, $\bar{\sigma}$ is a positive constant, $W$ is a one-dimensional Brownian Motion and the terminal condition (in the weak sense of (2.6)) for $Y$ is given by $\phi(E_T)$, recall (1.2). By Theorem 2.1 above, it has a unique solution and
we denote by $v$ its decoupling field. In particular, we have $v(0, p, e) = Y_0$. Now, we introduce, as in [4], the following process

$$ \bar{E}_t = E_t + (T-t)P_t, \qquad (2.35) $$

and observe then that $(\bar{E}, Y)$ satisfies

$$ \left\{ \begin{array}{l} d\bar{E}_t = -Y_t dt + (T-t) dW_t \\ dY_t = Z_t dW_t \end{array} \right. \qquad (2.36) $$

with the same terminal condition since $\bar{E}_T = E_T$. Independently from (2.34), the system (2.36) is thoroughly studied in Section 4 of [4], this system has also a unique solution with a decoupling field $Y_0 = \bar{v}(0, \bar{E}_0)$, and from Proposition 6 in [4] again, $\bar{v}$ is $\mathcal{C}^{1,2}([0, T) \times \mathbb{R})$. We deduce that $v(0, p, e) = \bar{v}(0, e + Tp)$ and more generally $v(t, p, e) = \bar{v}(t, e + (T - t)p)$ for $t \in [0, T)$. This yields that $v$ is $\mathcal{C}^{1,2,2}([0, T) \times \mathbb{R}^d \times \mathbb{R})$.

Let us observe that Proposition 2.2 states estimates for the first order derivatives of $\mathcal{V}$. And, though $\mathcal{V}$ is assumed to be smooth in Assumption 2.1, nothing is said on the control of higher order derivatives. As we will see later on, the control of the second order derivative and third order derivative with respect to $p$, will be needed. This should not come as a surprise as we will approximate the $P$ process by a discrete in time and space process, see Section 3.1.1. We shall not make further assumption on the behavior of $\mathcal{V}$ and its derivatives. Instead, we introduce a proxy to $\mathcal{V}$ that will be used in the proof of convergence of our numerical algorithm.

Let then $\mathcal{V}^\epsilon$ be a smoothing of the decoupling field $\mathcal{V}$, namely: For $\epsilon \in (0, 1)$,

$$ \mathcal{V}^{\epsilon}(t, p, e) := \int_{\mathbb{R}^d} \mathcal{V}(t, p+q, e)\varphi_{\epsilon}(q) \, dq \quad \text{where} \quad \varphi_{\epsilon}(q) := \frac{1}{\epsilon^d} \varphi\left(\frac{q}{\epsilon}\right), \quad q \in \mathbb{R}^d, \qquad (2.37) $$

and with $\varphi(\cdot)$ a smooth compactly supported probability density function.

We now show that $\mathcal{V}^\epsilon$ satisfies the same PDE as $\mathcal{V}$ up to an error term that we can control. From the previous Assumption 2.1, we first observe that $\mathcal{V}^\epsilon$, defined in (2.37) belongs also to $\mathcal{C}^{1,2,1}([0, T) \times \mathbb{R}^d \times \mathbb{R})$ and that it satisfies the same bounds as $\mathcal{V}$, namely

$$ 0 \le \partial_t \mathcal{V}^\epsilon(t, p, e) \le \frac{C}{T-t}, \quad |\partial_t \mathcal{V}^\epsilon(t, p, e)| \le C \qquad (2.38) $$

$$ \text{and} \quad \int_{\mathbb{R}} |\partial_p \mathcal{V}^\epsilon(t, p, z)| \, dz \le C(T-t), \quad \forall (t, p, e) \in [0, T) \times \mathbb{R}^d \times \mathbb{R}. \qquad (2.39) $$

**Proposition 2.3.** Under Assumption 2.1, the function $\mathcal{V}^\epsilon$ satisfies on $[0, T) \times \mathbb{R}^d \times \mathbb{R}$,

$$ \partial_t \mathcal{V}^\epsilon + b(p) \partial_p \mathcal{V}^\epsilon + \frac{1}{2} \bar{\sigma}^2 \Delta_p \mathcal{V}^\epsilon + \mu(p, \mathcal{V}(t, p, e)) \partial_e \mathcal{V}^\epsilon = \theta^\epsilon(t, p, e) \qquad (2.40) $$

with

$$ |\theta^\epsilon(t, p, e)| \le C\epsilon \left( \int_{\mathbb{R}^d} |\partial_p \mathcal{V}(t, p+q, e)| \varphi_\epsilon(q) \, dq + \partial_e \mathcal{V}^\epsilon(t, p, e) \right). \qquad (2.41) $$

*Proof.* Recall that $\mathcal{V}$ satisfies (2.21), and that, for $(t, p, e) \in [0, T) \times \mathbb{R}^d \times \mathbb{R}$,

$$ \partial_t \mathcal{V}^\epsilon(t, p, e) = \int_{\mathbb{R}^d} \partial_t \mathcal{V}(t, p+q, e) \varphi_\epsilon(q) \, dq \quad \text{and} \quad \partial_{pp}^2 \mathcal{V}^\epsilon(t, p, e) = \int_{\mathbb{R}^d} \partial_{pp}^2 \mathcal{V}(t, p+q, e) \varphi_\epsilon(q) \, dq. \qquad (2.42) $$

We thus compute, by linearity of $\partial_t$, $\Delta_p$,

$$ \partial_t \mathcal{V}^\epsilon(t, p, e) + \frac{1}{2} \bar{\sigma}^2 \Delta_p \mathcal{V}^\epsilon(t, p, e) = - \int_{\mathbb{R}^d} b(p+q) \partial_p \mathcal{V}(t, p+q, e) \varphi_\epsilon(q) \, dq \qquad (2.43) $$
$$ - \int_{\mathbb{R}^d} \mu(p+q, \mathcal{V}(t, p+q, e)) \partial_e \mathcal{V}(t, p+q, e) \varphi_e(q) \, dq \quad (2.44) $$

which leads to

$$ \partial_t \mathcal{V}^\epsilon(t, p, e) + b(p) \partial_p \mathcal{V}^\epsilon(t, p, e) + \frac{1}{2} \bar{\sigma}^2 \Delta_p \mathcal{V}^\epsilon(t, p, e) \quad (2.45) $$

$$ + \mu(p, \mathcal{V}(t, p, e)) \partial_e \mathcal{V}^\epsilon(t, p, e) = \theta^\epsilon(t, p, e) \quad (2.46) $$

with

$$ \theta^\epsilon(t, p, e) := \int_{\mathbb{R}^d} \{b(p) - b(p+q)\} \partial_p \mathcal{V}(t, p+q, e) \varphi_\epsilon(q) \, dq \quad (2.47) $$

$$ + \int_{\mathbb{R}^d} \{\mu(p, \mathcal{V}(t, p, e)) - \mu(p+q, \mathcal{V}(t, p+q, e))\} \partial_e \mathcal{V}(t, p+q, e) \varphi_\epsilon(q) \, dq. \quad (2.48) $$

We compute, using the uniform Lipschitz property of $b$, $\mu$ and $\mathcal{V}$ in the $p$ variable and the non-negativity of $\partial_e \mathcal{V}$, that

$$ |\theta^\epsilon(t, p, e)| \le C \int_{\mathbb{R}^d} |q| |\partial_p \mathcal{V}(t, p+q, e)| \varphi_\epsilon(q) \, dq + C \int_{\mathbb{R}^d} |q| \partial_e \mathcal{V}(t, p+q, e) \varphi_\epsilon(q) \, dq \quad (2.49) $$

$$ \le C\epsilon \left( \int_{\mathbb{R}^d} |\partial_p \mathcal{V}(t, p+q, e)| \varphi_\epsilon(q) \, dq + \partial_e \mathcal{V}^\epsilon(t, p, e) \right) \quad (2.50) $$

where we use the fact that $\varphi$ has compact support to get the last inequality. $\square$

To conclude this section, we collect some useful control in the $L^1$-sense of the function $\mathcal{V}^\epsilon$ and its derivatives. Their proof is very classical and we skip it.

**Lemma 2.1.** Under Assumption 2.1, the followings hold, for $(t, p) \in [0, T) \times \mathbb{R}^d$,

$$ \epsilon \int_{\mathbb{R}} |\partial_{p_i p_j}^2 \mathcal{V}^\epsilon(t, p, e)| \, de + \epsilon^2 \int_{\mathbb{R}} |\partial_{p_i p_j p_k}^3 \mathcal{V}^\epsilon(t, p, e)| \, de \le C, \quad (2.51) $$

for $i, j, k \in \{1, \dots, d\}$. And thus, for $(t, p, p') \in [0, T) \times \mathbb{R}^d \times \mathbb{R}^d$,

$$ \int_{\mathbb{R}} |\partial_p \mathcal{V}^\epsilon(t, p, e) - \partial_p \mathcal{V}^\epsilon(t, p', e)| \, de + \epsilon \int_{\mathbb{R}} |\partial_{pp'}^2 \mathcal{V}^\epsilon(t, p, e) - \partial_{pp'}^2 \mathcal{V}^\epsilon(t, p', e)| \, de \le \frac{C}{\epsilon} |p - p'|, \quad (2.52) $$

where $C$ is a positive constant.

### 3. NUMERICAL ALGORITHM

In this section, we first consider various versions of the splitting scheme that will be implemented in practice. We then introduce the numerical errors that have to be taken into account in order to obtain the convergence of the scheme. The main error decomposition is given at the end of this section, see Proposition 3.1, together with our main convergence result, see Theorem 3.1. The precise study of each error is postponed to the next section.
### 3.1. Fully implementable schemes

The schemes we consider are all based on the approximation of the Wiener measure by a discrete and finite set of path. Indeed, on the grid $\pi$, the Brownian motion is approximated by a discrete time process denoted $\widehat{W}$. At each date $t_n$, $\mathcal{S}_n^d \subset \mathbb{R}^d$ denotes the support of the law of the random variable $\widehat{W}_{t_n}$. To define $\widehat{W}$, we first define $(\Delta \widehat{W}_n := \widehat{W}_{t_{n+1}} - \widehat{W}_{t_n})_{0 \le n \le N-1}$ which stands for discrete approximation of the Brownian increments $(W_{t_{n+1}} - W_{t_n})_{0 \le n \le N-1}$.

**Definition 3.1** (Brownian increments approximation). We use the cubature formula introduced in Section A.2 of [9] which requires only $2d$ points. Denoting $(\mathfrak{e}^\ell)_{1 \le \ell \le d}$ the canonical basis of $\mathbb{R}^d$, we set, for $1 \le i \le I = 2d$, $\mathbb{P}(\Delta \widehat{W}_n = \omega_h^i) = \frac{1}{2d}$ and $\omega_h^i = -\sqrt{d\mathfrak{h}\mathfrak{e}^\ell}$, if $i = 2\ell$ or $\omega_h^i = \sqrt{d\mathfrak{h}\mathfrak{e}^\ell}$ if $i = 2\ell - 1$.

By construction, we observe that

$$ \text{for } x \in \mathcal{S}_n^d, x + \Delta \widehat{W}_n \in \mathcal{S}_{n+1}^d. \qquad (3.1) $$

Moreover, the set $\mathcal{S}_n^d$ are discrete and finite.

For later use, we make the following observation.

**Lemma 3.1.** The number of points in $\mathcal{S}_n^d$, denoted $\sharp\mathcal{S}_n^d$, is $O(n^d)$.

*Proof.* By definition, $\mathcal{S}_n^d$ is a subgrid of $\sqrt{d\mathfrak{h}}\mathbb{Z}^d$. Without loss of generality and for ease of presentation, we assume for this proof that $\sqrt{d\mathfrak{h}} = 1$. For $z \in \mathbb{Z}^d$, we denote $|z|_1 = \sum_{\ell=1}^d |z^\ell|$. We now present two ways of describing the set $\mathcal{S}_n^d$, for $n \ge 0$ and $d \ge 1$.

(1) Assume that $n \ge 2$. We first observe that the points in $\mathcal{S}_n^d$ are the points that can be attained after $n$ steps following at each step any direction in $\{\omega_h^i | 1 \le i \le 2d\}$ according to Definition 3.1. To go from step $n-1$ to step $n$, there are two possibilities: the direction that is followed ends up in $\mathcal{S}_{n-2}^d$ or it allows to reach $\{z \in \mathbb{Z}^d | |z|_1 = n\}$ (the points that are exactly at distance $n$ from the origin). We thus have, for $n \ge 2$,

$$ \mathcal{S}_n^d = \mathcal{S}_{n-2}^d \cup \{z \in \mathbb{Z}^d | |z|_1 = n\}, \qquad (3.2) $$

and $\mathcal{S}_1^d = \{\omega_h^i | 1 \le i \le 2d\}$, $\mathcal{S}_0^d = \{0\}$. Thus,

$$ \mathcal{S}_{2m}^d = \bigcup_{k=0}^{m} \{z \in \mathbb{Z}^d \mid |z|_1 = 2k\} \quad \text{and} \quad \mathcal{S}_{2m+1}^d = \bigcup_{k=0}^{m} \{z \in \mathbb{Z}^d \mid |z|_1 = 2k+1\}. $$

In particular, for $d=1$, we have $\sharp\{z \in \mathbb{Z}^d | |z|_1 = n\} = 2$ for $n > 0$, which yields then that

$$ \sharp \mathcal{S}_n^1 = n + 1, \quad n \ge 0. \qquad (3.3) $$

(2) For $d \ge 2, n \ge 1$, the set $\mathcal{S}_n^d$ can be described by induction on the dimension. Indeed, we observe that

$$ \mathcal{S}_n^d = \bigcup_{k \in \{-n, \dots, n\}} \{(z, k) \in \mathbb{Z}^d | z \in \mathcal{S}_{n-k}^{d-1}\}. $$

We obtain then

$$ \sharp \mathcal{S}_n^d = \sharp \mathcal{S}_n^{d-1} + 2 \sum_{k=1}^{n} \sharp \mathcal{S}_{n-k}^{d-1}. \qquad (3.4) $$

An easy induction, starting with (3.3), yields then that $\sharp\mathcal{S}_n^d = O(n^d)$. Note that for $d=2$, one computes $\sharp\mathcal{S}_n^2 = (n+1)^2$, but one has $\sharp\mathcal{S}_1^3 = 6$.

□
### 3.1.1. Schemes for generic diffusion

In this section, we treat the general case where $P$ is solution to a SDE with Lipschitz coefficients as in (1.1). The discrete approximation on the grid $\pi$ is given by the classical Euler Scheme

$$ \hat{P}_0 = P_0 \quad \text{and} \quad \hat{P}_{t_{n+1}} = \hat{P}_{t_n} + b(\hat{P}_{t_n})\mathfrak{h} + \sigma(\hat{P}_{t_n})\Delta\hat{W}_n, \qquad (3.5) $$

where $(\Delta\hat{W}_n)$ is introduced in Definition 3.1. For later use, we also define, for $p \in \mathbb{R}^d$ and any $0 \le n \le N - 1$,

$$ \hat{P}_{t_{n+1}}^{t_n, p} := p + b(p)\mathfrak{h} + \sigma(p)\Delta\hat{W}_n, \qquad (3.6) $$

and, for $1 \le i \le 2d$,

$$ (\hat{P}_{t_{n+1}}^{t_n, p})^i := p + b(p)\mathfrak{h} + \sigma(p)\omega_{\mathfrak{h}}^i. \qquad (3.7) $$

Let us denote by $\mathcal{P}_n$ the discrete and finite support of $\hat{P}_{t_n}$ for $0 \le n \le N - 1$. We observe, due to the very definition of $(\hat{P}_{t_n})$, that

$$ \text{for } p \in \mathcal{P}_n, \quad \hat{P}_{t_{n+1}}^{t_n, p} \in \mathcal{P}_{n+1}. \qquad (3.8) $$

In this context, we first introduce a discrete version of the operator $\mathcal{T}$, recall Definition 2.4, that will compute an approximation to (2.15) written in forward form: We shall use the celebrated Sticky Particle Dynamics (SPD) [2] see also Section 1.1 of [11]. The SPD is particularly simple to implement in our case, since, due to the monotonicity assumption on $(\mu, \psi)$, there is no particle colliding!

For $M \ge 1$, let

$$ D_M := \{e = (e_1, \dots, e_m, \dots, e_M) \in \mathbb{R}^M \mid e_1 \le \dots \le e_m \le \dots \le e_M\} \qquad (3.9) $$

$$ \text{and } \mathcal{I}_M := \left\{ \theta \in \mathcal{J} \mid \theta(\cdot) := H * \left( \frac{1}{M} \sum_{m=1}^{M} \delta_{e_m} \right), e \in D_M \right\} \qquad (3.10) $$

where $H$ is the Heaviside function $x \mapsto \mathbf{1}_{\{x \ge 0\}}$, * the convolution operator and $\delta_e$ is the Dirac mass at $e \in \mathbb{R}$.

For later use, we introduce, for $M \ge 0$ and $n \in \{0, \dots, N\}$, the set of functions

$$ \mathcal{K}_n^M = \{ \psi : \mathcal{P}_n \times \mathbb{R} \to [0, 1] \mid \text{ for } p \in \mathcal{P}_n, \psi(p, \cdot) \in \mathcal{I}_M \}. \qquad (3.11) $$

The discrete version of $\mathcal{T}$, denoted $\mathfrak{T}$, acts on empirical CDF belonging to $\mathcal{I}_M$ and is given by

$$ (0, \infty) \times \mathbb{R}^d \times \mathcal{I}_M \ni (h, p, \theta) \mapsto \mathfrak{T}_h^M(p, \theta) := H * \left( \frac{1}{M} \sum_{m=1}^{M} \delta_{e_m(h)} \right) \in \mathcal{I}_M. \qquad (3.12) $$

Above $(e_m^p(h))_{1 \le m \le M}$ is a set (of positions) of particles computed as follows. Given the initial position $e(0) \in D_M$ (representing $\theta$) and velocities $(\bar{F}_m(p))_{1 \le m \le M}$ set to $\bar{F}_m(p) = - \int_{(m-1)/M}^{m/M} \mu(p, y) dy$, we consider $M$ particles $(e_m^p)_{1 \le m \le M}$, whose positions at time $t \in [0, h]$ are simply given by

$$ e_m^p(t) = e_m(0) + \bar{F}_m(p)t. \qquad (3.13) $$

We observe that $(e_m^p(t))_{1 \le m \le M} \in D_M$, for all $t \in [0, h]$, as $-\mu$ is non-decreasing.

We now present the approximation used for the approximation of $\mathcal{D}_\mathfrak{h}(\psi)$, recall Definition 2.5 and (2.18).

Let $\psi \in \mathcal{K}_{n+1}^M$ and fix $p \in \mathcal{P}_n$. We observe

$$ \bar{v}_n(p, e) := \mathbb{E} \left[ \psi \left( \hat{P}_{t_{n+1}, p}^{t_n, p}, e \right) \right] = \frac{1}{2d} \sum_{i=1}^{2d} \psi \left( \left( \hat{P}_{t_{n+1}, p}^{t_n, p} \right)^i, e \right) \qquad (3.14) $$
recall (3.7).

Since $\psi \in \mathcal{K}_{n+1}^M$, recall (3.11), we know that, for $1 \le i \le 2d$, there is some $e^i \in D_M$, depending implicitly on $\psi((\hat{P}_{t_{n+1}}^{t_n, \mathbf{P}})^i, \cdot)$, such that

$$ \psi\left(\left(\hat{P}_{t_{n+1}}^{t_n, \mathbf{P}}\right)^i, e\right) = H * \left(\frac{1}{M} \sum_{m=1}^{M} \delta_{e_m^i}\right)(e). $$

By linearity, (3.14) reads

$$ \bar{v}_n(\mathbf{p}, e) = H * \left( \frac{1}{2dM} \sum_{i=1}^{I} \sum_{m=1}^{M} \delta_{e_m^i} \right) (e). \quad (3.15) $$

The approximation of the diffusion operator is thus given by

$$ \mathcal{K}_{n+1}^M \ni \psi \mapsto \mathcal{D}_n^M(\psi) := \bar{v}_n \in \mathcal{K}_n^{2Md}. \quad (3.16) $$

Finally, the scheme will have essentially two versions:

**CASE 1:** We keep all the particles at each step. The overall procedure will then be the iteration of the two operators $\mathfrak{T}$ and $\mathfrak{D}$, see Definition 3.2 below.

**CASE 2:** We keep the number of particles constant between two steps. To reduce the number of particles after the use of $\mathfrak{T}$, we apply another operator $\mathfrak{R}$, namely, for $M \ge m \ge 1$,

$$ \mathcal{I}^M \ni \psi \mapsto \mathfrak{R}^{M,m}(\psi) \in \mathcal{I}^m. \quad (3.17) $$

Various implementations are possible, we refer to the numerical section for a precise description. The overall procedure for this case is given in Definition 3.3 below.

Let us now finally introduce the scheme formally.

**Definition 3.2 (Scheme in CASE 1).** Fix $M \ge 1$ and set for $0 \le n \le N$, $M_n := (2d)^{N-n}M$.

(1) At $n = N$: Set $e_N := (\Lambda, \dots, \Lambda) \in D_{M_N}$ whose empirical CDF is the terminal condition $\phi = \mathbf{1}_{\{e \ge \Lambda\}}$, recall (1.2). Setting simply $u_N^{N,M}(p, \cdot) = \phi$, for all $p \in \mathcal{P}_N$, we observe that $u_N^{N,M} \in \mathcal{K}_N^{M_N}$.

(2) For $n < N$: Given $u_{n+1}^{N,M} \in \mathcal{K}_{n+1}^{M_{n+1}}$, define

$$ \bar{u}_n^{N,M} = \mathfrak{D}_n^{M_{n+1}}(u_{n+1}^{N,M}) \in \mathcal{K}_n^{M_n}, \quad (3.18) $$

recall (3.16), and then $u_n^{N,M}$ by, for each $p \in \mathcal{P}_n$,

$$ u_n^{N,M}(p, \cdot) = \mathfrak{T}_h^{M_n}(p, \bar{u}_n^{N,M}(p, \cdot)), \quad (3.19) $$

recall (3.12).

The approximation of $\mathcal{V}(0, P_0, \cdot)$ is then given by $u_0^{N,M}(P_0, \cdot)$.

**Definition 3.3 (Scheme in CASE 2).** Fix $M \ge 1$.

(1) At $n = N$: Set $e_N := (\Lambda, \dots, \Lambda) \in D_M$ whose empirical CDF is the terminal condition $\phi = \mathbf{1}_{\{e \ge \Lambda\}}$, recall (1.2). Setting simply $v_N^{N,M}(p, \cdot) := \phi$, for all $p \in \mathcal{P}_N$, we observe that $v_N^{N,M} \in \mathcal{K}_N^M$.
(2) For $n < N$: Given $v_{n+1}^{N,M} \in \mathcal{K}_{n+1}^M$, define $\bar{v}_n^{N,M} = \mathfrak{D}_n^M(v_{n+1}^{N,M}) \in \mathcal{K}_n^{2dM}$, recall (3.16), and then for each $p \in \mathcal{P}_n$,

$$ \tilde{v}_n^{N,M}(p, \cdot) = \mathfrak{R}^{2dM,M}(\bar{v}_n^{N,M}(p, \cdot)) \qquad (3.20) $$

recall (3.17), and finally $v_n^{N,M}(p, \cdot)$ by,

$$ v_n^{N,M}(p, \cdot) = \mathfrak{F}_h^M(p, \tilde{v}_n^{N,M}(p, \cdot)), \qquad (3.21) $$

recall (3.12).

The approximation of $\mathcal{V}(0, P_0, \cdot)$ is then given by $v_0^{N,M}(P_0, \cdot)$.

**Remark 3.1.** Let us insist on the fact that, unless some recombination occurs at the tree level for $P$, the above schemes become rapidly untractable. The worst case complexity is easily obtained for the tree $\hat{P}$: at each step $n$, there is at most $(2d)^n$ branches, $0 \le n \le N-1$. Implementing these schemes, one faces two issues: the number of operations to perform at each step when applying the discrete transport operator, see (3.12), and the number of particles to keep to be able to apply the approximation of the diffusion operator (3.15). Indeed, if one is only interested in the approximation of $\mathcal{V}(0, P_0, \cdot)$ then the sum appearing in (3.15) is only computed at time 0 but again, in order to do that, one needs to keep track of all the particles until time 0.

In **CASE 1**, see Definition 3.2, we obtain by backward induction that the number of particles at each step for a given branch is $M_n = (2d)^{N-n}M$. At each step $n$, the discrete transport operator must be called $(2d)^n$ times (the number of branches at most), each call requiring to recompute the position of $M_{n+1}$ particles, see (3.13): this will cost $O((2d)^{N-1}M)$ for step $n < N$. And at step 0, one will evaluate (3.15) that will cost $O((2d)^N M)$. Overall, the complexity is $O((2d)^NMN)$.

In **CASE 2**, see Definition 3.3, the only improvement comes from the fact that the number of particles at each step for a given branch is always $M$. The cost of applying the discrete transport operator at step $n$ is now $(2d)^n M$, summing on the total number of steps we obtain a complexity of $O((2d)^N M)$. And at step 0, one will evaluate (3.15) that will cost $O(2dM)$. Overall, the complexity is reduced to $O((2d)^NM)$.

### 3.1.2. A simplified functional Brownian setting

We now consider a special case for the process $P$ for which the numerical implementation is less computationally demanding, see Remark 3.3 below.

**Assumption 3.1.** The process $P$ is given as a function of the Brownian motion, namely

$$ P_t = \mathfrak{P}(t, W_t) \text{ for } \mathfrak{P}: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d. $$

**Remark 3.2.** The main application of this case is when $P$ is the Brownian motion itself (also known as the Bachelier model in financial application) or a Geometric Brownian Motion (also known as Black-Scholes model in financial application).

Compared to the previous section, the only difference is the approximation of $P$. In this context, it is naturally given by

$$ \hat{P}_t := \mathfrak{P}(t, \widehat{W}_t), \quad t \in \pi. \qquad (3.22) $$

Note that we slightly abuse the notation here by keeping the same ones as the ones introduced in the previous section.

So, for $0 \le n \le N$,

$$ \mathcal{P}_n = \{p = \mathfrak{P}(t_n, w), w \in \mathfrak{S}_n\}, \qquad (3.23) $$
which is the (discrete) support of $\hat{P}_{t_n}$. We also define, for $p = \mathfrak{P}(t_n, w)$,

$$ \hat{P}_{t_{n+1}}^{t_{n,p}} = \mathfrak{P}(t_{n+1}, w + \Delta\hat{W}_n) \quad \text{and} \quad (\hat{P}_{t_{n+1}}^{t_{n,p}})^i = \mathfrak{P}(t_{n+1}, w + \omega_h^i), \qquad (3.24) $$

for $1 \le i \le 2d$, recall Definition 3.1.

Now that the approximation of $P$ has been clarified in this context, we can use the various schemes given in Definitions 3.2 and 3.3.

**Remark 3.3.** (i) The main difference between the generic diffusion case and the functional Brownian setting comes from the numerical complexity associated to their implementation. Indeed, it is the case that the growth size of $\mathfrak{S}_n^d$ with respect to $n$ is tamed (for low dimensional problem) because the “tree” associated to the brownian motion is naturally recombining, which is a priori not the case for generic diffusion. According to Lemma 3.1, $\#\mathfrak{S}_n^d = O(n^d)$. Following Remark 3.1, we can estimate the complexity of the algorithms in the functional Brownian setting.

In **CASE 1**, we have $M_{n+1} = O(M(N-n+1)^d)$ and $O(n^d)$ branches at each step $n$, we thus compute that the complexity is $O(M \sum_{n=0}^{N-1} (N-n)^d n^d) = O(MN^{2d+1})$ for large $N$. In **CASE 2**, we have $M_{n+1} = M$ and still $O(n^d)$ branches at each step $n$, we thus compute that the complexity is $O(M \sum_{n=0}^{N-1} n^d) = O(MN^{d+1})$ for large $N$.

Thus, the methods we study here are impacted by the “curse of dimensionality”. However, we should also note that for dimension around 5 and in the case of the brownian setting, the problems are tractable numerically, in particular when working with **CASE 2**, see Section 5 for illustration.

(ii) For implementation purposes, we just use the Brownian setting. A way to control the growth for generic diffusion is to introduce some interpolation on a grid see e.g. [5] or to use recombination techniques developed for cubature methods, see e.g. [12].

## 3.2. Error decomposition and main result

We first present the various sources of errors introduced by the schemes above. In the sequel, we shall mainly conduct our analysis for the scheme given in **CASE 1**.

The error we seek to control is, for $\gamma_n = u_n^{N,M}$ or $\gamma_n = v_n^{N,M}$, recall Definitions 3.2 and 3.3,

$$ \text{err}(N, M) := \int_{\mathbb{R}} |\gamma_0(P_0, e) - \mathcal{V}(0, P_0, e)| \, de. \qquad (3.25) $$

And we first observe

$$ \text{err} \le \mathcal{E}_0(P_0) + \mathcal{E}_0^r, \qquad (3.26) $$

with for $0 \le n \le N$ and $p \in \mathcal{P}_n$,

$$ \mathcal{E}_n^r(p) := \int_{\mathbb{R}} |\mathcal{V}^\epsilon(t_n, p, e) - \mathcal{V}(t_n, p, e)| \, de, \qquad (3.27) $$

$$ \mathcal{E}_n(p) := \int_{\mathbb{R}} |\gamma_n(p, e) - \mathcal{V}^\epsilon(t_n, p, e)| \, de. \qquad (3.28) $$

In **CASE 1**, we define, for $p \in \mathcal{P}_n$,

$$ \mathcal{E}_n^T(p) := \int_{\mathbb{R}} |\tilde{\mathfrak{F}}_h^{M_n}(p, \bar{u}_n^{N,M}(p, \cdot))(e) - \tilde{T}_h(p, \bar{u}_n^{N,M}(p, \cdot))(e)| \, de, \qquad (3.29) $$

$$ \bar{\mathcal{E}}_n(p) := \int_{\mathbb{R}} |\tilde{T}_h(p, \bar{u}_n^{N,M}(p, \cdot))(e) - \tilde{T}_h(p, \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, \hat{P}_{t_{n+1},p}^{t_{n,p}})])(e)| \, de, \qquad (3.30) $$
recall Definition 3.2 and Remark 2.1.

Moreover, for $p \in \mathbb{R}^d$, we set

$$ \mathcal{E}_n^D(p) := \int_{\mathbb{R}} |\tilde{T}_h(p, \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, \hat{P}_{t_{n+1},p}, \cdot)])](e) - \tilde{T}_h(p, \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, P_{t_{n+1},p}^{t_{n,p}}, \cdot)])](e) |\mathrm{d}e, \quad (3.31) $$

$$ \mathcal{E}_n^S(p) := \int_{\mathbb{R}} |\tilde{T}_h(p, \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, P_{t_{n+1},p}^{t_{n,p}}, \cdot)])](e) - \mathcal{V}^\epsilon(t_n, p, e)|\mathrm{d}e. \quad (3.32) $$

We have the following.

**Lemma 3.2.** Under CASE 1, for $0 \le n \le N-2$,

$$ \mathcal{E}_n(p) \le \mathcal{E}_n^T(p) + \bar{\mathcal{E}}_n(p) + \mathcal{E}_n^D(p) + \mathcal{E}_n^S(p) \quad (3.33) $$

and

$$ \mathcal{E}_{N-1}(p) \le \mathcal{E}_{N-1}^{\mathrm{T}}(p) + \mathcal{E}_{N-1}^{\mathrm{r}}(p) + \int_{\mathbb{R}} |\mathcal{S}_{\mathfrak{h}}(\phi(\cdot))(p, e) - \Theta_{\mathfrak{h}}(\phi(\cdot))(p, e)| \mathrm{d}e \quad (3.34) $$

observing that

$$ \mathcal{E}_{N-1}^{\mathrm{T}}(p) = \int_{\mathbb{R}} |\mathcal{T}_{\mathfrak{h}}^{M}(p, \phi(\cdot))(e) - \mathcal{T}_{\mathfrak{h}}(p, \phi(\cdot))(e)| \mathrm{d}e. \quad (3.35) $$

*Proof.* The first inequality is straightforward. For the second, we observe that $u_{N,M} = \phi$ and does not depend on $p$ variable so that $\bar{u}_{N-1} = \phi$.

$$ \begin{aligned} \mathcal{E}_{N-1}(p) &:= \int_{\mathbb{R}} |u_{N-1}^{N,M}(p, e) - \mathcal{V}^{\epsilon}(t_{N-1}, p, e)| \mathrm{d}e \\ &= \int_{\mathbb{R}} |\mathcal{T}_{\mathfrak{h}}^{M}(p, \phi(\cdot))(e) - \mathcal{V}^{\epsilon}(t_{N-1}, p, e)| \mathrm{d}e \\ &\le \int_{\mathbb{R}} |\mathcal{T}_{\mathfrak{h}}^{M}(p, \phi(\cdot))(e) - \mathcal{T}_{\mathfrak{h}}(p, \phi(\cdot))(e)| \mathrm{d}e + \int_{\mathbb{R}} |\mathcal{T}_{\mathfrak{h}}(p, \phi(\cdot))(e) - \mathcal{V}(t_{N-1}, p, e)| \mathrm{d}e \\ &\quad + \int_{\mathbb{R}} |\mathcal{V}(t_{N-1}, p, e) - \mathcal{V}^{\epsilon}(t_{N-1}, p, e)| \mathrm{d}e. \end{aligned} $$

We then observe that $\mathcal{T}_{\mathfrak{h}}(p, \phi(\cdot))(e) = \mathcal{S}_{\mathfrak{h}}(\phi(\cdot))(p, e)$ as again $\phi$ does not depend on $p$. $\square$

We now comment quickly the various errors appearing above: The term $\mathcal{E}_n^T(p)$ is the local error due to the approximation of the transport operator by the SPD. The term $\bar{\mathcal{E}}_n(p)$ will allow to propagate the local errors thanks to the $L^1$ stability of $\mathcal{T}_h(\cdot)$. The term $\mathcal{E}^D$ is linked to the approximation of the diffusion by the tree and $\mathcal{E}^S$ is the splitting error applied to the proxy $\mathcal{V}^\epsilon$.

**Proposition 3.1** (Stability). Under CASE 1, the following hold

$$ \mathcal{E}_0(P_0) \le \sum_{n=0}^{N-2} \mathbb{E} \left[ \mathcal{E}_n^T(\hat{P}_{t_n}) + \mathcal{E}_n^D(\hat{P}_{t_n}) + \mathcal{E}_n^S(\hat{P}_{t_n}) \right] + \mathbb{E} \left[ \mathcal{E}_{N-1}(\hat{P}_{t_{N-1}}) \right]. \quad (3.36) $$

*Proof.* Since $\tilde{T}_h$ is $L^1$-stable, (3.30) yields for $p \in \mathcal{P}_n$,

$$ \bar{\mathcal{E}}_n(p) \le \int_{\mathbb{R}} |\bar{u}^{N,M}(t_n, p, e) - \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, \hat{P}_{t_{n+1},p}^{t_{n,p}}, e)]| \mathrm{d}e. $$
Under CASE 1, recall (3.18), we then get

$$ \bar{\mathcal{E}}_n(p) \le \mathbb{E} \left[ \left| \int_{\mathbb{R}} u^{N,M} \left( t_{n+1}, \hat{P}_{t_{n+1},p}^t, e \right) - \mathcal{V}^\epsilon \left( t_{n+1}, \hat{P}_{t_{n+1},p}^t, e \right) \right| \mathrm{d}e \right]. \quad (3.37) $$

Thus,

$$ \mathbb{E}[\bar{\mathcal{E}}_n(\hat{P}_{t_n})] \le \mathbb{E}[\mathcal{E}_{n+1}(\hat{P}_{t_{n+1}})]. \quad (3.38) $$

We thus deduce from Lemma 3.2

$$ \mathbb{E}[\mathcal{E}_n(\hat{P}_{t_n})] \le \mathbb{E}[\mathcal{E}_n^T(\hat{P}_{t_n}) + \mathcal{E}_n^D(\hat{P}_{t_n}) + \mathcal{E}_n^S(\hat{P}_{t_n})] + \mathbb{E}[\mathcal{E}_{n+1}(\hat{P}_{t_{n+1}})]. \quad (3.39) $$

The proof is concluded by induction observing $\mathbb{E}[\mathcal{E}_N(\hat{P}_{t_N})] = 0$. $\square$

The next section is dedicated to the study of the various error appearing above. We prove in Section 4.5 below the main theoretical result of this work:

**Theorem 3.1.** Let Assumption 2.1 hold. Then, under CASE 1,

$$ \mathrm{err}(N, M) \le C \left( \frac{1}{M} + \frac{\sqrt{\mathfrak{h}}}{\epsilon^2} + \epsilon \right), \quad (3.40) $$

recall (3.25). Moreover, setting $\epsilon = \mathfrak{h}^{\frac{1}{6}}$ and $M = \frac{1}{\epsilon}$, we have

$$ \int_{\mathbb{R}} |\mathcal{V}(0, P_0, e) - u_0^{N,M}(P_0, e)| \mathrm{d}e \le C \mathfrak{h}^{\frac{1}{6}}. \quad (3.41) $$

# 4. STUDY OF THE ERRORS

## 4.1. Approximation of transport operator

We discuss here the error introduced by the use of the SPD approximation in our framework. The numerical analysis of this method is now well known see e.g. [11].

**Lemma 4.1.** Under Assumption 2.1 and in CASE 1, the following holds

$$ \mathcal{E}_n^T(p) \le C \frac{\mathfrak{h}}{M_n}, \quad p \in \mathcal{P}_n. \quad (4.1) $$

*Proof.* By Definition 3.2, for $p \in \mathcal{P}_n$ and $n \le N - 2$, we have

$$ \mathcal{E}_n^T(p) = \int_{\mathbb{R}} |\mathfrak{T}_{\mathfrak{h}}^{M_n}(\bar{u}_n^{N,M}(p, \cdot), p) - \tilde{\mathfrak{T}}_{\mathfrak{h}}(\bar{u}_n^{N,M}(p, \cdot), p)| \mathrm{d}e \quad (4.2) $$

with $\bar{u}_n^{N,M}(p, \cdot) \in \mathcal{I}_{M_n}$. For fixed $p$, the property of $\mathfrak{T}_{\mathfrak{h}}^M$ are discussed in Section 1.1 of [11]. We use here the estimate given in Theorem 3.1 of [11], observing that, in our context, there is no error due to the discretization of the initial condition. Indeed, $\bar{u}_n^{N,M}(p, \cdot)$ appears already as the CDF of an empirical distribution. Thus, straightforwardly, we get

$$ \mathcal{E}_n^T(p) \le C \frac{\mathfrak{h}}{M_n}, \quad (4.3) $$

where $C$ does not depend on $p$ but depends on the Lipschitz constant of $\partial_y\mu$, recall Assumption 2.1.

For the step $n = N - 1$, we have,

$$ \mathcal{E}_{N-1}^{T}(p) = \int_{\mathbb{R}} |\mathfrak{T}_{\mathfrak{h}}^{M}(p, \phi(\cdot))(e) - \mathfrak{T}_{\mathfrak{h}}(p, \phi(\cdot)(e))| \mathrm{d}e \quad (4.4) $$

recall (3.35). It turns out that for our application $\phi$ is trivially the CDF of an empirical distribution, recall (1.2).

Thus, we indeed have, $\mathcal{E}_{N-1}^T(p) \le C \frac{\mathfrak{h}}{M_{N-1}}$. $\square$
## 4.2. Regularization error

**Proposition 4.1.** Under Assumption 2.1, the following holds, for $n \le N$,

$$ \mathcal{E}_n^r(p) \le C\epsilon, \quad p \in \mathbb{R}^d. $$

*Proof.* We first observe that $\mathcal{E}_N^r(p) = 0$ as $\phi$ does not depend on $p$. For $n < N$, we have that

$$ \mathcal{E}_n^r(p) = \int_{\mathbb{R}} |\mathcal{V}^\epsilon(t_n, p, e) - \mathcal{V}(t_n, p, e)| \, de \quad (4.5) $$

$$ = \int_{\mathbb{R}} \left| \int_{\mathbb{R}^d} \{\mathcal{V}(t_n, p+q, e) - \mathcal{V}(t_n, p, e)\} \varphi_\epsilon(q) \, dq \right| \, de \quad (4.6) $$

$$ \le \int_{\mathbb{R}^d} \int_{\mathbb{R}} |\mathcal{V}(t_n, p+q, e) - \mathcal{V}(t_n, p, e)| \, de \varphi_\epsilon(q) \, dq. \quad (4.7) $$

Thus, using (2.24), we compute

$$ \mathcal{E}_0^r(p) \le C \int_{\mathbb{R}^d} |q| \varphi_\epsilon(q) \, dq \le C\epsilon, \quad (4.8) $$

which concludes the proof. $\square$

## 4.3. Splitting error

Recall that the splitting error is given by

$$ \mathcal{E}_n^S(p) = \int_{\mathbb{R}} \left| \tilde{T}_h \left( p, \mathbb{E} \left[ \mathcal{V}^\epsilon \left( t_{n+1}, P_{t_{n+1}, p}^{t_n, p}, \cdot \right) \right] \right) (e) - \mathcal{V}^\epsilon(t_n, p, e) \right| \, de. \quad (4.9) $$

for $0 \le n \le N-1$, $p \in \mathbb{R}^d$.

In [6], the error due the theoretical splitting has already been studied and the results obtained there can be used in our setting. However, we should point out that here $\mathcal{V}^\epsilon$ appears instead of $\mathcal{V}$. We thus first observe the following.

**Lemma 4.2.** Under Assumption 2.1, the following holds, for $p \in \mathbb{R}^d$,

$$ \mathcal{E}_n^S(p) \le C(1 + |p|^2)h^{\frac{3}{2}} + \mathfrak{C}_n(p) \quad (4.10) $$

with

$$ \mathfrak{C}_n(p) := \int_{\mathbb{R}} |\Theta_h(\mathcal{V}^\epsilon(t_{n+1}, \cdot))(p, e) - \mathcal{V}^\epsilon(t_n, p, e)| \, de. \quad (4.11) $$

*Proof.* We compute

$$ \mathcal{E}_n^S(p) \le \int_{\mathbb{R}} \left| \tilde{T}_h \left( p, \mathbb{E} \left[ \mathcal{V}^\epsilon \left( t_{n+1}, P_{t_{n+1}, p}^{t_n, p}, \cdot \right) \right] \right) (e) - \Theta_h(\mathcal{V}^\epsilon(t_{n+1}, \cdot))(p, e) \right| \, de \\ + \int_{\mathbb{R}} |\Theta_h(\mathcal{V}^\epsilon(t_{n+1}, \cdot))(p, e) - \mathcal{V}^\epsilon(t_n, p, e)| \, de. $$

The first term in the RHS above is then controlled by invoking Proposition 2.1. $\square$

It remains to study $\mathfrak{C}_n$, recall (4.11). The upper bound for this term is obtained in Proposition 4.2 below, which allows to conclude for the splitting error in Corollary 4.1. We first need the following result.
**Lemma 4.3.** Let Assumption 2.1 hold. Then, for $n \le N - 2$,

$$ \sup_{(t,p,e) \in [t_n, t_{n+1}] \times \mathbb{R}^d \times \mathbb{R}} |\Theta_{t_{n+1}-t}(\mathcal{V}^\epsilon(t_{n+1}, \cdot))(p, e) - \mathcal{V}^\epsilon(t, p, e)| \le C\epsilon, \quad (4.12) $$

where C is a positive constant.

*Proof.* Without loss of generality, we do the proof for $n = 0$, working on $[0, t_1 = 0]$. We denote $w(s, \cdot) = \Theta_{t_1-s}(\mathcal{V}^\epsilon(t_1, \cdot))$ for $s \in [0, t_1]$. Recall that the following FBSDE is well posed: for $s \in [0, t_1]$,

$$ \left\{ \begin{array}{l} P_s = p + \int_0^s b(P_t) \, dt + \bar{\sigma} W_s, \\ E_s = e + \int_0^s \mu(P_t, Y_t) \, dt, \\ Y_s = \mathcal{V}^\epsilon(t_1, P_{t_1}, E_{t_1}) - \int_s^{t_1} Z_t \, dW_t, \end{array} \right. \quad (4.13) $$

with $Y_s = w(s, P_s, E_s)$ for $s \in [0, t_1]$. Without loss of generality, we will prove (4.12) at the point $(0, p, e)$ (and for $n = 0$).

Let $V_t = Y_t - \mathcal{V}^\epsilon(t, P_t, E_t)$. Applying Ito's formula and using the martingale property of Y we get,

$$ dV_t = d\mathcal{M}_t - \mu(P_t, Y_t)\partial_e\mathcal{V}^\epsilon(t, P_t, E_t)\,dt \quad (4.14) $$

$$ - \left\{ \partial_t \mathcal{V}^\epsilon(t, P_t, E_t) + b(P_t) \partial_p \mathcal{V}^\epsilon(t, P_t, E_t) + \frac{1}{2} \Delta_p \mathcal{V}^\epsilon(t, P_t, E_t) \right\} dt \quad (4.15) $$

where $\mathcal{M}$ is a square integrable martingale with $\mathcal{M}_0 = 0$.

Using the PDE (2.40) satisfied by $\mathcal{V}^\epsilon$ we get

$$ dV_t = d\mathcal{M}_t - (\{\mu(P_t, Y_t) - \mu(P_t, \mathcal{V}^\epsilon(t, P_t, E_t))\}\partial_e\mathcal{V}^\epsilon(t, P_t, E_t) + \theta^\epsilon(t, P_t, E_t))\,dt. \quad (4.16) $$

Observe that

$$ \mu(P_t, Y_t) - \mu(P_t, \mathcal{V}^\epsilon(t, P_t, E_t)) = c_t V_t \quad \text{with} \quad c_t := \int_0^1 \partial_y \mu(P_t, \mathcal{V}^\epsilon(t, P_t, E_t) + \lambda V_t) \, d\lambda \quad (4.17) $$

and from the property of $\mu$, recall (2.4), we have

$$ c_t \le -\ell_1 < 0. \quad (4.18) $$

We get, for $s \in [0, t_1]$

$$ V_s - V_0 = - \int_0^s (c_t V_t \partial_e \mathcal{V}^\epsilon(t, P_t, E_t) + \theta^\epsilon(t, P_t, E_t)) \, dt + \mathcal{M}_s. \quad (4.19) $$

We set, for $0 \le t \le t_1$, $\mathcal{E}_t = e^{\int_0^t c_s \partial_e \mathcal{V}^\epsilon(s, P_s, E_s) \, ds}$ and, we have

$$ 0 \le \mathcal{E}_t \le 1, \quad \text{for all } 0 \le t \le t_1, \quad (4.20) $$

since $\partial_e \mathcal{V}^\epsilon$ is non negative, recall also (4.18). We then compute

$$ \mathcal{E}_s V_s - V_0 = - \int_0^s \theta^\epsilon(t, P_t, E_t) \mathcal{E}_t \, dt + \mathcal{N}_s \quad (4.21) $$
with $\mathcal{N}$ a square integrable martingale such that $\mathcal{N}_0 = 0$. In particular, recalling (2.41), we get

$$V_0 \le |\mathbb{E}[\mathcal{E}_s V_s]| + C\mathfrak{h}\epsilon + C\epsilon\mathbb{E}\left[\int_0^s \partial_e \mathcal{V}^\epsilon(t, P_t, E_t)\mathcal{E}_t dt\right] \quad (4.22)$$

where we use the fact that $|\partial_p \mathcal{V}|$ is bounded and $\partial_e \mathcal{V}^\epsilon(t, P_t, E_t) \ge 0$. The above inequality reads also

$$|V_0| \le |\mathbb{E}[\mathcal{E}_s V_s]| + C\mathfrak{h}\epsilon - C\epsilon\mathbb{E}\left[\int_0^s \frac{c_t}{|c_t|} \partial_e \mathcal{V}^\epsilon(t, P_t, E_t) \mathcal{E}_t dt\right]. \quad (4.23)$$

Observing that

$$\dot{\mathcal{E}}_t = c_t \partial_e \mathcal{V}^\epsilon(t, P_t, E_t) \mathcal{E}_t,$$

since $|c_t| \ge \ell_1 > 0$, we obtain

$$|V_0| \le |\mathbb{E}[\mathcal{E}_s V_s]| + C(1 + \mathfrak{h})\epsilon \quad (4.24)$$

and since $V_{t_1} = 0$

$$|w(0, p, e) - \mathcal{V}^\epsilon(0, p, e)| = |V_0| \le C\epsilon,$$

which concludes the proof. $\square$

**Proposition 4.2.** Let Assumption 2.1 hold. Then, for $n \le N - 2$, $p \in \mathbb{R}^d$,

$$\mathfrak{E}_n(p) = \int_{\mathbb{R}} |\Theta_{\mathfrak{h}}(\mathcal{V}^\epsilon(t_{n+1}, \cdot))(p, e) - \mathcal{V}^\epsilon(t_n, p, e)| \, de \le C\mathfrak{h}\epsilon. \quad (4.25)$$

*Proof.* Without loss of generality, we prove the statement in the same setting as the one used for the proof of Lemma 4.3. We just stress here the dependence upon the initial condition: $V_t^e = Y_t^e - \mathcal{V}^\epsilon(t, P_t, E_t^e)$ since $E_0 = e$. We consider the tangent process $\partial_e E$ given by

$$\partial_e E_t = 1 + \int_0^t \partial_y \mu(P_s, Y_s) \partial_e w(s, P_s, E_s) \partial_e E_s \, ds, \quad (4.26)$$

$$= e^{\int_0^t \partial_y \mu(P_s, Y_s) \partial_e w(s, P_s, E_s) \, ds}. \quad (4.27)$$

And we observe that $0 \le \partial_e E_t \le 1$, for all $0 \le t \le t_1$, since $\partial_e w \ge 0$ and $\mu$ is decreasing in $y$, recall (2.4).

In order to bound the error $\int_{\mathbb{R}} |V_0^e| \, de$, we will study the dynamics of $t \mapsto \int_{\mathbb{R}} |\mathbb{E}[V_t^e \partial_e E_t^e]| \, de$. Recalling (4.19), we compute

$$d(V_t^e \partial_e E_t) = d\mathcal{N}_t + V_t^e \partial_y \mu(P_t, Y_t) \partial_e w(t, P_t, E_t) \partial_e E_t \, dt \quad (4.28)$$

$$-(c_t V_t^e \partial_e \mathcal{V}^\epsilon(t, P_t, E_t) \partial_e E_t + \theta^\epsilon(t, P_t, E_t^e) \partial_e E_t) \, dt, \quad (4.29)$$

where $\mathcal{N}$ is a square integrable martingale satisfying $\mathcal{N}_0 = 0$.

$$\mathbb{E}[V_{\mathfrak{h}}^e \partial_e E_{\mathfrak{h}} - V_0^e] = \mathbb{E}\left[\int_0^{\mathfrak{h}} V_t^e (\partial_y \mu(P_t, Y_t) \partial_e [w(t, P_t, E_t^e)] - c_t V_t^e \partial_e [\mathcal{V}^\epsilon(t, P_t, E_t^e)]) dt\right] \quad (4.30)$$

$$+ \mathbb{E}\left[\int_0^{\mathfrak{h}} \theta^\epsilon(t, P_t, E_t) \partial_e E_t dt\right]. \quad (4.31)$$

Since $V_{\mathfrak{h}} = 0$ and $\partial_e w \ge 0$, $\partial_e \mathcal{V}^\epsilon \ge 0$, $\partial_e E_t \ge 0$, we obtain

$$\int_{\mathbb{R}} |V_0^e| \, de \le \int_0^{\mathfrak{h}} \mathbb{E}\left[\int_{\mathbb{R}} |V_t^e \partial_y \mu(P_t, Y_t)| \partial_e [w(t, P_t, E_t^e)] \, de \right] \, dt \quad (4.32)$$
$$ + \int_0^\beta \mathbb{E} \left[ \int_{\mathbb{R}} |c_t V_t^e| \partial_e [\mathcal{V}^\epsilon(t, P_t, E_t^e)] \mathrm{d}e \right] \mathrm{d}t \\ + \int_0^\beta \mathbb{E} \left[ \int_{\mathbb{R}} |\theta^\epsilon(t, P_t, E_t)| \partial_e E_t \mathrm{d}e \right] \mathrm{d}t. $$

We obtain, using Lemma 4.3 and the bound on $\|\partial_y\mu\|_\infty$,

$$ \int_{\mathbb{R}} |V_t^e \partial_y \mu(P_t, Y_t)| \partial_e[w(t, P_t, E_t^e)] \mathrm{d}e \le C\epsilon \int_{\mathbb{R}} \partial_e[w(t, P_t, E_t^e)] \mathrm{d}e, \quad (4.33) \\ \le C\epsilon, $$

where we used the fact that $w$ is bounded in the last inequality. Similar arguments, since $c$ is bounded, lead to

$$ \int_{\mathbb{R}} |c_t V_t^e| \partial_e [\mathcal{V}^\epsilon(t, P_t, E_t^e)] \mathrm{d}e \le C\epsilon. \quad (4.34) $$

From (2.41), we know that

$$ |\theta^\epsilon(t, P_t, E_t)| \le C\epsilon \left( \int_{\mathbb{R}^d} |\partial_p \mathcal{V}(t, P_t + q, E_t^e)| \varphi_\epsilon(q) \mathrm{d}q + \partial_e \mathcal{V}^\epsilon(t, P_t, E_t) \right). \quad (4.35) $$

We compute, using the change of variable $e \mapsto \zeta = E_t^e$

$$ \int_{\mathbb{R}} \int_{\mathbb{R}^d} |\partial_p \mathcal{V}(t, P_t + q, E_t^e)| \varphi_\epsilon(q) \mathrm{d}q \partial_e E_t^e \mathrm{d}e = \int_{\mathbb{R}} \int_{\mathbb{R}^d} |\partial_p \mathcal{V}(t, P_t + q, \zeta)| \varphi_\epsilon(q) \mathrm{d}q \mathrm{d}\zeta, \quad (4.36) $$

$$ \le C, \quad (4.37) $$

where we used the $L^1$-bound on $\partial_p V$ given in (2.23). We then compute

$$ \mathbb{E} \left[ \left| \int_{\mathbb{R}} |\theta^\epsilon(t, P_t, E_t)| \partial_e E_t \mathrm{d}e \right| \right] \le C\epsilon \left( 1 + \mathbb{E} \left[ \left| \int_{\mathbb{R}} \partial_e \mathcal{V}^\epsilon(t, P_t, E_t) \partial_e E_t \mathrm{d}e \right| \right] \right), \quad (4.38) $$

$$ \le C\epsilon. \quad (4.39) $$

The proof is concluded by combining the previous inequality, the estimates given in (4.33) and (4.34) with (4.32). $\square$

Combining Lemma 4.2 with the result Proposition 4.2, we obtain straightforwardly the following.

**Corollary 4.1.** Let Assumption 2.1 hold. Then, for $n \le N - 2$, $p \in \mathbb{R}^d$,

$$ \mathcal{E}_n^S(p) \le C(1 + |p|^2)\mathfrak{h}^{\frac{3}{2}} + C\epsilon\mathfrak{h}. \quad (4.40) $$

### 4.4. Diffusion error

We now study the term given in (3.31) and we straightforwardly observe

$$ \mathcal{E}_n^D(p) = \int_{\mathbb{R}} |\tilde{T}_h(p, \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, \hat{P}_{t_{n+1},p}, \cdot)])(e) - \tilde{T}_h(p, \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, P_{t_{n+1},p}, \cdot)])(e)| \mathrm{d}e, \quad (4.41) $$

$$ \le \int_{\mathbb{R}} |\mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, \hat{P}_{t_{n+1},p}, e)] - \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, P_{t_{n+1},p}, e)]| \mathrm{d}e, \quad (4.42) $$

where we use the $L^1$-stability of $\tilde{T}_h$.

This motivates the introduction of the following auxiliary result.
**Proposition 4.3.** Let Assumption 2.1 hold. For $0 \le n \le N-1$, let $w_n$ be the solution on $[t_n, t_{n+1}] \times \mathbb{R}^d$ to the following PDE:

$$ \partial_t w + b(p)\partial_p w + \frac{1}{2}\bar{\sigma}^2 \Delta_{pp} w = 0 \text{ and } w(t_{n+1}, \cdot) := \mathcal{V}^\epsilon(t_{n+1}, \cdot, e), \quad (4.43) $$

for a fixed $e \in \mathbb{R}$. Then, under Assumption 2.1, the followings hold, for $q, q' \in \mathbb{R}^d$ and $t \in [0, T)$,

$$ \int_{\mathbb{R}} |\partial_{p^\ell} w_n(t, q, e) - \partial_{p^\ell} w_n(t, q', e)| \, de \le \frac{C}{\epsilon} |q - q'| \quad (4.44) $$

and

$$ \int_{\mathbb{R}} |\partial_{p^\ell p^{\ell'}}^2 w_n(t, q, e) - \partial_{p^\ell p^{\ell'}}^2 w_n(t, q', e)| \, de \le \frac{C}{\epsilon^2} |q - q'|, \quad (4.45) $$

for all $\ell, \ell' \in \{1, \dots, d\}$. Moreover,

$$ \int |\partial_p w_n(t, p, e)| \, de \le C. \quad (4.46) $$

*Proof.* For ease of presentation, the proof is done in the one-dimensional case.

By Theorem 2.3.5 in Zhang [13], we have the following expressions for the first and second order derivative with respect to vector $p$, for $(t, p, e) \in [t_n, t_{n+1}] \times \mathbb{R} \times \mathbb{R}$,

$$ \partial_p w_n(t, p, e) = \mathbb{E} \left[ \partial_p \mathcal{V}^\epsilon \left( t_{n+1}, P_{t_{n+1}}^{t,p}, e \right) \partial_p P_{t_{n+1}}^{t,p} \right], \quad (4.47) $$

$$ \hat{\partial}_{pp}^2 w_n(t, p, e) = \mathbb{E} \left[ \hat{\partial}_{pp}^2 \mathcal{V}^\epsilon \left( t_{n+1}, P_{t_{n+1}}^{t,p}, e \right) \left( \hat{\partial}_p P_{t_{n+1}}^{t,p} \right)^2 + \hat{\partial}_p \mathcal{V}^\epsilon \left( t_{n+1}, P_{t_{n+1}}^{t,p}, e \right) \hat{\partial}_{pp}^2 P_{t_{n+1}}^{t,p} \right], \quad (4.48) $$

where, for $t \le s \le t_{n+1}$,

$$ P_s^{t,p} = p + \int_t^s b(P_r^{t,p}) \, dr + \bar{\sigma}(W_s - W_t), \quad (4.49) $$

$$ \hat{\partial}_p P_s^{t,p} = 1 + \int_t^s b'(P_r^{t,p}) \hat{\partial}_p P_r^{t,p} \, dr, \quad (4.50) $$

$$ \hat{\partial}_{pp}^2 P_s^{t,p} = \int_t^s \left( b''(P_r) |\hat{\partial}_p P_r^{t,p}|^2 + b'(P_r) \hat{\partial}_{pp}^2 P_r^{t,p} \right) dr. \quad (4.51) $$

In this setting, classical arguments lead to,

$$ \mathbb{E} \left[ \sup_{s \in [t, t_{n+1}]} |\hat{\partial}_p P_s^{t,p}|^\kappa + \sup_{s \in [t, t_{n+1}]} |\hat{\partial}_{pp}^2 P_s^{t,p}|^\kappa \right] \le C_\kappa, \quad (4.52) $$

and for $q, q' \in \mathbb{R}$, $s \in [t, T]$, $\kappa \ge 1$,

$$ \mathbb{E} \left[ \left| P_s^{t,q} - P_s^{t,q'} \right|^\kappa + \left| \hat{\partial}_p P_s^{t,q} - \hat{\partial}_p P_s^{t,q'} \right|^\kappa + \left| \hat{\partial}_{pp}^2 P_s^{t,q} - \hat{\partial}_{pp}^2 P_s^{t,q'} \right|^\kappa \right] \le C |q - q'|^\kappa. \quad (4.53) $$

We now study the second order derivative given in (4.48).

We compute, for $(t, q, q') \in [t_n, t_{n+1}] \times \mathbb{R} \times \mathbb{R}$,

$$ \int_{\mathbb{R}} |\hat{\partial}_{pp}^2 w_n(t, q, e) - \hat{\partial}_{pp}^2 w_n(t, q', e)| \, de \quad (4.54) $$
$$ \leq \mathbb{E} \left[ \left| \int_{\mathbb{R}} \left| \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q'}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| de \right] := B_1 \quad (4.55) $$

$$ + \mathbb{E} \left[ \left| \int_{\mathbb{R}} \left| \partial_p \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \partial_{pp}^2 P_{t_{n+1}}^{t,q} - \partial_p \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q'}, e) \partial_{pp}^2 P_{t_{n+1}}^{t,q'} \right| de \right] := B_2. \quad (4.56) $$

For the first term $B_1$,

$$ B_1 \leq \mathbb{E} \left[ \left| \int_{\mathbb{R}} \left| \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| de \right] + \mathbb{E} \left[ \left| \int_{\mathbb{R}} \left| \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 - \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q'}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| de \right]. $$

We then compute

$$ \mathbb{E} \left[ \left| \int_{\mathbb{R}} \left| \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| de \right], \quad (4.57) $$

$$ \leq \mathbb{E} \left[ \left| \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| \int_{\mathbb{R}} \left| \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \right| de \right], \quad (4.58) $$

$$ \leq \frac{C}{\epsilon} \mathbb{E} \left[ \left| \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| \right], \quad (4.59) $$

where we used for the last inequality Lemma 2.1. Observing that

$$ \left| \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| \leq \left| \partial_p P_{t_{n+1}}^{t,q} + \partial_p P_{t_{n+1}}^{t,q'} \right| \cdot \left| \partial_p P_{t_{n+1}}^{t,q} - \partial_p P_{t_{n+1}}^{t,q'} \right|, $$

we combine Cauchy-Schwarz inequality with (4.53) to obtain

$$ \mathbb{E} \left[ \left| \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| \right] \leq C |q - q'|. \quad (4.60) $$

Combining the previous inequality with (4.59), we get

$$ \mathbb{E} \left[ \left| \int_{\mathbb{R}} \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q} \right)^2 - \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| de \right] \leq \frac{C}{\epsilon} |q - q'|. \quad (4.61) $$

We also compute

$$ \begin{aligned} \mathbb{E} \left[ \left| \int_{\mathbb{R}} \left| \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 - \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q'}, e) \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \right| de \right], \\ \leq \mathbb{E} \left[ \left| \left( \partial_p P_{t_{n+1}}^{t,q'} \right)^2 \int_{\mathbb{R}} \left| \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) - \partial_{pp}^2 \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q'}, e) \right| de \right], \\ \leq \frac{C}{\epsilon^2} |q - q'|, \end{aligned} $$

where we used Lemma 2.1 and (4.52), for the last inequality. Combining (4.61) with the previous inequality, we finally obtain that

$$ B_1 \leq \frac{C}{\epsilon^2} |q - q'| + \frac{C}{\epsilon} |q - q'|. \quad (4.62) $$

For the term $B_2$ given in (4.56), we observe

$$ B_2 \leq \mathbb{E} \left[ \left| \int_{\mathbb{R}} \left| \partial_p \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \partial_{pp}^2 P_{t_{n+1}}^{t,q} - \partial_p \mathcal{V}^\epsilon (t_{n+1}, P_{t_{n+1}}^{t,q}, e) \partial_{pp}^2 P_{t_{n+1}}^{t,q'} \right| de \right] $$
$$ + \mathbb{E} \left[ \left| \int_{\mathbb{R}} |\partial_p \mathcal{V}^\epsilon(t_{n+1}, P_{t_{n+1}}^{t,q}, e) \partial_{pp}^2 P_{t_{n+1}}^{t,q'} - \partial_p \mathcal{V}^\epsilon(t_{n+1}, P_{t_{n+1}}^{t,q'}, e) \partial_{pp}^2 P_{t_{n+1}}^{t,q'} | \right| \mathrm{d}e \right]. $$

Combining Lemma 2.1, the estimates given in (4.53) and (4.52), with similar computations as done previously, we obtain

$$ B_2 \le \frac{C}{\epsilon} |q - q'| + C |q - q'|. \quad (4.63) $$

Inserting (4.62) and (4.63) back into (4.56), we obtain the

$$ \int_{\mathbb{R}} |\hat{\partial}_{pp}^2 w_n(t, q, e) - \hat{\partial}_{pp}^2 w_n(t, q', e)| \mathrm{d}e \le \frac{C}{\epsilon} |q - q'|, \quad (4.64) $$

as $\epsilon \in (0, 1)$.

Similar arguments allow to retrieve the estimates for the first order derivatives. $\square$

**Proposition 4.4.** Under Assumption 2.1, the following holds

$$ \mathcal{E}_n^D(p) \le C \frac{\mathfrak{h}^{\frac{3}{2}}}{\epsilon^2}, $$

for $0 \le n \le N-1$ and $p \in \mathbb{R}^d$.

*Proof.* For ease of presentation, we do the proof in dimension $d=1$.

(1) We first recall Definition 3.1 and (3.6). In this context, we introduce two new discrete processes:

$$ \begin{aligned} \hat{P}_t^{t_{n,p}} &= p + b(p)(t - t_n) + \bar{\sigma} \frac{\sqrt{t - t_n}}{\sqrt{\hbar}} \Delta \hat{W}_n, \\ \bar{P}_t^{t_{n,p}, \lambda} &= p + b(p)(t - t_n) + \lambda \bar{\sigma} \frac{\sqrt{t - t_n}}{\sqrt{\hbar}} \Delta \hat{W}_n, \end{aligned} $$

for $p \in \mathbb{R}$, $t \in [t_n, t_{n+1}]$, $\lambda \in [0, 1]$ and $0 \le n \le N-1$.

Now, recalling (4.42), we first observe

$$ \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, \hat{P}_{t_{n+1}}^{t_{n,p}, e})] - \mathbb{E}[\mathcal{V}^\epsilon(t_{n+1}, P_{t_{n+1}}^{t_{n,p}, e})] = \mathbb{E}[w_n(t_{n+1}, \hat{P}_{t_{n+1}}^{t_{n,p}, e})] - w_n(t_n, p, e). \quad (4.65) $$

We apply the discrete Ito formula given of Proposition 14 in [7] to our much simpler framework, to obtain

$$ \mathbb{E}[w_n(t_{n+1}, \hat{P}_{t_{n+1}}^{t_{n,p}, e})] - w_n(t_n, p, e) = \int_{t_n}^{t_{n+1}} \mathbb{E} \left[ \hat{\partial}_t w_n(t, \hat{P}_t^{t_{n,p}, e}) + b(p) \hat{\partial}_p w_n(t, \hat{P}_t^{t_{n,p}, e}) \right] \mathrm{d}t \\ + \frac{\bar{\sigma}^2}{2} \mathbb{E} \left[ \int_{t_n}^{t_{n+1}} \int_0^1 \partial_{pp}^2 w_n(t, \bar{P}_t^{t_{n,p}, \lambda}, e) \mathrm{d}\lambda \mathrm{d}t \right]. $$

Using the PDE (4.43) satisfied by $w_n$, the previous equality leads to

$$ \begin{aligned} \mathbb{E}[w_n(t_{n+1}, \hat{P}_{t_{n+1}}^{t_{n,p}, e})] - w_n(t_n, p, e) &= \int_{t_n}^{t_{n+1}} \mathbb{E} \left[ (b(p) - b(\hat{P}_t^{t_{n,p}, e})) \hat{\partial}_p w_n(t, \hat{P}_t^{t_{n,p}, e}) \right] \mathrm{d}t \end{aligned} \quad (4.66) $$

$$ + \frac{\bar{\sigma}^2}{2} \mathbb{E} \left[ \int_{t_n}^{t_{n+1}} \int_0^1 \left[ \hat{\partial}_{pp}^2 w_n(t, \bar{P}_t^{t_{n,p}, \lambda}, e) - \hat{\partial}_{pp}^2 w_n(t, \hat{P}_t^{t_{n,p}, e}) \right] \mathrm{d}\lambda \mathrm{d}t \right]. \quad (4.67) $$
(2) For $e \in \mathbb{R}$, we denote by $A^1(e)$, the term in (4.66) and by $A^2(e)$, the term in (4.67). We compute

$$ egin{aligned} |A^1(e)| &\le \int_{t_n}^{t_{n+1}} \mathbb{E}[|(b(p) - b(\hat{P}_t^{t_n,p}))| \|\partial_p w_n(t, \hat{P}_t^{t_n,p}, e)\|] dt, \\ &\le C\mathfrak{h}^{\frac{1}{2}} \int_{t_n}^{t_{n+1}} \mathbb{E}[|\partial_p w_n(t, \hat{P}_t^{t_n,p}, e)|] dt, \end{aligned} $$

and thus

$$ egin{aligned} \int_{\mathbb{R}} |A^1(e)| \, de &\le C\mathfrak{h}^{\frac{1}{2}} \int_{t_n}^{t_{n+1}} \mathbb{E} \left[ \left| \int_{\mathbb{R}} |\partial_p w_n(t, \hat{P}_t^{t_n,p}, e)| \, de \right| \right] dt, \\ &\le C\mathfrak{h}^{\frac{3}{2}}, \end{aligned} \quad (4.68) $$

where we used (4.46) for the last inequality.

Using Proposition 4.3, we compute

$$ egin{aligned} \int_{\mathbb{R}} |\partial_{pp}^2 w_n(t, \bar{P}_t^{t_n,p,\lambda}, e) - \partial_{pp}^2 w_n(t, \hat{P}_t^{t_n,p}, e)| \, de &\le \frac{C}{\epsilon^2} |\bar{P}_t^{t_n,p,\lambda} - \hat{P}_t^{t_n,p}|, \\ &\le \frac{C}{\epsilon^2} \sqrt{\mathfrak{h}}. \end{aligned} $$

This yields

$$ \int_{\mathbb{R}} |A_2(e)| \, de \le \frac{C}{\epsilon^2} \mathfrak{h}^{\frac{3}{2}}. \quad (4.69) $$

Combining the previous inequality with (4.68) and (4.66), (4.67), we get

$$ \int_{\mathbb{R}} |\mathbb{E}[w_n(t_{n+1}, \hat{P}_{t_{n+1}}^{t_n,p}, e)] - w_n(t_n, p, e)| \, de \le \frac{C}{\epsilon^2} \mathfrak{h}^{\frac{3}{2}}. $$

The proof is concluded recalling (4.65) and (4.42).

□

## 4.5. Proof of Theorem 3.1

The proof of our main result is now almost straightforward. From Lemma 4.1, Corollary 4.1 and Proposition 4.4, we observe that

$$ \mathbb{E}[\mathcal{E}_n^T(\hat{P}_{t_n})] \le C \frac{\mathfrak{h}}{M}, \quad \mathbb{E}[\mathcal{E}_n^D(\hat{P}_{t_n})] \le C \frac{\mathfrak{h}^{\frac{3}{2}}}{\epsilon^2} \quad \text{and} \quad \mathbb{E}[\mathcal{E}_n^S(\hat{P}_{t_n})] \le C \mathfrak{h}(\mathfrak{h}^{\frac{1}{2}} + \epsilon). \quad (4.70) $$

Summing over $n$, we get

$$ \sum_{n=0}^{N-2} \mathbb{E}[\mathcal{E}_n^T(\hat{P}_{t_n}) + \mathcal{E}_n^D(\hat{P}_{t_n}) + \mathcal{E}_n^S(\hat{P}_{t_n})] \le C \left( \frac{1}{M} + \frac{\sqrt{\mathfrak{h}}}{\epsilon^2} + \epsilon \right). \quad (4.71) $$

For the last step, we obtain

$$ \mathbb{E}[\mathcal{E}_{N-1}(\hat{P}_{t_{N-1}})] \le C \left( \frac{h}{M} + \epsilon + \mathfrak{h}(\mathfrak{h}^{\frac{1}{2}} + \epsilon) \right) \quad (4.72) $$

combining (3.34) with Proposition 4.1 and Proposition 2.1. The previous inequality combined with (4.71) and Proposition 3.1 yields (3.40).

Balancing optimally the errors in (3.40) allows to conclude by proving (3.41).

□
# 5. NUMERICS

In this section, we realize numerical experiments to test in practice the schemes introduced in Section 3.

## 5.1. Setting

We first introduce the models we will use. They have already been considered in [6], which facilitates the comparison with previous numerical results.

Let us first define the following toy model where the process P corresponds to a Brownian motion and is a multidimensional version of the model given in Remark 2.2:

**Example 5.1 (Linear model).**

$$dP_t = \sigma dW_t, dE_t = \left( \frac{1}{\sqrt{d}} \sum_{\ell=1}^{d} P_t^{\ell} - Y_t \right) dt \quad (5.1)$$

with terminal function $(p, e) \mapsto \phi(p, e) = \mathbf{1}_{\{e \ge 0\}}$ and where W is a d-dimensional Brownian motion and $\sigma > 0$.

We will also consider a multiplicative model:

**Example 5.2 (Multiplicative model).** Let W be a d-dimensional Brownian motion. For all $\ell \in \{1, \dots, d\}$, we set

$$dP_t^\ell = \mu P_t^\ell dt + \sigma P_t^\ell dW_t^\ell, \quad P_0^\ell = 1, \text{ and } dE_t = \tilde{\mu}(P_t, Y_t) dt \quad (5.2)$$

with $(p, y) \mapsto \tilde{\mu}(p, y) = \left(\prod_{\ell=1}^{d} p^{\ell}\right)^{\frac{1}{\sqrt{d}}} e^{-\theta y}$, for some $\theta > 0$. The terminal condition is given by $(p, e) \mapsto \phi(p, e) = \mathbf{1}_{\{e \ge 0\}}$.

As mentioned in [6], one advantage of these models is that they can be reduced to models with lower dimension (at most 2) whatever the value of $d$. For comparison, we will also use some numerical proxy: for the first model in Example 5.1, we use the method considered in [1] applied to the one dimensional reduction of the model, for the second model in Example 5.2 we use the NN & Upwind scheme. Both methods are introduced and discussed in [6].

We also observe that the above models fit into the setting of Section 3.1.2, recall Assumption 3.1. As shown below and as expected, they are thus numerically tractable. We shall use in this regard the scheme introduced in Definition 3.2, linked to CASE 1. We will also consider CASE 2, where the number of particles is capped. To this end, we need to define precisely the operator $\mathfrak{R}^{M,m}$ appearing in Definition 3.3, which is the purpose of the next section.

## 5.2. Definition of $\mathfrak{R}$

There are various possible implementations for the operator $\mathfrak{R}$ whose goal is to tame the computational cost, recall Definition 3.3. Using this operator $\mathfrak{R}$ introduces another numerical error. In view of the convergence study in the previous sections that focuses on $L^1$-norm, one could aim to optimise this error measurement also in the implementation of $\mathfrak{R}$. Namely, with the notations of Definition 3.3 (step 2), this means that we would like to have, for each $0 \le n \le N-1$, $p \in \mathbb{R}^d$,

$$\int_{\mathbb{R}} |\tilde{v}_n^{N,M}(p, e) - \bar{v}_n^{N,M}(p, e)| \, de \quad (5.3)$$

minimal. The solution to this problem is known, see Section 2 of [11]. Indeed, $\tilde{v}_n^{N,M}(p, \cdot)$ and $\bar{v}_n^{N,M}(p, \cdot)$ are CDF of probability measures on $\mathbb{R}$, respectively denoted $\hat{\mu}^{N,M}(p)$ and $\bar{\mu}^{N,M}(p)$ for later use. It turns out that minimising (5.3) corresponds exactly to minimise the $W_1$-distance between the two distributions. Let us recall the general approximation result stated in Section 2 of [11].
**Lemma 5.1.** Let $F$ be the CDF associated to a probability distribution $\nu \in \mathcal{P}_1(\mathbb{R})$. We note $F^{-1}$ its generalized inverse, given by

$$F^{-1}(v) = \inf\{y \in \mathbb{R} | F(y) \ge v\}, v \in (0, 1). \qquad (5.4)$$

Let $\hat{\nu}_M = \frac{1}{M} \sum_{m=1}^{M} \delta_{e_m}$, then the minimum of

$$\mathcal{W}_1(\mu, \hat{\nu}_M) = \int_{\mathbb{R}} |H * \nu(x) - H * \hat{\nu}_M(x)| dx$$

is achieved for $e_m = F^{-1}(\frac{2m-1}{2M})$, $1 \le m \le M$.

In our framework, the above result can be specialized since the starting distribution $\bar{\mu}^{N,M}(p)$ is already atomic (with $2dM$ atoms).

**Corollary 5.1.** For $\xi \in \mathcal{D}_{2dM}$, let $\nu := \frac{1}{2dM} \sum_{m=1}^{2dM} \delta_{e_m}$. Then the optimal $\hat{\nu}_M$, recall Lemma 5.1 is given by

$$\hat{\nu}_M = \frac{1}{M} \sum_{m=1}^{M} \delta_{e_m} \quad \text{with} \quad e_m = \xi_{(2m-1)d}. \qquad (5.5)$$

*Proof.* According to Lemma 5.1 and the definition of $F^{-1}$ in (5.4), one should find, for $1 \le m \le M$,

$$i^*(m) := \min\{i \in \{1, \dots, M\} \mid F(\xi_i) \ge \frac{2m-1}{2M}\}. \qquad (5.6)$$

Observing that $F(\xi_i) = \frac{i}{2dM}$, we obtain that $i^*(m) = (2m-1)d$, which concludes the proof of the corollary. $\square$

**Definition 5.1** (Operator $\mathfrak{R}^{M,m}$). For $M \ge m \ge 1$, we first denote $S^M$ the operator of sorting $M$ particles in ascending order, namely

$$\mathbb{R}^M \ni \xi \mapsto S^M(\xi) \in \mathcal{D}_M. \qquad (5.7)$$

At each step $n$, for a set of $2dM$ particles namely $e \in \mathbb{R}^{2dM}$, $S^{2dM}(e)$ could be written as follows:

$$S^{2dM}(\xi) = (\tilde{S}^{2d,1}(\xi), \dots, \tilde{S}^{2d,M}(\xi)), \qquad (5.8)$$

where for $1 \le i \le M$, $\tilde{S}^{2d,i}(\xi) := (S_j^{2dM}(\xi))_{2d(i-1)<j \le 2d}$, $i \in \mathcal{I}_{2d}$, with $S_i^{2dM}(\xi)$ denoting the $i$-th coordinate of vector $S^{2dM}(\xi)$. For $\psi^\xi \in \mathcal{I}^{2dM}$ associated to $\xi \in \mathcal{D}_{2dM}$, we recall that

$$\psi^\xi = H * \left( \frac{1}{2dM} \sum_{m=1}^{2dM} \delta_{S_m^{2dM}(\xi)} \right),$$

we then define the following operator

$$\psi^\xi \ni \mathcal{I}^{2dM} \mapsto \mathfrak{R}^{2dM,M}(\psi^\xi) := H * \left( \frac{1}{M} \sum_{m=1}^{M} \delta_{\tilde{e}_m} \right) \in \mathcal{I}^M, \qquad (5.9)$$

where four different choices for $\tilde{e}_m$ are considered, leading to four different implementations. Namely, for each $1 \le m \le M$,
<table><caption>TABLE 1. L1-error and L∞-error for model Example 5.1 with different numbers of particles with respect to leftmost, mean, rightmost, optim, note that the time steps N = 20.</caption><thead><tr><th>Sigma</th><th>Number of particles</th><th>Method</th><th>L1-error</th><th>L∞-error</th></tr></thead><tbody><tr><td rowspan="4">1.0</td><td rowspan="4">1000</td><td>leftmost</td><td>5.71e-3</td><td>1.67e-2</td></tr><tr><td>mean</td><td>4.96e-3</td><td>1.50e-2</td></tr><tr><td>rightmost</td><td>5.39e-3</td><td>1.80e-2</td></tr><tr><td>optim</td><td>3.72e-3</td><td>1.02e-2</td></tr><tr><td rowspan="4">5000</td><td rowspan="4"></td><td>leftmost</td><td>5.05e-3</td><td>1.45e-2</td></tr><tr><td>mean</td><td>4.92e-3</td><td>1.46e-2</td></tr><tr><td>rightmost</td><td>4.92e-3</td><td>1.53e-2</td></tr><tr><td>optim</td><td>3.47e-3</td><td>1.01e-2</td></tr><tr><td rowspan="4">0.3</td><td rowspan="4">1000</td><td>leftmost</td><td>1.34e-3</td><td>0.85e-2</td></tr><tr><td>mean</td><td>0.88e-3</td><td>0.61e-2</td></tr><tr><td>rightmost</td><td>1.25e-3</td><td>0.91e-2</td></tr><tr><td>optim</td><td>1.19e-3</td><td>0.69e-2</td></tr><tr><td rowspan="4">5000</td><td rowspan="4"></td><td>leftmost</td><td>0.94e-3</td><td>0.68e-2</td></tr><tr><td>mean</td><td>0.88e-3</td><td>0.61e-2</td></tr><tr><td>rightmost</td><td>0.92e-3</td><td>0.65e-2</td></tr><tr><td>optim</td><td>0.56e-3</td><td>0.46e-2</td></tr><tr><td rowspan="4">0.01</td><td rowspan="4">1000</td><td>leftmost</td><td>0.27e-3</td><td>0.39e-2</td></tr><tr><td>mean</td><td>0.13e-3</td><td>0.20e-2</td></tr><tr><td>rightmost</td><td>0.48e-3</td><td>0.43e-2</td></tr><tr><td>optim</td><td>0.11e-3</td><td>0.17e-2</td></tr><tr><td rowspan="4">5000</td><td rowspan="4"></td><td>leftmost</td><td>0.09e-3</td><td>0.13e-2</td></tr><tr><td>mean</td><td>0.08e-3</td><td>0.17e-2</td></tr><tr><td>rightmost</td><td>0.13e-3</td><td>0.18e-2</td></tr><tr><td>optim</td><td>0.06e-3</td><td>0.11e-2</td></tr></tbody></table>

mean: The mean position of 2d particles of $\tilde{S}^{2d,m}(\xi)$, namely

$$ \tilde{e}_m := \frac{1}{2d} \sum_{j=2d(m-1)+1}^{2dm} S_j^{2dM}(\xi). \qquad (5.10) $$

leftmost: The minimum position among the 2d particles of $\tilde{S}^{2d,m}(\xi)$ is kept,

$$ \tilde{e}_m := \min\{S_j^{2dM}(\xi), 2d(m-1) + 1 \le j \le 2dm\} = S_{2d(m-1)+1}^{2dM}(\xi). \qquad (5.11) $$

rightmost: The maximum position among the 2d particles of $\tilde{S}^{2d,m}(\xi)$ is kept,

$$ \tilde{e}_m := \max\{S_j^{2dM}(\xi), 2d(m-1) + 1 \le j \le 2dm\} = S_{2dm}^{2dM}(\xi). \qquad (5.12) $$

optim: The optimal position according to Corollary 5.1 is kept, $\tilde{e}_m := S_{d(2m-1)}^{2dM}(\xi)$.

## 5.3. Numerical results

First of all, to validate empirically our different implementations of operator $\mathfrak{R}$, we report the L1-error and L∞-error against the “proxy” solution for different numbers of particles, based on the different strategies to reduce the number of particles, see Definition 5.1 for $\sigma = 0.01, 0.3, 1.0$, see Table 1. We observe that the four implementations of $\mathfrak{R}$: leftmost, mean, rightmost, optim of positions of particles, yield practically the same
<table><caption>TABLE 2. Computational cost in Example 5.1 for different dimension *d* (for the *P*-variable) for CASE 1 and CASE 2 schemes.</caption><thead><tr><th>Dimension</th><th>d = 1</th><th>d = 4</th></tr></thead><tbody><tr><td>Time for CASE 1</td><td>46.96 s</td><td>784 s</td></tr><tr><td>Time for CASE 2</td><td>0.3 s</td><td>1.3 s</td></tr></tbody></table>

<figure><img src="image_1.png" /><figcaption>FIGURE 1. Model of Example 5.1: Comparison of the two methods CASE 1 & CASE2 with *d* = 4. The Proxy solution is given by the same particle method on the one-dimensional PDE. For *BT&SPD*: both CASE 1 and CASE 2, the number of particles is *M* = 3500 and the number of time steps *N* = 20. (a) *σ* = 0.01. (b) *σ* = 0.3. (c) *σ* = 1.0.</figcaption></figure>

numerical results for different levels of volatilities. The *optim* strategy appears to be better, as expected. Hence in the following, we only consider the *CASE 2* implemented with operator $\Re \text{ optim}$. The clear advantage of considering the *CASE 2* scheme is that the number of particles does not grow exponentially during the iterations, hence it takes less computational time compared to *CASE 1*, see Table 2.

To validate empirically our numerical schemes, we plot the value function $\mathcal{V}$ of model Example 5.1 using scheme *CASE 1* and *CASE 2* against the *proxy* solution presented in [1], see Figure 1. As expected, both *CASE 1* and *CASE 2* scheme could reproduce correctly the entropy solution of the PDE (2.1). We also observe that the solutions obtained by *CASE 1* and *CASE 2* scheme are quite close.

Apart from the linear model, we have also tested numerically our *CASE 1* and *CASE 2* schemes in model Example 5.2 against a NN & Upwind method based on splitting schemes presented in [6], see Figure 2. We
<figure><img src="image_2.png" /><figcaption>FIGURE 2. Model of Example 5.2: Comparison of the two methods CASE 1 & CASE2 with $d = 4$. The Proxy solution is given by the Neural nets & Upwind in [6]. For *BT&SPD*: both CASE 1 and CASE 2, the number of particles is $M = 3500$ and the number of time steps $N = 20$. (a) $\sigma = 0.01$. (b) $\sigma = 0.3$. (c) $\sigma = 1.0$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 3. Convergence rate on $h$ for the model in Example 5.1 with parameters $d = 4, \sigma = 1.0$.</figcaption></figure>

note that the function $\mu$ in Example 5.2 is always non negative, thus we use the Upwind scheme which is less “diffusive” than the Lax-Friedrichs scheme as argued in [6]. As expected, both CASE 1 and CASE 2 schemes give satisfying numerical results for different levels of volatility.

At last, we want to empirically estimate the convergence rate of the error introduced by our numerical scheme. We consider the model Example 5.1 where $\sigma = 1.0$. We consider a set of number of time steps $N :=$
{2, 4, 8, 16, 32, 64} and time step $\Delta t := \frac{T}{N}$, and compute the L1-error for BT & SPD method. Note that the proxy solution is always given by method in [1] applied to one dimensional equivalent model. The empirical convergence rate with respect to the time step is close to one see Figure 3, which is much better than the order $\frac{1}{6}$ obtained in Theorem 3.1.

## ACKNOWLEDGMENTS

The authors would like to thank Benjamin Jourdain for fruitful discussions.

## DATA AVAILABILITY STATEMENT

No new data/codes were created or analyzed in this study.

## REFERENCES

[1] M. Bossy and D. Talay, Convergence rate for the approximation of the limit law of weakly interacting particles: application to the burgers equation. *Ann. Appl. Probab.* **6** (1996) 818–861.

[2] Y. Brenier and E. Grenier, Sticky particles and scalar conservation laws. *SIAM J. Numer. Anal.* **35** (1998) 2317–2328.

[3] R. Carmona and F. Delarue, Singular FBSDEs and scalar conservation laws driven by diffusion processes. *Probab. Theory Relat. Fields* **157** (2013) 333–388.

[4] R. Carmona, F. Delarue, G.-E. Espinosa and N. Touzi, Singular forward-backward stochastic differential equations and emissions derivatives. *Ann. Appl. Probab.* **23** (2013) 1086–1128.

[5] J.-F. Chassagneux and C. Garcia Trillos, Cubature method to solve BSDEs: error expansion and complexity control. *Math. Comput.* **89** (2020) 1895–1932.

[6] J.-F. Chassagneux and M. Yang, Numerical approximation of singular forward-backward SDEs. *J. Comput. Phys.* **468** (2022) 111459.

[7] J.-F. Chassagneux, D. Crisan and F. Delarue, Numerical method for FBSDEs of Mckean-Vlasov type. *Ann. Appl. Probab.* **29** (2019) 1640–1684.

[8] J.-F. Chassagneux, H. Chotai and D. Crisan, Modelling multiperiod carbon markets using singular forward-backward SDEs. *Math. Oper. Res.* **48** (2023) 463–497.

[9] L.G. Gyurkó and T.J. Lyons, Efficient and practical implementations of cubature on wiener space, in Stochastic Analysis 2010. Springer (2011) 73–111.

[10] S. Howison and D. Schwarz, Risk-neutral pricing of financial instruments in emission markets: a structural approach. *SIAM J. Finan. Math.* **3** (2012) 709–739.

[11] B. Jourdain and J. Reygner, Optimal convergence rate of the multitype sticky particle approximation of one-dimensional diagonal hyperbolic systems with monotonic initial data. *Discrete Contin. Dyn. Syst.* – *A* **36** (2016).

[12] C. Litterer and T. Lyons, High order recombination and an application to cubature on Wiener space. *Ann. Appl. Probab.* **22** (2012) 1301–1327.

[13] J. Zhang, *Some fine properties of backward stochastic differential equations*. Ph.D. thesis, Purdue University (2001).

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the *Subscribe to Open model* (S2O). We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.