2.1. Fixed-point equation for $\bar{\gamma}$

We introduce a closely related SDE to (2), which takes function $\gamma$ as input. Namely, given $\gamma\in L^{\infty}([0,T],\mathbb{R}^K) $ denote by $Z^{\gamma}$ the solution to the Markovian SDE,

(5) \begin{equation} \textrm{d} Z_t = \gamma(t)\big( \alpha(t,Z_t) \, \textrm{d} t + \beta(t,Z_t) \, \textrm{d} W_t \big),\quad Z_0=\xi.\end{equation}

The following result is a direct application of Lemma 3.

Proposition 3. Let Assumptions 1 and 2 hold true for some $R>0$ . Then (5) is strongly well-posed for any $\gamma\in L^{\infty}([0,T],\mathbb{R}^K) $ . Furthermore, for any $\kappa>0$ , there exists $C>0$ dependent on $\kappa, T,R$ , and K, such that

\begin{align*} \lVert Z^{\gamma} - Z^{\widetilde{\gamma}} \rVert_{\mathcal{S}^2_{[0,T]}} \leq C \|\gamma - \widetilde{\gamma}\|_{\infty}, \end{align*}

for any $\gamma,\widetilde{\gamma} \in L^{\infty}([0,T],\mathbb{R}^K)$ such that $\|\gamma\|_{\infty},\|\widetilde{\gamma}\|_{\infty} \leq \kappa$ .

Now, define the operator $\Psi:C\big([0,T],\mathbb{R}^K) \to C\big([0,T],\mathbb{R}^K)$ given as

\begin{equation*} \Psi(\gamma)(t) \,:\!=\, \mathrm{E} [ \varphi(Z^{\gamma}_t)] , \quad t\in[0,T],\end{equation*}

which is well-defined under Assumptions 1 and 2.

Let us start by observing that solving the MV-SDE (2) is equivalent to finding a fixed point for $\Psi$ . Indeed, if $\bar{\gamma}$ is a solution to

(6) \begin{equation}\gamma = \Psi( \gamma) ,\end{equation}

then the process $Z^{\bar{\gamma}}$ solves (2) and vice versa. Also note that, if $\bar{\gamma}$ solves (6), then $\bar{\gamma}$ solves the minimization problem

(7) \begin{equation}\min_{\gamma\in C([0,T])} F^2(\gamma),\end{equation}

where

(8) \begin{equation}F \,:\, C\big([0,T],\mathbb{R}^K\big) \to \mathbb{R}, \quad F(\gamma)\,:\!=\, \| \Psi(\gamma) - \gamma \|_{L^2_{([0,T])}}.\end{equation}

Therefore, under Assumptions 1 and 2, Proposition 1 implies that there exists a unique $\bar{\gamma} \in C\big([0,T],\mathbb{R}^K\big)$ that solves (7), in particular,

(9) \begin{equation}\min_{\gamma\in C([0,T])} F^2(\gamma) = F^2(\bar{\gamma}) = 0.\end{equation}

Lemma 1. Let Assumptions 1 and 2 hold true for some $R>0$ . Then, the functional F in (7) is locally Lipschitz continuous in sup-norm $\|\cdot\|_\infty$ . Precisely, for any $\kappa > 0$ there exists $C > 0$ such that

(10) \begin{align}| F(\gamma) - F(\hat\gamma)| \leq C \lVert \gamma-\hat\gamma \rVert_{\infty}, \quad \gamma, \hat\gamma \in C([0,T],\mathbb{R}^K), \quad \textrm{with } \|\gamma\|_{\infty},\|\hat\gamma\|_{\infty} \leq \kappa,\end{align}

where $C>0$ depends on $\kappa, T,R$ , and K.

Proof. Given $\kappa>0$ , take $\gamma, \hat\gamma \in C([0,T],\mathbb{R}^K) \ \textrm{with } \|\gamma\|_{\infty},\|\hat\gamma\|_{\infty} \leq \kappa$ . By the reverse triangle inequality we have

(11) \begin{equation}\Big|\lVert \Psi(\gamma) - \gamma \rVert_{L^2([0,T])} - \lVert \Psi(\hat\gamma) - \hat\gamma \rVert_{L^2([0,T])} \Big| \leq \lVert \Psi(\gamma) - \Psi(\hat\gamma) \rVert_{L^2([0,T])} + \lVert \gamma - \hat\gamma \rVert_{L^2([0,T])}.\end{equation}

The second term in the right-hand side of (11) is dominated by $\sqrt{T} \lVert \gamma- \hat\gamma\rVert_{\infty}.$ For the first term, by the definition of $\Psi$ and by the Lipschitz continuity of $\varphi$ , we obtain

\begin{align*} \lVert \Psi(\gamma) - \Psi(\hat\gamma) \rVert^2_{L^2([0,T])} =& \int_0^T\Big| \mathrm{E}\big[\varphi(Z^\gamma_t) - \varphi(Z^{\hat\gamma}_t)\big]\Big|^2 \, \textrm{d} t \leq C \int_0^T\mathrm{E}\Big[\big|Z^\gamma_t - Z^{\hat\gamma}_t\big|^2\Big] \, \textrm{d} t\\ \leq& C\, T \lVert \gamma-\hat\gamma \rVert^2_{\infty}, \end{align*}

where the last inequality follows from Proposition 3. Thus, by combining the two bounds for the terms in the right-hand side of (11), we get the result in (10).

A quick inspection of the above proof shows that (10) also holds with $\|\cdot\|_\infty$ replaced by $\|\cdot\|_{L^2([0,T])}$ since the same is also true for Proposition 3.

Remark 5. Depending on the numerical accuracy requirements, instead of $F(\gamma)$ defined above, we could also use an objective function $\widehat F(\gamma)$ given as

\begin{align*}\widehat F^2(\gamma) = \int_0^T w(s)|\Psi(\gamma)(s) - \gamma(s)|^2 \, \textrm{d}s,\end{align*}

with a kernel (weight) function w over [0,T]. For example, we can take $w(t)=C_1e^{C_2 t}$ for $C_1,C_2>0$ . Such a choice of objective function $\widehat F(\gamma)$ in the minimization problem (7) would force the SGD algorithm to be more accurate towards the end of the interval [0,T].

2.2. Projection on a finite-dimensional domain

In order to solve (7) numerically using a gradient-descent method, it is imperative that we first narrow down the infinite-dimensional domain $C([0,T],\mathbb{R}^{K})$ to a finite-dimensional one. This approach is a natural choice in our setting of solving (7) as the solution $\bar \gamma$ is a continuous map. Hereafter, we fix $n\in\mathbb{N}$ and a family $g_0,\ldots,g_n$ of linearly independent vectors of $C([0,T],\mathbb{R})$ . Define the finite-dimensional sub-space of $C([0,T],\mathbb{R}^{K})$ as

(12) \begin{equation} S_n \,:\!=\, \{ a_0 g_0 + \cdots + a_n g_n\,:\, a=(a_0, \ldots, a_n) \in \mathbb{R}^{(n+1) K } \},\end{equation}

and the distance between $\gamma \in C([0,T],\mathbb{R}^{K})$ and $S_n$ as

\begin{equation*}d_n(\gamma) \,:\!=\, \inf\{ \| \gamma - p \|_{\infty} \,:\, p \in S_n \}.\end{equation*}

Also introduce a linear lifting operator

(13) \begin{equation} \mathscr{L} \,:\, \mathbb{R}^{(n+1)\times K} \to S_n, \quad \mathscr{L} a = \mathscr{L} (a_0, \ldots, a_n) \,:\!=\, a_0g_0+ \cdots + a_n g_n.\end{equation}

Remark 6. The operator $\mathscr{L}$ is an isomorphism between finite-dimensional normed vector spaces. In particular, both $\mathscr{L}$ and its inverse $\mathscr{L}^{-1}$ are continuous and bounded operators. Furthermore, by a compactness argument, it can be easily seen that for any $\gamma \in C([0,T],\mathbb{R}^K),$ there exists $p_n(\gamma) \in S_n$ such that

(14) \begin{equation} \| \gamma - p_n(\gamma) \|_{\infty} = d_n(\gamma).\end{equation}

In other words, $d_n$ is well-defined as a minimum and such a minimum is attained at $p=p_n(\gamma)$ . It is important to observe that both $p_n(\gamma)$ and $d_n(\gamma)$ only depend on $S_n$ and not on the particular choice of the basis $g_i$ , $i=0,\ldots,n$ .

Example 3. (Polynomial of best approximation.) Consider $n\in\mathbb{N}$ and $g_i (t) = t^i$ for $i=0,\ldots, n$ . Then the elements of $S_n$ read as

\begin{equation*} p(t) = \sum_{i=0}^n a_i t^i, \quad t\in[0,T].\end{equation*}

This choice clearly leads to $S_n = \mathcal{P}_n(\mathcal{R}^K)$ , the space of $\mathbb{R}^K$ -valued polynomials of order n. We stress that $S_n = \mathcal{P}_n(\mathcal{R}^K)$ could be also obtained by fixing a basis $g_i$ , $i=0,\ldots,n$ , which is not the canonical one. Possible choices include Legendre, Chebyshev, Laguerre, and Hermite polynomials [Reference Süli and Mayers43]. In particular, the numerical results presented in Section 4 refer to the specific choice of $g_i$ given by the Lagrange polynomials centred at Chebyshev points, i.e.

(15) \begin{equation}g_i(t) = \Bigg( \prod_{\substack{ 0\leqslant l \leqslant n \\ l\neq i}} \frac{t-t_l}{t_i - t_l} \Bigg), \quad t_i = \frac{T}{2} + \frac{T}{2} \cos\Bigl( \frac{2i+1}{2n+2}\pi\Bigr),\quad i=0,1,\ldots,n.\end{equation}

Finally note that, in particular, the distance $d_n(\gamma)$ and the projection $p_n(\gamma)$ (see Remark 6) are the same for all these choices of basis, for they all generate the same finite-dimensional space $S_n = \mathcal{P}_n(\mathcal{R}^K)$ . In this case, $p_n$ in (14) is called the polynomial of best approximation of $\gamma$ (see [Reference Isaacson and Keller34, Chapter 5]).

Since the maps $\gamma$ are bounded by $\|\varphi\|_\infty$ , we introduce the auxiliary map $\textbf{ h}\,:\,\mathbb{R}^{K}\to\mathbb{R}^K$ so that we can ensure that $\textbf{h}(\mathscr{L} a)$ is also bounded. Concretely, let $\textbf{h}$ be a smooth bounded function $\mathbb{R}^{K}\ni x \mapsto \textbf{h}(x)=\big(\textbf{h}_1(x_1),\ldots,\textbf{h}_K(x_K)\big)$ defined as

(16) \begin{align}\mathbb{R}^K \ni x \mapsto \textbf{h}(x)=\begin{cases}x ,& \textrm{if}\ |x| \leq \| \varphi \|_{\infty} + 2 d_n (\bar{\gamma})\\[4pt]\| \varphi \|_{\infty} + 2 d_n (\bar{\gamma}),& \textrm{if}\ |x| \geq \| \varphi \|_{\infty} + 2 d_n (\bar{\gamma})\end{cases},\end{align}

where it is clear that $\|\textbf{h}\|_{\infty} \leq \| \varphi \|_{\infty} + 2 d_n (\bar{\gamma})$ .

Let the constant $\tau_n>0$ be such that

\begin{equation*}B_{\tau_n}\supset \mathscr{L}^{-1}\big( B_{\| \varphi \|_{\infty} + d_n (\bar{\gamma})} \big),\end{equation*}

where $B_{\tau_n}$ and $B_{\| \varphi \|_{\infty} + d_n (\bar{\gamma})}$ denote the balls in $\mathbb{R}^{(n+1) K}$ and $S_n$ , respectively, with radius $\tau_n$ and $\| \varphi \|_{\infty} + d_n (\bar{\gamma})$ , centred at zero. Note that such $\tau_n$ exists as $\mathscr{L}^{-1}$ is a linear, and thus a bounded operator (see Remark 6). Finally, we consider a non-negative real-valued function $H\in C^1(\mathbb{R}^{(n+1) K})$ such that

(17) \begin{equation}\textrm{supp }H \subset B^c_{\tau_n}\quad\textrm{and}\quad \lim_{|a|\to + \infty} H (a)=+\infty.\end{equation}

Example 5. A function H satisfying hypothesis (17) above is given by

\begin{equation*}H(a)\,:\!=\,\begin{cases}0, & \textrm{if } a \in B_{\tau_n}\\(|a|-\tau_n)^2, & \textrm{otherwise}\end{cases} , \quad {\nabla H(a)\,:\!=\,\begin{cases}0, & \textrm{if } a \in B_{\tau_n}\\\frac{2 (|a| - \tau_n)}{|a|}a, & \textrm{otherwise}\end{cases} ,}\end{equation*}

Remark 7. Let $\bar{\gamma} \in C([0,T],\mathbb{R}^K)$ be the unique minimizer to (7) and consider the element $p_n(\bar{\gamma})$ introduced in Remark 6. Since $\|\bar{\gamma} \|_{\infty} \leq \| \varphi \|_{\infty}$ , it is clear that $\|p_n(\bar{\gamma})\|_{\infty}\leq \| \varphi \|_{\infty} + d_n(\bar{\gamma})$ . Therefore, by (16) and (17), we obtain that

\begin{equation*}\textbf{h}\big( p_n(\bar{\gamma}) \big) = p_n(\bar{\gamma})\quad\textrm{and}\quad H\big(\mathscr{L}^{-1} p_n(\bar{\gamma})\big)\equiv 0.\end{equation*}

In place of (7), we consider the following optimization problem:

(18) \begin{equation} \min_{a\in \mathbb{R}^{(n+1) K}} G(a) \quad\textrm{with}\quad G(a)\,:\!=\, F^2(\textbf{h}\circ \mathscr{L} a) + H (a).\end{equation}

The role of map $\mathbf{h}$ is to ensure that the estimate of $\bar{\gamma}$ , defined in (9), remains bounded when computing the loss function F. This is reflected when we compute $F^2(\textbf{h}\circ \mathscr{L} a)$ in (18). Moreover, the function H penalizes the coefficient vector a when the norm of our approximation $\mathcal{L}a$ grows beyond $\| \varphi \|_{\infty} + d_n (\bar{\gamma}).$ This is reflected in the definition of H in (17) and also motivated in Remark 7.

Lemma 2. Let Assumptions 1 and 2 hold true for some $R>0$ . Let $\bar{\gamma} \in C([0,T],\mathbb{R}^K)$ be the unique solution to (7) and assume that $\| \varphi \|_{\infty} + d_n(\bar{\gamma})\leq \kappa$ for some $\kappa>0$ . Then the minimization problem (18) admits a solution. Furthermore, for any

\begin{equation*}\bar a_n \in \mathrm{arg\,min}_{a\in \mathbb{R}^{(n+1) K}} G(a)\end{equation*}

we have

(19) \begin{equation}F(\mathscr{L} \bar a_n) \leq C\, d_n(\bar \gamma) ,\end{equation}

where $C>0$ depends on $\kappa, T,R$ , and the dimensions.

Proof. The boundedness of $\varphi$ implies that $F^2(\textbf{h}\circ \mathscr{L} a)$ is bounded. Thus, by (17), G(a) tends to $+\infty$ as $|a| \to \infty$ . On the other hand, Lemma 1 together with the continuity of $\textbf{h}$ and H imply that G is continuous on $\mathbb{R}^{(n+1) K}$ . Therefore, the fact that G has at least a minimum $\bar a_n$ on $\mathbb{R}^{(n+1) K}$ stems from Weierstrass’ extreme value theorem.

Furthermore, let $\bar a_n$ be a minimizer of G we then have

\begin{align*} F^2(\textbf{h}\circ \mathscr{L} \bar a_n)\leq F^2(\textbf{h}\circ \mathscr{L}\bar a_n) + H_n(\mathscr{L}\bar a_n)& \leq F^2\big(\textbf{h}\circ p_n(\bar{\gamma})\big) + H \big(\mathscr{L}^{-1} p_n(\bar{\gamma})\big)\\ & = F^2\big(p_n(\bar{\gamma})\big) = \big|F\big( p_n(\bar{\gamma})\big) - F(\bar \gamma)\big|^2\\ & \leq C \, \lVert p_n(\bar{\gamma}) - \bar \gamma \rVert^2_{\infty} = C \, d^2_n(\bar \gamma). \end{align*}

In the above, we used Remark 7 and Lemma 1. The result follows since $\| \bar{\gamma} \|_{\infty}\leq \| \varphi \|_{\infty}$ and $\| p_n(\bar{\gamma}) \|_{\infty}\leq \| \varphi \|_{\infty} + d_n(\bar{\gamma})\leq \kappa$ .

From the above result we observe that if we can compute the distance $d_n$ explicitly, we have an explicit upper bound for the error encountered when solving (18) instead of (7). For instance, when choosing $S_n$ as the space of the nth-order polynomials (see Example 3), an explicit error estimate can be obtained as shown next.

Proposition 4. Let Assumptions 1–3 hold true for some $N\in\mathbb{N}$ , $R>0$ . Let $g_0,\ldots,g_n$ be a family of linearly independent vectors of $C([0,T],\mathbb{R})$ generating the space of the $\mathbb{R}^K$ -valued polynomials of order n, i.e. $S_n = \mathcal{P}_n(\mathcal{R}^K)$ in (12). Then the minimization problem (18) admits a solution. Furthermore, for any

\begin{equation*}\bar a_n \in \mathrm{arg\,min}_{a\in \mathbb{R}^{(n+1) K}} G(a)\end{equation*}

we have

\begin{equation*}F(\mathscr{L}\bar a_n) \leq \frac{C}{n^N},\end{equation*}

where $C>0$ depends on T,N,R, and the dimensions (and is independent of n).

Proof. By Lemma 2 the minimization problem (18) admits a solution and (19) holds true. By Proposition 2, $\bar{\gamma}$ is N-times differentiable on [0, T]. In particular, (4) implies that its derivatives are bounded on [0, T] by a constant that only depends on N, R, and the dimensions.

A theorem of Jackson (see [Reference Timan49, Section 5.1.5 (p. 261)]) then yields

\begin{equation*}d_n(\bar{\gamma}) \leq \frac{C}{n^N} ,\end{equation*}

which concludes the proof.

3.1. Some preliminaries

To explain the SGD algorithm we begin with some heuristics. Hereafter we write a generic element of $a\in\mathbb{R}^{(n+1)K}$ as

\begin{equation*}a=(a_{0,1}, \ldots, a_{0,K},\ldots,a_{n,1}, \ldots, a_{n,K} ).\end{equation*}

• • For any $k = 0, \ldots, n,$ and $j = 1,\ldots, K$ , by formally differentiating (21) with respect to the element $a_{k,j}$ , we obtain that the pair $(Z^a,\partial_{a_{k,j}} Z^a)$ solves the system given by (21) and

  (22) \begin{align}\textrm{d} Y^{k,j}_t &= g_k(t) \nabla \textbf{h}_j\big( \mathscr{L} a(t) \big) \Big( \alpha(t, Z_t )\, \textrm{d} t + \beta(t, Z_t ) \,\textrm{d} W_t \Big) \nonumber\\&+ \sum^d_{i=1} Y^{k,j,i}_t \textbf{h}\big( \mathscr{L} a(t) \big) \Big( \partial_{z_i} \alpha(t, Z_t ) \, \textrm{d} t + \partial_{z_i} \beta (t, Z_t)\, \textrm{d} W_t \Big), \quad Y^{k,j}_0 = 0, \end{align}
  with $\textbf{h}$ introduced in (16). Note that $Y^{k,j}$ and Z take values in $\mathbb{R}^d$ , with components denoted by $Y^{k,j,i}$ and $Z^i$ , for $i=1,\ldots,d$ , respectively. In the following, we will sometimes reinforce the notation by denoting $(Z^{a}_t (\xi,W),Y^{a;k,j}_t (\xi,W))$ as the strong solution to (21), (22) with respect to a given initial random variable $\xi$ and Brownian motion W.

• • Suppose that we can always exchange derivatives with time integrals and expected values (we justify this exchange of operations below). For any $a\in\mathbb{R}^{(n+1) K}$ , and for any $k=0,\ldots,n$ and $j=1,\ldots, K$ , we obtain

  (23) \begin{align}\nonumber &\partial_{a_{k,j}} \bigg( \int_0^T \big| \mathrm{E} \big[ \varphi(Z^{a}_t) - \textbf{h}\big((\mathscr{L} a) (t)\big) \big] \big|^2 \,\textrm{d} t \bigg) \\&= 2 \int_0^T \Big\langle \mathrm{E}\Big[ \varphi(Z^{a}_t) - \textbf{h}\big((\mathscr{L} a) (t)\big) \Big] , \mathrm{E} \Big[ \partial_{a_{k,j}} \Big( \varphi(Z^{a}_t) - \textbf{h}\big((\mathscr{L} a) (t)\big) \Big) \Big] \Big\rangle \, \textrm{d} t.\end{align}

Let now W, $\tilde W$ be two independent Brownian motions and $\xi$ , $\tilde \xi$ be two identically distributed independent random variables, which are also independent of $(W,\tilde W)$ . Assuming weak-uniqueness for the solutions to (21) we have

(24) \begin{align}&\Big\langle \mathrm{E}\Big[ \varphi(Z^{ a}_t)- \textbf{h}\big((\mathscr{L} a) (t)\big) \Big] , \mathrm{E} \Big[ \partial_{a_{k,j}} \Big( \varphi(Z^{a}_t) - \textbf{h}\big((\mathscr{L} a) (t)\big) \Big) \Big] \Big\rangle \nonumber\\ &= \Big\langle \mathrm{E}\Big[ \varphi \big(Z^{a}_{t}(\xi, W)\big) - \textbf{h}\big((\mathscr{L} a) (t)\big)\Big] , \mathrm{E} \Big[ \partial_{a_{k,j}} \Big( \varphi \big(Z^{a}_{t}(\tilde\xi,\tilde W)\big) - \textbf{ h}\big((\mathscr{L} a) (t)\big) \Big) \Big] \Big\rangle \nonumber\\ &= \mathrm{E}\Big[ \Big\langle \varphi \big(Z^{a}_{t}(\xi, W)\big) - \textbf{h}\big((\mathscr{L} a) (t)\big) , \partial_{a_{k,j}} \Big( \varphi \big(Z^{a}_{t}(\tilde\xi,\tilde W)\big) - \textbf{h}\big((\mathscr{L} a) (t)\big) \Big) \Big\rangle \Big].\end{align}

• • Therefore, combining (23) and (24), and since $\partial_{a_{k,j}} Z^a$ solves (22), we formally obtain

  \begin{equation*} \nabla_a G(a) = \mathrm{E}[ v(a; \xi, {W} ; \, \tilde\xi, \tilde W) ] + \nabla_a H(a),\end{equation*}
  with $v(a; \xi, {W} ; \, \tilde\xi, \tilde W)\in\mathbb{R}^{(n+1)K}$ , whose components are given by
  (25) \begin{align}\nonumber&v_{k,j}(a; \xi, {W} ; \, \tilde\xi, \tilde W) \\& \,:\!=\, 2 \int_0^T \Big\langle \varphi \big(Z^{a}_{t}(\xi, W)\big) - \textbf{h}\big((\mathscr{L} a) (t)\big) , \nabla_x \varphi \big(Z^{a}_{t}(\tilde\xi,\tilde W)\big) Y^{a;k,j}_{t}(\tilde\xi,\tilde W) - \partial_{a_{k,j}} \textbf{ h}\big( (\mathscr{L} a) (t) \big) \Big\rangle \, \textrm{d} t,\end{align}
  where $\nabla_x \varphi$ stands for the Jacobian matrix of $\varphi\,:\,\mathbb{R}^d \to \mathbb{R}^K$ .

3.2. Algorithm

In the previous section, we presented heuristics on how certain computations can be performed in order to implement the stochastic gradient algorithm. With the help of those explanations, we can now introduce the algorithm. First, let:

• • $(\eta_m)_{m\in\mathbb{N}}$ be a deterministic sequence of positive scalars (learning rates) such that

  (26) \begin{equation}\sum_{m=0}^{\infty} \eta_m = + \infty\quad \textrm{and} \quad \sum_{m=0}^{\infty} \eta^2_m < + \infty ;\end{equation}

• • $( \bar{\Omega} , \bar{\mathcal{F}} , (\bar{\mathcal{F}}_m)_{m\in\mathbb{N}} , \bar{\mathrm{P}})$ be a filtered probability space rich enough to support four independent samples $(W_m)_{m\in\mathbb{N}}$ , $(\tilde W_m)_{m\in\mathbb{N}}$ , $(\xi_{m})_{m \in\mathbb{N}}$ , and $(\tilde\xi_{m})_{m \in\mathbb{N}}$ such that

• • $(W_m)_{m\in\mathbb{N}}$ and $(\tilde W_m)_{m\in\mathbb{N}}$ are sequences of independent q-dimensional standard Brownian motions adapted to $(\bar{\mathcal{F}}_m)_{m\in\mathbb{N}}$ ; namely, for any $m\in\mathbb{N}$ , $W_m = (W_{m,t})_{t\in[0,T]}$ is a q-dimensional standard Brownian motion such that $\mathcal{F}^{W_m}_T \subset \bar{\mathcal{F}}_m$ , and the same holds for $\tilde W_m$ ;

• • $(\xi_{m})_{m \in\mathbb{N}}$ and $(\tilde\xi_{m})_{m \in\mathbb{N}}$ of $\xi$ are sequences of independent $\mathbb{R}^d$ -valued random variables distributed like $\xi$ .

For a deterministic $\textbf{a}_0 \in \mathbb{R}^{(n+1) K},$ the estimate update is obtained iteratively as

\begin{equation*}\textbf{a}_{m+1} \,:\!=\, \textbf{a}_m - \eta_m \textbf{v}_{m+1}, \quad \textbf{v}_{m+1} = v(\textbf{a}_m; \xi_{m+1}, W_{m+1};\tilde\xi_{m+1}, \tilde W_{m+1}) + \nabla H(\textbf{a}_{m}), \quad m\in \mathbb{N}_0.\end{equation*}

The function v is defined component-wise by (25). Note that the $\mathbb{R}^{(n+1) K}$ -valued stochastic processes $(\textbf{a}_m)_{m\in\mathbb{N}}$ and $(\textbf{v}_m)_{m\in\mathbb{N}}$ are recursively defined, as in Algorithm 1, and are adapted to the filtration $(\bar{\mathcal{F}}_m)_{m\in\mathbb{N}}$ , which then represents the history of the algorithm.

Algorithm 1: SGD-MVSDE

Figure 1 below contains the flow-chart of Algorithm 1. We can also formulate a mini-batch version of the SGD-MVSDE algorithm by computing the gradient of the loss function v via averaging its value over M replications $\bigl(W^i_{m+1},\tilde W^i_{m+1},\xi^i_{m+1},\tilde\xi^i_{m+1}\bigr)^M_{i=1}$ of $(W_{m+1},\tilde W_{m+1},\xi_{m+1},\tilde\xi_{m+1})$ for a given iteration index m.

Algorithm 2: Mini-batch-SGD-MVSDE

Note that Algorithm 2 for $M=1 $ is equivalent to Algorithm 1.

Figure 1. Flow-chart of Algorithm 1.

3.3. Convergence of Algorithm 1

We recall that H is a function in $C^1(\mathbb{R}^{(n+1) K})$ and that it satisfies (17). In order to obtain a convergence result for Algorithm 1, we introduce the following assumption.

For example, the H specified in Example 5 satisfies Assumption 4. This requirement of Lipschitz regularity is essential to prove convergence of any iterative method with a sequential step (see Proposition 6 below).

Our main convergence result shows the convergence of the iterates $\textbf{a}_{m}$ in Algorithm 1 to a stationary point of G.

Theorem 1. Let Assumptions 1 and 4 be in force. Assume also that $\alpha(t,\cdot),\beta(t,\cdot)\in C^2(\mathbb{R}^d)$ , for any $t\in[0,T]$ , and $\varphi\in C^2(\mathbb{R}^K)$ , with derivatives bounded by some $R>0$ . The map G is differentiable and there exists an $\mathbb{R}^{(n+1)K}$ -valued random variable $\textbf{a}_{\infty}$ such that, in Algorithm 1, we have

\begin{equation*}\lim_{m\to + \infty} \textbf{a}_{m} = \textbf{a}_{\infty} \quad \bar{\mathrm{P}}-\textrm{almost surely}\quad \textrm{and} \quad \nabla_a G(\textbf{a}_{\infty}) = 0.\end{equation*}

To prove Theorem 1, we first prove two intermediate results. These results enable us to employ the developments in [Reference Bertsekas and Tsitsiklis6, Reference Patel and Zhang41] for the convergence of SGD algorithms in non-convex settings. The first preliminary result contains: (i) a standard consistency condition that allows us to understand the process $(\textbf{v}_{m})_{m\in\mathbb{N}}$ as a sequence of stochastic gradients of G; and (ii) a bound on the second conditional moment of the gradient noise.

Proposition 5. Under the assumptions of Theorem 1, we have

(27) \begin{equation}\mathrm{E}_{\bar{\mathrm{P}}}\big[ \textbf{v}_{m+1} - \nabla_a G( \textbf{a}_m) | \bar{\mathcal{F}}_m \big] = 0, \quad m \in \mathbb{N}_0,\end{equation}

where we set $\bar{\mathcal{F}}_0\,:\!=\,\{ \bar{\Omega}, \emptyset \}$ . Furthermore, the random variables

\begin{equation*}\mathrm{E}_{\bar{\mathrm{P}}}\big[ | \textbf{v}_{m+1} - \nabla_a G( \textbf{a}_m) |^2 \big| \bar{\mathcal{F}}_m \big] \end{equation*}

are bounded, uniformly with respect to $m\in \mathbb{N}_0$ .

The next proposition is a regularity result for G and its gradient.

We are now in the position to prove our main result in Theorem 1.

Proof of Theorem 1. Differentiability of G follows from Proposition 6. By Propositions 5 and 6, the hypothesis of [Reference Bertsekas and Tsitsiklis6, Proposition 3] and [Reference Patel and Zhang41, Theorem 1] are fulfilled. In particular, Theorem 1 in [Reference Patel and Zhang41] (see the remark on p. 5 there) implies that $\bar{\mathrm{P}}$ -almost surely either ${\textbf{a}_m}$ converges in $\mathbb{R}^{(n+1) K}$ or $|{\textbf{a}_m}|\to +\infty$ . On the other hand, by Proposition 3 in [Reference Bertsekas and Tsitsiklis6], either $G(\textbf{a}_m)\to -\infty$ or $G(\textbf{a}_m)$ converges to a finite value and $\nabla G(\textbf{a}_m)\to 0$ , $\bar{\mathrm{P}}$ -almost surely. Since

\begin{equation*}G\geq 0 \quad\textrm{and}\quad \lim_{|a|\to+\infty} G(a) = +\infty,\end{equation*}

we have that $\textbf{a}_m$ converges $\bar{\mathrm{P}}$ -almost surely to a stationary point of G.

4.1. Kuramoto–Shinomoto–Sakaguchi MV-SDE

Consider the following MV-SDE, which is related to the famous Kuramoto–Shinomoto–Sakaguchi model (see [Reference Frank24, Section 5.3.2]):

\begin{equation*} \textrm{d}X_t = \bigl(\mathrm{E}[\!\sin\!(X_t)] \cos(X_t) - \mathrm{E}[\!\cos(X_t)] \sin\!(X_t) \bigr) \textrm{d} t + \sigma \, \textrm{d} W_t, \quad X_0 = x_0.\end{equation*}

The above MV-SDE can be seen in the form of (2) with $K=3$ , $d=1$ , $q=1$ , and

\begin{align*}\varphi(x) = \bigl( \sin\!(x), \cos\!(x), 1\bigr),&&\alpha(t,x) = ( \!\cos(x), -\sin\!(x), 0)^\top,&& \beta (t,x) =( 0, 0, \sigma )^{\top}.\end{align*}

Here, $\bar{\gamma} (t) = \mathrm{E}[\varphi(X_t)]$ takes values on $\mathbb{R}^3$ . However, as $\bar{\gamma}_3(t)\equiv 1$ , we run the SGD algorithm with $K=2$ to approximate $\bar{\gamma}_1(t)=\mathrm{E}[\!\sin\!(X_t)]$ and $\bar{\gamma}_2(t) = \mathrm{E}[\!\cos\!(X_t)]$ , which are the real unknown functions. Note that $\varphi, \alpha$ , and $\beta$ in this example satisfy the assumptions required in the convergence result of Section 3.3.

We test different choices of algorithm hyperparameters: $r_0 = 1, 5, 10$ , $\rho = 0.6, 0.7, 0.8, 0.9$ , $n=3,4,5,6{, 7, 8}$ , and $M= 1, 10, 100, 1000, 10{,}000$ . The benchmark ${\bar{\gamma}}^{\textrm{MC}}$ in (28), based on the Monte Carlo particle-system approximation, is obtained with $N_{0}=10^6$ particles. For each combination of hyperparameters, we perform 1000 independent runs of the algorithm, and record the average number of iterations of the SGD algorithm required to achieve the relative error accuracy $\varepsilon_m<1\%$ .

In Tables 3–5 we report the lowest average number of iterations obtained over different $(r_0,\rho)$ combinations for a given M and n for time horizons $T=0.5, 1, 2,$ respectively. We observe that the average number of iterations decreases with the size of the mini-batch sample M, and increases with the time horizon T. By taking into account the computational time, it is clear that $M=1000$ gives the best performance. It is evident that our mini-batch-SGD algorithm converges very quickly to the Monte Carlo particle-system-based benchmark.

Table 3. Kuramoto–Shinomoto–Sakaguchi MV-SDE. Average number of iterations m (execution time in seconds) over 1000 independent runs of the algorithm to achieve relative error accuracy $\varepsilon_m<1\%$ with the best combination of $(r_0,\rho)$ for each pair (n,M). Here $T=0.5$ , $x_0=0.5$ , and $\sigma=0.5$ , and the Monte Carlo benchmark exhibited an execution time of 1.98 seconds.

Table 4. Kuramoto–Shinomoto–Sakaguchi MV-SDE. Average number of iterations m (execution time in seconds) over 1000 independent runs of the algorithm to achieve relative error accuracy $\varepsilon_m<1\%$ with the best combination of $(r_0,\rho)$ for each pair (n,M). Here $T=1$ , $x_0=0.5$ , and $\sigma=0.5$ , and the Monte Carlo benchmark exhibited an execution time of 3.60 seconds.

Table 5. Kuramoto–Shinomoto–Sakaguchi MV-SDE. Average number of iterations m (execution time in seconds) over 1000 independent runs of the algorithm to achieve relative error accuracy $\varepsilon_m<1\%$ with the best combination of $(r_0,\rho)$ for each pair (n,M). Here $T=2$ , $x_0=0.5$ , and $\sigma=0.5$ , and the Monte Carlo benchmark exhibited an execution time of 6.59 seconds.

In Figures 2–4 we plot the output curves $\mathscr{L} \textbf{a}_m=(\mathscr{L} \textbf{a}_m^1, \mathscr{L} \textbf{a}_m^2)$ of the SGD algorithm versus the benchmark curves ${\bar{\gamma}}^{\textrm{MC}}=({\bar{\gamma}}^{\textrm{MC}}_1 , {\bar{\gamma}}^{\textrm{MC}}_2),$ for all values of the polynomial degree $n=3, 4, \ldots, 8$ and for time horizons $T=0.5, 1, 2$ respectively. Note that the approximation does not deteriorate as the time horizon T increases. In particular, the slight loss of precision that can be observed towards the right end of the time interval can be explained by the choice of polynomial approximations. For instance, the gap between $(\mathscr{L} \textbf{a}_m)^1$ and ${\bar{\gamma}}^{\textrm{MC}}_1$ , which can be seen in Figure 2 in the region $t\in[0.4,0.5]$ , does not appear in the analogous plot of Figure 3 (or Figure 4) after the same number of steps. This shows that the error does not propagate over time: this is mainly due to the fact that ours is a full history-type minimization algorithm, in that it approximates the minimum in (7) over the set of the continuous curves on [0, T].

Figure 2. Kuramoto–Shinomoto–Sakaguchi MV-SDE. Comparison of the output curves $\mathscr{L} \textbf{a}_m=(\mathscr{L} \textbf{a}_m^1, \mathscr{L} \textbf{a}_m^2)$ of the SGD algorithm versus the benchmark curves ${\bar{\gamma}}^{\textrm{MC}}=({\bar{\gamma}}^{\textrm{MC}}_1 , {\bar{\gamma}}^{\textrm{MC}}_2)$ for all values of n, for time-step size $h=10^{-2}$ , $T=0.5$ , $M=1000$ , $x_0=0.5$ , and $\sigma=0.5$ .

Figure 3. Kuramoto–Shinomoto–Sakaguchi MV-SDE. Comparison of the output curves $\mathscr{L} \textbf{a}_m=(\mathscr{L} \textbf{a}_m^1, \mathscr{L} \textbf{a}_m^2)$ of the SGD algorithm versus the benchmark curves ${\bar{\gamma}}^{\textrm{MC}}=({\bar{\gamma}}^{\textrm{MC}}_1 , {\bar{\gamma}}^{\textrm{MC}}_2)$ for all values of n, for time-step size $h=10^{-2}$ , $T=1$ , $M=1000$ , $x_0=0.5$ , and $\sigma=0.5$

Finally observe that the algorithm works best for low values of the polynomial order, namely $n=3,4$ . Indeed, typically the number of iterations needed to reach the target relative error accuracy $\varepsilon_m<1\%$ (cf. Tables 3–5) increase as higher values of n are considered. For the best values of the mini-batch size, namely $M=100,1000$ , this phenomenon seems less evident. However, we can see from Figures 2–4 that the quality of the approximation deteriorates when the number of iterations does not increase with n. This can be interpreted as an overfitting phenomenon, which is expected as the benchmark curves ${\bar{\gamma}}^{\textrm{MC}}=({\bar{\gamma}}^{\textrm{MC}}_1 , {\bar{\gamma}}^{\textrm{MC}}_2)$ are relatively flat.

4.2. Polynomial drift MV-SDE

Consider the MV-SDE with polynomial drift given as follows:

\begin{equation*}\textrm{d}X_t = \left( \mathrm{E}[X_t] - X_t \mathrm{E}[X_t^2] + \delta X_t \right) \, \textrm{d}t + X_t \,\textrm{d}W_t , \quad X_0=x_0.\end{equation*}

The latter can be cast in the form of (2) with $K=3$ , $d=1$ , $q=1$ , and

\begin{equation*}\varphi(x)=(x, x^2, 1), \quad \alpha(t,x)=(1, -x, \delta x)^T, \quad \beta(t,x)=(0 , 0, x)^T.\end{equation*}

Once more, the function $\bar{\gamma} (t) = \mathrm{E}[\varphi(X_t)]$ takes values on $\mathbb{R}^3$ but $\bar{\gamma}_3(t)\equiv 1$ . Thus, we run the SGD algorithm with $K=2$ to approximate $\bar{\gamma}_1(t)=\mathrm{E}[(X_t)^2]$ and $\bar{\gamma}_2(t) = \mathrm{E}[X_t]$ , which are the real unknown functions. Also note that $\varphi, \alpha$ , and $\beta$ in this example are unbounded and thus do not satisfy the assumptions in the results of Section 3.3. In particular, the drift coefficient has super-linear growth in the measure component.

We test different choices of algorithm hyperparameters: $r_0 = 1, 5, 10$ , $\rho = 0.6, 0.7,$ and $M= 1, 10, 100, 1000$ . We set the polynomial degree $n=3,$ and test the performance of the SGD algorithm for different time horizons, $T=0.1, 0.5, 1$ . The benchmark ${\bar{\gamma}}^{\textrm{MC}}$ in (28) based on the Monte Carlo particle-system approximation is obtained with $N_{0}=10^{ 7}$ particles. For each combination of the hyperparameters, we perform 1000 independent runs of the algorithm, and record the average number of iterations of the SGD algorithm required to achieve the relative error accuracy $\varepsilon_m<1\%$ .

In Table 6 we report the lowest average number of iterations over different $(r_0,\rho)$ combinations for a given M and T. Once more, as the value of the mini-batch increases the number of steps to reach convergence decreases. Compared to the model tested in Section 4.1, increasing the time horizon T has a stronger impact on the convergence, which is due to the quadratic growth of the coefficient. We stress that for this model even the MC particle method requires an increasing number of particles to obtain reliable benchmark curves for increasing values of T.

Table 6. Polynomial drift MV-SDE. Average number of iterations (execution time in seconds) over 1000 independent runs of the algorithm to achieve relative error accuracy $\varepsilon_m<1\%$ with the best combination of $(r_0,\rho)$ for each pair (T,M). Here $x_0=1$ , $\delta=0.8$ and the Monte Carlo benchmark exhibited an execution time of 2.57 seconds for $T=0.1$ , 9.78 seconds for $T=0.5$ , and 19.53 seconds for $T=1$ .

Figure 4. Kuramoto–Shinomoto–Sakaguchi MV-SDE. Comparison of the output curves $\mathscr{L} \textbf{a}_m=(\mathscr{L} \textbf{a}_m^1, \mathscr{L} \textbf{a}_m^2)$ of the SGD algorithm versus the benchmark curves ${\bar{\gamma}}^{\textrm{MC}}=({\bar{\gamma}}^{\textrm{MC}}_1 , {\bar{\gamma}}^{\textrm{MC}}_2)$ for all values of n, for time-step size $h=10^{-2}$ , $T=2$ , $M=1000$ , $x_0=0.5$ , and $\sigma=0.5$

In Figure 5 we plot the output curves $\mathscr{L} \textbf{a}_m=(\mathscr{L} \textbf{a}_m^1, \mathscr{L} \textbf{a}_m^2)$ of the SGD algorithm versus the benchmark curves ${\bar{\gamma}}^{\textrm{MC}}=({\bar{\gamma}}^{\textrm{MC}}_1 , {\bar{\gamma}}^{\textrm{MC}}_2)$ for all time horizons $T=0.1$ , $T=0.5$ , and $T=1$ . In the plots, particularly for $T=0.5$ , notice the considerable loss of accuracy in our approximation $\mathscr{L} \textbf{a}_m^2$ of $\bar{\gamma}_2(t) = \mathrm{E}[X_t]$ near the right boundary of the interval (similar effects were visible in the plots relative to the Kuramoto model). After running additional tests, we do not believe that this effect is due to the polynomial projection, similar to Runge-type phenomena in polynomial interpolation over compact sets, but we claim that it is rather due to the norm used to define the cost function $F^2$ (cf. (8)), in combination with the chosen threshold level we set for $\varepsilon_m$ in (29). Indeed, lowering this threshold from to $10^{-2}$ to $10^{-3}$ reduces considerably the error, also near the boundaries of the interval. However, one could look for strategies to distribute the error in a more uniform way across the interval, and obtain a better qualitative result while keeping the same error threshold. One possibility of improvement is to first employ the alternative objective function presented in Remark 5, i.e. using a weight function to force the SGD algorithm to be more accurate towards the end of the interval. One can also use a similar weight function to compute a weighted $L^2({[0,T]})$ norm for error $\varepsilon_m$ in the stopping criterion. Another possibility is to experiment with different choices of basis functions listed in Examples 3 and 4 and choose the one with the best error performance.

Figure 5. Polynomial drift MV-SDE. Comparison of the output curves $\mathscr{L} \textbf{a}_m=(\mathscr{L} \textbf{a}_m^1, \mathscr{L} \textbf{a}_m^2)$ of the SGD algorithm versus the benchmark curves ${\bar{\gamma}}^{\textrm{MC}}=({\bar{\gamma}}^{\textrm{MC}}_1 , {\bar{\gamma}}^{\textrm{MC}}_2)$ for all values of T, for time-step size $h=10^{-2}$ , $n=3$ , $M=1000$ , $x_0=1$ , and $\delta=0.8$ .

4.3. Convolution-type MV-SDE

We consider the following convolution-type MV-SDE:

(30) \begin{equation}\textrm{d} X_t = \left. \mathrm{E}\left[\exp\left(-{(X_t- x)^2}/{2}\right)\right]\right|_{x=X_t} \,\textrm{d} t + \sigma\, \textrm{d} W_t, \quad X_0 = \xi \sim \mathcal{N}_{(0,1)}.\end{equation}

As shown in [Reference Belomestny and Schoenmakers5], the equation above can be approximated by employing a projection technique based on generalized Fourier series, which we sketch out here for the reader’s convenience. Noticing that the drift coefficient in (30) can be written as

\begin{equation*}\mathrm{E}[b(x,X_t)]|_{x=X_t}, \quad b(x,y) = \exp\left(-{(y- x)^2}/{2}\right),\end{equation*}

we can expand $b(\cdot,x)$ in a suitable weighted $L^2$ space as follows. Setting the measure $\mu(\textrm{d}x) = e^{-x^2/2} \, \textrm{d}x$ on the Borel sets of $\mathbb{R}$ , we consider the following orthonormal basis of $L^2(\mathbb{R},\mu)$ of the normalized Hermite polynomials

\begin{equation*}\overline{H}_k(x) = c_k (\!-1)^k \mathrm{e}^{x^2} \frac{\textrm{d}^k}{\textrm{d} x^k} \bigl( \mathrm{e}^{-x^2}\bigr), \quad c_k = \bigl(2^k k! \sqrt{\pi} \bigr)^{-1/2}, \quad k\in\mathbb{N}_0.\end{equation*}

Therefore, for $K\in\mathbb{N}$ sufficiently large, $b(x,\cdot)$ can be approximated in $L^2(\mathbb{R},\mu)$ as

(31) \begin{equation}b(x,y) \approx \sum_{k=0}^{K} \alpha_k(x) \varphi_k(y),\end{equation}

where

\begin{equation*}\varphi_k(x) = \overline{H}_k(x) \mathrm{e}^{-x^2/2}, quad \alpha_k(x)= \int_{\mathbb{R}} b(x,y) \varphi_k(y) \, \textrm{d}y = \pi^{1/4}\Bigl( \frac{1}{2}\Bigr)^{k/2} \frac{x^k}{\sqrt{k!}} \mathrm{e}^{-x^2/4}.\end{equation*}

Therefore, the MV-SDE (4.3) can be approximated by the MV-SDE with separable coefficients

\begin{equation*}\textrm{d} X^{(K)}_t = \sum^K_{k=0}\mathbb {E}\big[\varphi_k(X^{(K)}_t)\big] \alpha_k\big(X^{(K)}_t\big) \, \textrm{d} t + \sigma \,\textrm{d} W_t,\quad X_0 = \xi = \mathcal{N}_{(0,1)},\end{equation*}

which can be clearly cast in the form of (2). For a given $K\in\mathbb{N}$ , we run the SGD algorithm in order to approximate $\bar{\gamma}^{(K)}_k(t) = \mathbb {E}[\varphi_k(X^{(K)}_t)]$ for $k=0,\ldots,K$ .

Figure 6. Monte Carlo method densities $\tilde{w}^{(K),\textrm{MC}}_T(x)$ over the interval $[\!-3,4]$ , for $T=1$ and $K = 3, 5, 10, 20$ . The benchmark vector ${\bar{\gamma}}^{\textrm{MC}}(T)$ was computed with $N_{0}=10^7$ particles. Here $X_0 \sim \mathcal{N}_{(0,1)}$ and $\sigma = 0.1$ . In the model (30) we had $X_0 \sim \mathcal{N}_{(0,1)}$ and $\sigma = 0.1$ .

For this model, we test the accuracy of the SGD algorithm by checking the quality of the density approximation of the solution $X_T$ to Section 4.3 that can be derived from the $\bar{\gamma}^{(K)}_k$ curves as follows (see also [Reference Belomestny and Schoenmakers5]). For a fixed $T>0$ denote by $w_T(x)$ and $w^{(K)}_T(x)$ the densities of $X_T$ and $X^{(K)}_T$ , respectively. Following the same argument that led to (31), for $K\in\mathbb{N}$ sufficiently large, $w_T(x)$ can be approximated in $L^2(\mathbb{R},\mu)$ as

\begin{equation*}w_T(x) \approx w^{(K)}_T(x) \approx \sum_{k=0}^{K} \int_{\mathbb{R}} \varphi_k(y) w^{(K)}_T(y) \, \textrm{d}y\, \varphi_k(x) = \sum_{k=0}^{K} \bar{\gamma}^{(K)}_k(T) \varphi_k(x) =: \tilde{w}^{(K)}_T(x).\end{equation*}

Denote by $\tilde{w}^{(K),\textrm{SGD}}_T$ and $\tilde{w}^{(K),\textrm{MC}}_T$ , the approximations of $\tilde{w}^{(K)}_T$ obtained by replacing the vector $\bar{\gamma}^{(K)}(T)$ with the output of the SGD algorithm, $\mathscr{L} \textbf{a}_m(T)$ , and with the particle-system benchmark, ${\bar{\gamma}}^{\textrm{MC}}(T)$ , respectively.

Figure 7. Comparison of $\tilde{w}^{(K),\textrm{SGD}}_T(x)$ and $\tilde{w}^{(K),\textrm{MC}}_T(x)$ over the interval $[\!-3,4]$ , for $T=1$ and $K = 10$ . The benchmark vector ${\bar{\gamma}}^{\textrm{MC}}(T)$ was computed with $N_{0}=10^7$ particles. The SGD output $\mathscr{L} \textbf{a}_m(T)$ was obtained after $m=172$ iterations (35 seconds of computation time), with parameters $n = 3$ , $M = 100$ , $r_0 = 5$ , and $\rho = 0.9$ . In the model (30) we had $X_0 \sim \mathcal{N}_{(0,1)}$ and $\sigma = 0.1$ .

In our tests we set $\sigma = 0.1$ and $T = 1$ as in [Reference Belomestny and Schoenmakers5]. The number of particles to derive the benchmark ${\bar{\gamma}}^{\textrm{MC}}(T)$ is set as $N_{0}=10^7$ . Note that in [Reference Belomestny and Schoenmakers5] the authors used $N_{0}=5\times 10^2$ , which is why our density plots differ considerably from theirs.

In Figure 6 we plot the approximate density $\tilde{w}^{(K),\textrm{MC}}_T(x)$ for $K = 3, 5, 10, 20$ , with x ranging on the interval $[\!-3,4]$ . From this we conclude that the value $K=10$ is high enough in order for $\tilde{w}^{(K)}$ to accurately approximate ${w}_T$ over the interval $[\!-3,4]$ . In Figure 7 we test the accuracy of the SGD algorithm by plotting $\tilde{w}^{(K),\textrm{SGD}}_T(x)$ versus $\tilde{w}^{(K),\textrm{MC}}_T(x)$ , with $K=10$ , over the same interval. The parameters of the SGD algorithm were chosen as $n = 3$ , $M = 100$ , $r_0 = 5$ , $\rho = 0.9$ , and the algorithm halted after $m=172$ iterations (35 seconds of computation time) when the relative error reached $\varepsilon_m<1\%$ .