1 Introduction↩︎

Stochastic partial differential equations (SPDEs) have emerged as essential modeling frameworks in a wide range of fields, including finance, environmental sciences, biology, engineering, and spatial statistics, where phenomena exhibit both spatial-temporal dynamics and inherent stochasticity [1], [2].

The numerical analysis of parabolic SPDEs has been the subject of extensive investigation in the literature, with particular attention devoted to the construction and analysis of high-order time discretization schemes [1], [3]. Notably, in [4] the authors proposed a Milstein finite difference scheme for equations in one spatial dimension, while in [5] the authors extended these studies to two spatial dimensions.

In this work, we consider a different class of SPDEs, specifically degenerate Langevin-type SPDEs in two spatial dimensions, which are formulated as follows: \[\label{SPDE} d_\mathbf{Y} u_t(x,v) = \frac{a_t(x,v)}{2} \partial_{vv} u_t(x,v) dt + \sigma_t(x,v) \partial_v u_t(x,v) dW_t, \qquad \mathbf{Y} := \partial_t + v \partial_x,\tag{1}\]

where \(W\) denotes a Wiener process defined on a complete probability space \((\Omega, \mathscr{F}, P)\) endowed with a filtration \((\mathscr{F}_t)_{t \geqslant 0}\) satisfying the usual conditions.

This class of equations arises naturally in mathematical finance, where 1 describes the price dynamics of path-dependent financial contracts, such as Asian options [6], driven by a common stochastic factor \(W\) representing the market uncertainty.

Since, even in the deterministic setting (i.e., \(\sigma \equiv 0\)), analytic solutions for this type of equation are generally not available, various numerical schemes have been proposed over the years. For example, explicit and implicit schemes for a bounded value problem related to the Kolmogorov equation are studied in [7], while the authors in [8] provided both theoretical results for Kolmogorov PDEs and associated numerical methods.

More recently, in the context of degenerate SPDEs, the authors in [9] demonstrated how the Itô-stochastic Magnus expansion can be employed to efficiently solve such equations in two spatial dimensions numerically.

In the present work, we will focus on the constant coefficient case, i.e., when \(a\) and \(\sigma\) are constant, where an analytical solution to 1 is available (see [10], [11]). The availability of a closed-form solution in this case helps us to check the validity of our numerical scheme and its convergence rate, whereas the scheme itself is (hopefully) more widely applicable.

First of all, let us clarify the concept of solution of the SPDE that we will consider. Since the symbol \(d_\mathbf{Y}\) indicates that the equation is solved in the Itô (or strong) sense, a solution to 1 is a continuous process \(u_t = u_t(x,v)\) that is twice differentiable in \(v\) and such that \[\label{SPDE95gamma} u_t \left ( \gamma^B_t (x,v) \right ) = u_0(x,v) + \frac{1}{2} \int_0^t{( a_s \partial_{vv} u_s ) \left ( \gamma^B_s (x,v) \right ) ds} + \int_0^t{(\sigma_s \partial_v u_s) \left ( \gamma^B_s (x,v) \right ) dW_s},\tag{2}\] where \(t \mapsto \gamma^B_t (x,v)\) denotes the integral curve of the advection field \(v \partial_x\) starting from \((x,v)\), that is \[\gamma^B_t (x,v) = e^{tB}(x,v) = (x+tv,v), \qquad B = \left ( \begin{matrix} 0 \quad 1 \\ 0 \quad 0 \end{matrix} \right ).\]

Let \(\overline{W}, W\) be independent real Brownian motions, \(a>0\) and \(\sigma \in [0,\sqrt{a}]\). The Langevin model is defined in terms of the system of SDEs \[\label{langevin95xv} \begin{cases} dX_t = V_t dt \\ dV_t = \sqrt{a - \sigma^2} d\overline{W}_t - \sigma dW_t. \end{cases}\tag{3}\] We can interpret these equations as those that govern the displacement of a particle in a phase space, where the velocity of the particle, \(V_t\), is only partially observed through the process \(W\), where \(\overline{W}\) represents the noise in the observation.

It can be shown that the SPDE 1 is the forward Kolmogorov equation of the SDE 3 conditioned to the Brownian observation given by \(\mathscr{F}_t^W = \sigma(W_s, s \leqslant t)\), see [10], [11], for further details.

2 Numerical method for the Langevin SPDE↩︎

As implementing a numerical scheme to work directly on 1 can be quite challenging, then we follow the alternative approach proposed in [7] and [8]. For this purpose, recalling the definition of a solution to 1 , for \(\delta> 0\), we have \[\label{SPDE95t95d} u_{t+\delta} \left ( \gamma^B_\delta(x,v) \right ) = u_t(x,v) + \frac{1}{2} \int_t^{t+\delta}{( a_s \partial_{vv} u_s ) \left ( \gamma^B_{s-t} (x,v) \right ) ds} + \int_t^{t+\delta}{(\sigma_s \partial_v u_s) \left ( \gamma^B_{s-t} (x,v) \right ) dW_s}.\tag{4}\]

If we now take the drift term in implicit form and the diffusion term in explicit form, we get a semi-implicit scheme. More precisely, we approximate \(u\) in the first integral of 4 as \(u_s(x,v) \approx u_{t+\delta}(x,v)\) and thus \((a_s \partial_{vv}u_s) \left ( \gamma^B_{s-t} (x,v) \right ) \approx (a_{t+\delta} \partial_{vv} u_{t+\delta}) \left ( \gamma^B_\delta(x,v) \right ) = (a_{t+\delta} \partial_{vv} u_{t+\delta}) (x+\delta v,v)\). Then we can use a zero order approximation for the diffusion term (leading to a semi-implicit Euler-Maruyama scheme) or a first order one (leading to a semi-implicit Milstein scheme).

Now, let us split the time interval \([0,T]\) into \(0=t_0<t_1<\dots<t_N=T\), where \(t_n = n \delta\), and let us write \(U^n\) to indicate the solution at time \(t_n\) on a space grid.

In the semi-implicit Euler-Maruyama scheme, we approximate \(u\) in the second integral in 4 as \(u_s(x,v) \approx u_{t}(x,v)\) and thus \((\sigma_s \partial_{v}u_s) \left ( \gamma^B_{s-t} (x,v) \right ) \approx (\sigma_t \partial_{v}u_t) \left ( \gamma^B_0 (x,v) \right ) = (\sigma_t \partial_{v}u_t) (x,v)\). Therefore, the SPDE could be approximated by \[S_{\delta}U^{n+1} = U^n + \frac{\delta}{2} S_{\delta}K_1 U^{n+1} + \sqrt{\delta} Z_{n} K_2 U^n,\] where \(K_1\) and \(K_2\) are the discretizations of the operators \(a_t \partial_{vv}\) and \(\sigma_t \partial_v\) respectively, \(Z_{n} \sim \mathcal{N}_{0,1}\) are independent normal random variables, for every \(n = 0,\dots,N-1\). \(S_{\delta}\) is the shift matrix such that \(u_t(x+\delta v,v) = S_{\delta}u_t(x,v)\), we will be more precise in the definition of the operator \(S_{\delta}\) later on.

In the semi-implicit Milstein scheme, we approximate \(u\) in the second integral as \(u_s(x,v) \approx u_{t}(x,v) + \sigma_t \partial_v u_t (x,v) (W_s-W_t)\) and thus \((\sigma_s \partial_v u_s) \left ( \gamma^B_{s-t} (x,v) \right ) \approx (\sigma_t \partial_v u_t + \sigma_t^2 \partial_{vv} u_t) (x,v) (W_s-W_t)\). Recalling that \(\int_t^{t+\delta}{(W_s-W_t)dW_s} = \frac{1}{2} \left ( (\Delta W^{n+1})^2 - \delta\right ) = \frac{\delta}{2} (Z_{n+1}^2 - 1)\), the SPDE could be approximated by \[\label{Milstein95scheme} S_{\delta}U^{n+1} = U^n + \frac{\delta}{2} S_{\delta}K_1 U^{n+1} + \sqrt{\delta} Z_{n} K_2 U^n + \delta\frac{1}{2} (Z_{n}^2 - 1) K_3 U^n,\tag{5}\] where \(K_3\) is the discretization of the operator \(\sigma_t^2 \partial_{vv} u_t\). This scheme can be written as \[(I - \frac{\delta}{2} K_1)U^{n+1} = S_{-\delta}(I + \sqrt{\delta} Z_{n+1} K_2 + \delta\frac{1}{2} (Z_{n+1}^2 - 1) K_3) U^n,\]

Let us suppose that we are working on \(t \in [0,T], x \in [0,L_1], v \in [0,L_2]\), we need a grid \(G\) such that \((t_n,x_j+\delta v_k,v_k) = (t_n,x_{\bar{j}},v_k) \in G\), \(\forall (t_n,x_j,v_k) \in G\) with \(x_j+\delta v_k \leqslant L_1\).

A possible solution to this is again to follow [7], i.e. fixing \(\delta\) and \(h_v\), and setting \(h_x := \delta h_v\). In this setting, now we have that \(t_n = n\delta, x_j = jh_x, v_k = kh_v\) and thus \((t_n,x_j+\delta v_k,v_k) = (n\delta,jh_x+\delta kh_v,kh_v) = (n\delta,(j+k)h_x,kh_v) = (t_n,x_{j+k},v_k)\), for every \(n=0,\dots,N\), \(j=1,\dots,N_x\), \(k=1,\dots,N_v\) such that \(k+j \leqslant N_x\).

Therefore, we can define the shift matrix as follows \[\label{shift95matrix} S_{\delta}:= \begin{bmatrix} S & & & \\ & S^2 & & \\ & & \ddots & \\ & & & S^{N_v-1} \\ & & & & \mathrm{0}_{N_v} \end{bmatrix} \in \mathbb{R}^{N_x N_v \times N_x N_v},\tag{6}\] which can be rewritten as \(S_{\delta}= \sum_{j=1}^{N_v}{E^{(j)} \otimes S^j}\), where \(E^{(j)} \in \mathbb{R}^{N_v \times N_v}\) are zero matrices with only \(E^{(j)}_{j,j}=1\) and \[S := \begin{bmatrix} 0 & 1 & & & \\ & 0 & 1 & & \\ & & \ddots & \ddots & \\ & & & 0 & 1 \\ & & & & 0 \end{bmatrix} \in \mathbb{R}^{N_x \times N_x}.\] Notice that the matrix \(S_{-\delta}\) is the transpose of \(S_{\delta}\).

In general \(N_v \ll N_x\), since \(h_x = \delta h_v\) and we are working in an almost square spatial domain. Moreover, if the spatial domain is large enough, we know from the theoretical results that the SPDE solution is negligible close to the boundary, enabling us to neglect its behavior in that region.

3 Fourier stability and convergence analysis↩︎

3.1 Fourier analysis of mean-square stability↩︎

In order to perform the Fourier analysis of mean-square stability, let us first recall the definition of the Fourier transform and inverse Fourier transform pair \[\begin{align} \hat{u}_t(\xi,\eta) &= \int_{\mathbb{R}^2}{u_t(x,v) e^{-i \langle (x,v),(\xi,\eta) \rangle} dx dv}, \\ u_t(x,v) &= \frac{1}{4 \pi^2}\int_{\mathbb{R}^2}{\hat{u}_t(\xi,\eta) e^{i \langle (x,v),(\xi,\eta) \rangle} d\xi d\eta}. \end{align}\] Now, we compute the Fourier transform of the composition \(u_t \circ \gamma_\delta^B\) \[\begin{align} \mathcal{F}(u_t \circ \gamma_t^B)(\xi,\eta) &= \int_{\mathbb{R}^2}{(u_t \circ \gamma_t^B)(x,v) e^{-i \langle (x,v),(\xi,\eta) \rangle} dx dv} \\ &= \int_{\mathbb{R}^2}{u_t(x,v) e^{-i \langle (x - t v,v),(\xi,\eta) \rangle} dx dv} \\ &= \int_{\mathbb{R}^2}{u_t(x,v) e^{-i \langle (x,v),(\xi,\eta - t \xi) \rangle} dx dv} = \hat{u}_t(\xi,\eta - t \xi). \end{align}\] Analogously, we obtain \[\begin{align} \mathcal{F}((a_t \partial_{vv} u_t) \circ \gamma_t^B)(\xi,\eta) &= \mathcal{F}(a_t \partial_{vv} u_t) (\xi,\eta-t\xi), \\ \mathcal{F}((\sigma_t \partial_v u_t) \circ \gamma_t^B)(\xi,\eta) &= \mathcal{F}(\sigma_t \partial_v u_t) (\xi,\eta-t\xi). \end{align}\] In the case of constant coefficients satisfying \(a > 0\) and \(\sigma \in [0, \sqrt{a})\), we get \[\begin{align} \mathcal{F}((a \partial_{vv} u_t) \circ \gamma_t^B)(\xi,\eta) &= -(\eta-t\xi)^2 a \hat{u}_t (\xi,\eta-t\xi), \\ \mathcal{F}((\sigma \partial_v u_t) \circ \gamma_t^B)(\xi,\eta) &= i (\eta-t\xi) \sigma \hat{u}_t (\xi,\eta-t\xi). \end{align}\]

By using the previous computations, the Fourier transform of the SPDE 2 is given by \[d\hat{u}_t(\xi,\eta-t\xi) = \left ( -\frac{1}{2} a (\eta-t\xi)^2 dt + i \sigma (\eta-t\xi) dW_t \right ) \hat{u}_t (\xi,\eta-t\xi),\] so that its solution can be written as follows \[\hat{u}_t (\xi,\eta-t\xi) = X_t(\xi,\eta) e^{-i \langle (x_0,v_0),(\xi,\eta) \rangle},\] where \[\label{X95t} X_t(\xi,\eta) := \exp \left ( \frac{1}{2} (\sigma^2-a) \left ( \eta^2 t - \xi \eta t^2 + \frac{\xi^2}{3} t^3 \right ) + i \sigma (\eta - t \xi) W_t \right ).\tag{7}\] Thus, by a change of variable, we get \[\label{Fourier95sol} \hat{u}_t (\xi,\eta) = X_t(\xi,\eta + t \xi) e^{-i \langle (x_0,v_0),(\xi,\eta + t \xi) \rangle}.\tag{8}\] From 8 , we have \(\mathbb{E}[|X_t|^2]\xrightarrow[t \to +\infty]{} 0\), which implies that \(\mathbb{E}[\| \hat{u}_t\|_{L^2}^2] \xrightarrow[t \to +\infty]{} 0\) and, by isometry, we also have \(\mathbb{E}[\|u_t\|_{L^2}^2] \xrightarrow[t \to +\infty]{} 0\).

For the numerical solution, we use a discrete-continuous Fourier decomposition. If we denote \(U^n_{j,k} = u(n\delta,j h_x, k h_v)\), then \[U^0_{j,k} = \frac{1}{4\pi^2} \int_{-\frac{\pi}{h_v}}^{\frac{\pi}{h_v}}{\int_{-\frac{\pi}{h_x}}^{\frac{\pi}{h_x}}{\hat{U}^0(w,z) e^{i \left ( (j-j_0)h_x w + (k-k_0)h_v z \right )} dw dz}},\] where \(j_0 = x_0/h_x, k_0 = v_0/h_v\), and \[\hat{U}^0(\xi,\eta) := h_x h_v \sum_{j \in \mathbb{Z}} \sum_{k \in \mathbb{Z}} U^0_{j,k} e^{-i \left ( (j-j_0) h_x \xi + (k-k_0) h_v \eta \right )}.\] If we set \(U^0_{j,k} = h_x^{-1} h_v^{-1} \delta_{(j_0,k_0)}\), we have \(\hat{U}^0(w,z) = 1\) for all \((w,z) \in \mathbb{R}^2\).
Similarly, for the \(n\)-th step \[U^n_{j,k} = \frac{1}{4\pi^2} \int_{-\frac{\pi}{h_v}}^{\frac{\pi}{h_v}}{\int_{-\frac{\pi}{h_x}}^{\frac{\pi}{h_x}}{\hat{U}^n(\xi,\eta) e^{i \left ( (j-j_0) h_x \xi + (k-k_0) h_v \eta \right )} d\xi d\eta}}.\]

By analogy to the theoretical solution 8 , we make the ansatz \[\hat{U}^n(\xi,\eta) = X^n(\xi,\eta + n\delta\xi) \hat{U}^0(\xi,\eta) = X^n(\xi,\eta + n\delta\xi),\] since \(\hat{U}^0 \equiv 1\), where \(X^n(\xi,\eta)\) can be considered the numerical approximation of \(X(n\delta)\).

Analogously to the continuous case, we say that a numerical solution is asymptotically mean-square stable if for any \((\xi,\eta) \in [-\frac{\pi}{h_x},\frac{\pi}{h_x}] \times [-\frac{\pi}{h_v},\frac{\pi}{h_v}]\), we have \[\label{mean95square95stability} \lim_{n \to + \infty}{\mathbb{E}[|X^n(\xi,\eta)|^2]} = 0.\tag{9}\]

First, for convenience, let us rewrite the scheme 5 in the form \[(I - \frac{\delta}{2} K_1)U^{n+1} = S_{-\delta}(I + \sqrt{\delta} Z_{n+1} K_2 + \delta\frac{1}{2} (Z_{n+1}^2 - 1) K_3) U^n,\] where \(K_1 = a D_{vv}, K_2 = \sigma D_v\) and \(K_3 = \sigma^2 D_v^2\).

For simplicity, next we assume that \(j_0 = k_0 = 0\), and we compute \[\begin{align} \mathcal{F}(S_{-\delta}D_v U)(\xi,\eta) &= h_x h_v \sum_{j \in \mathbb{Z}} \sum_{k \in \mathbb{Z}} (S_{-\delta}\frac{1}{2 h_v} \left ( U_{j,k+1} - U_{j,k-1} \right ) ) e^{-i \left ( j h_x \xi + k h_v \eta \right )} \\ &= h_x h_v \sum_{j \in \mathbb{Z}} \sum_{k \in \mathbb{Z}} \frac{1}{2 h_v} \left ( U_{j-k-1,k+1} - U_{j-k+1,k-1} \right ) e^{-i \left ( j h_x \xi + k h_v \eta \right )} \\ &= h_x h_v \sum_{j \in \mathbb{Z}} \sum_{k \in \mathbb{Z}} U_{j,k} e^{-i \left ( (j+k) h_x \xi + k h_v \eta \right )} \frac{1}{2 h_v} \left ( e^{-i \eta h_v} - e^{i \eta h_v} \right ) \\ &= h_x h_v \sum_{j \in \mathbb{Z}} \sum_{k \in \mathbb{Z}} U_{j,k} e^{-i \left ( j h_x \xi + k h_v (\eta + \delta\xi) \right )} \\ &= i q_1(\eta)\hat{U}(\xi,\eta + \delta\xi), \qquad q_1(\eta) := \frac{\sin(\eta h_v)}{h_v}, \end{align}\] where we have applied the identity \(e^{-i \eta h_v} - e^{i \eta h_v} = i 2 \sin(\eta h_v)\).

Analogously, we obtain \[\mathcal{F}(S_{-\delta}D_v^2 U)(\xi,\eta) = - q_1^2(\eta) \hat{U}(\xi,\eta + \delta\xi).\] Moreover, we consider \[\begin{align} D_{vv} U_{j,k} &= \frac{1}{h_v^2} \left ( U_{j,k+1} - 2 U_{j,k} + U_{j,k-1} \right ) \\ &= \frac{1}{4\pi^2} \int_{-\frac{\pi}{h_v}}^{\frac{\pi}{h_v}}{\int_{-\frac{\pi}{h_x}}^{\frac{\pi}{h_x}}{\hat{U}(\xi,\eta) e^{i \left ( j\xi h_x + k\eta h_v \right )} \frac{1}{h_v^2} \left ( e^{i \eta h_v} - 2 + e^{-i \eta h_v} \right ) d\xi d\eta}}; \\ \end{align}\] so that we get \[\begin{align} \mathcal{F} (D_{vv} U)(\xi,\eta) &= - q_2(\eta) \hat{U}(\xi,\eta), \qquad q_2(\eta) := \frac{4\sin^2(\frac{\eta h_v}{2})}{h_v^2}; \\ \end{align}\] where we have used the identity \(e^{i \eta h_v} - 2 + e^{-i \eta h_v} = - 4\sin^2(\frac{\eta h_v}{2})\). Therefore, applyng the Fourier transform to 5 we obtain \[\left ( 1 + \frac{\delta}{2} a q_2(\eta) \right ) X(\xi,\eta + (n+1) \delta\xi)^{n+1} = \left ( 1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta) \right ) X(\xi,\eta + (n+1) \delta\xi)^n,\] from which we obtain the expression\[\label{C95n} X(\xi,\eta + (n+1) \delta\xi)^{n+1} = C_n X(\xi,\eta + (n+1) \delta\xi)^n, \qquad C_n := \frac{1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta)}{1 + \frac{\delta}{2} a q_2(\eta)},\tag{10}\] Recursively using the previous expression and that \(X^0 \equiv 1\), we get \[\label{X95N} X(\xi,\eta + (n+1) \delta\xi)^{n+1} = \prod_{l = 0}^n C_l X(\xi,\eta + (n+1) \delta\xi)^0 = \prod_{l = 0}^n C_l.\tag{11}\]

Next, since the \(C_n\) are independent, a sufficient condition for stability is to prove that \[\label{suff95cond95stability} \mathbb{E}\left [ \left | \frac{1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta)}{1 + \frac{\delta}{2} a q_2(\eta)} \right |^2 \right ] < 1,\tag{12}\] for any \((\xi,\eta) \in [-\frac{\pi}{h_x},\frac{\pi}{h_x}] \times [-\frac{\pi}{h_v},\frac{\pi}{h_v}] \backslash \{ (0,0) \}\).

In order to prove 12 , note that \[| 1 + \frac{\delta}{2} a q_2(\eta) |^2 = 1 + \delta a q_2(\eta) + \frac{\delta^2}{4} a^2 q_2^2(\eta),\] and \[\begin{align} &\left | 1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta) \right |^2 \\ &= \left ( 1 - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta) \right )^2 + \left ( \sqrt{\delta} \sigma Z_n q_1(\eta) \right )^2 \\ &= 1 - \delta\sigma^2 \left ( Z_n^2 - 1 \right ) q_1^2(\eta) + \delta^2 \frac{\sigma^4}{4} \left ( Z_n^2 - 1 \right )^2 q_1^4(\eta) + \delta\sigma^2 Z_n^2 q_1^2(\eta). \end{align}\] Next, since \(\mathbb{E}\left [ Z_n^2 \right ] = 1\) and \(\mathbb{E}\left [ Z_n^4 \right ] = 3\), by taking the expectation, we have \[\mathbb{E}\left [ \left | 1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta) \right |^2 \right ] = 1 + \delta^2 \frac{\sigma^4}{2} q_1^4(\eta) + \delta\sigma^2 q_1^2(\eta).\] Moreover, since \(\sigma^2 < a\), \(2 \sigma^4 < a^2\), then the following inequalities hold

• \(\sigma^2 q_1^2(\eta) < a q_2(\eta)\);

• \(\sigma^4 q_1^4(\eta) < \frac{a^2}{2} q_2^2(\eta)\).

and we end the proof.

3.2 Fourier analysis of convergence↩︎

In order to prove the convergence Theorem, we first need to define the following sets \[\Omega_{low} := \{ (\xi,\eta) \in \mathbb{R}^2 : |\eta| \leqslant{\mathord{{\rm min}}}\{h_v^{-2p},\delta^{-p}\} \}\] and \[\Omega_{high} := \Omega_{low}^c \cap [-\pi h_x^{-1},\pi h_x^{-1}] \times [-\pi h_v^{-1},\pi h_v^{-1}],\] for some \(0<p<\frac{1}{4}\).

From 7 , we have \[X_{t_{n+1}}(\xi,\eta + t_{n+1} \xi) = X_{t_n}(\xi,\eta + t_{n+1} \xi) \exp \left ( \frac{1}{2} (\sigma^2-a) \left ( \eta^2 \delta+ \xi \eta \delta^2 + \frac{\xi^2}{3} \delta^3 \right ) + i \sigma \eta \sqrt{\delta} Z_n - i \sigma \xi \delta W_{t_n} \right ).\] Now, if \(X_n\) is the numerical approximation of \(X_{\delta n}(\xi,\eta + (n+1) \delta\xi)\) we have \[X_{n+1} = C_n X_n,\] where \[C_n = \exp \left ( \frac{1}{2} (\sigma^2-a) \left ( \eta^2 \delta+ \xi \eta \delta^2 + \frac{\xi^2}{3} \delta^3 \right ) + i \sigma \eta \sqrt{\delta} Z_n - i \sigma \xi \delta W_{t_n} + e_n \right ),\] and \(e_n\) is the logarithmic error between the numerical solution and the exact solution introduced during \([n\delta,(n+1)\delta]\). Aggregating over \(N\) time steps, at \(t_N = \delta N = T\), \[X^N(\xi,\eta + N \delta\xi) = X(T,\xi,\eta+T\xi) \exp{\left ( \sum_{n=0}^{N-1}{e_n} \right ) },\] where \(X(T)\) is the exact solution at time \(T\).

Thus, we have \[e_n = \log{C95n} - \frac{1}{2} (\sigma^2-a) \left ( \eta^2 \delta+ \xi \eta \delta^2 + \frac{\xi^2}{3} \delta^3 \right ) - i \sigma \eta \sqrt{\delta} Z_n + i \sigma \xi \delta W_{t_n},\] and \[\sum_{n=0}^{N-1}{e_n} = \sum_{n=0}^{N-1}{\log{C95n}} - \frac{1}{2} (\sigma^2-a) \left ( \eta^2 + \xi \eta \delta+ \frac{\xi^2}{3} \delta^2 \right ) T - i \sigma \eta \sqrt{\delta} \sum_{n=0}^{N-1}{Z_n} + i \sigma \xi \delta\sum_{n=0}^{N-1}{W_{t_n}}.\]

By 10 , we have \[C_n = \frac{1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta)}{1 + \frac{\delta}{2} a q_2(\eta)},\]

By Taylor expansion \[\begin{align} \log{C95n} &= \log{\left ( 1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta) \right )} - \log{\left (1 + \frac{\delta}{2} a q_2(\eta) \right )} \\ &= i \sqrt{\delta} \sigma q_1(\eta) Z_n - \delta\frac{\sigma^2}{2} q_1^2(\eta) Z_n^2 + \delta\frac{\sigma^2}{2} q_1^2(\eta) + \delta\frac{\sigma^2}{2} q_1^2(\eta) Z_n^2 - \frac{\delta}{2} a q_2(\eta) \\ &+ O \left ( \delta\sqrt{\delta} \right ) i \phi_1(Z_n) + O \left ( \delta^2 \right ) \phi_2(Z_n) \\ &= i \sqrt{\delta} \sigma q_1(\eta) Z_n + \frac{\delta}{2} \left ( \sigma^2 q_1^2(\eta) - a q_2(\eta) \right ) + O \left ( \delta\sqrt{\delta} \right ) i \phi_1(Z_n) + O \left ( \delta^2 \right ) \phi_2(Z_n). \end{align}\] Thus \[\begin{align} \sum_{n=0}^{N-1}{e_n} &= \sum_{n=0}^{N-1}{\log{C95n}} - \frac{1}{2} (\sigma^2-a) \left ( \eta^2 + \xi \eta \delta+ \frac{\xi^2}{3} \delta^2 \right ) T - i \sigma \eta \sqrt{\delta} \sum_{n=0}^{N-1}{Z_n} + i \sigma \xi \delta\sum_{n=0}^{N-1}{W_{t_n}} \\ &= i \left ( \sqrt{\delta} \sigma (q_1(\eta) - \eta) \sum_{n=0}^{N-1}{Z_n} \right ) + \frac{1}{2} (a - \sigma^2) \left ( \eta^2 + \xi \eta \delta+ \frac{\xi^2}{3} \delta^2 \right ) T + \frac{1}{2} \left ( \sigma^2 q_1^2(\eta) - a q_2(\eta) \right ) T \\ &+ O \left ( \delta\sqrt{\delta} i \sum_{n=0}^{N-1}{\phi_1(Z_n)} \right ) + O \left ( \delta^2 \sum_{n=0}^{N-1}{\phi_2(Z_n)} \right ), \end{align}\] since \(\sum_{n=0}^{N-1}{W_{t_n}} = \sqrt{\delta} \sum_{n=0}^{N-1}{\sum_{l=0}^n{Z_l}}\).

Now \[\begin{align} &\exp{\left ( \sum_{n=0}^{N-1}{e_n} \right )} = \exp{\left ( \frac{1}{2} (a - \sigma^2) \left ( \eta^2 + \xi \eta \delta+ \frac{\xi^2}{3} \delta^2 \right ) T + \frac{1}{2} \left ( \sigma^2 q_1^2(\eta) - a q_2(\eta) \right ) T \right )} \\ &\cdot \exp \left ( i \left ( \sqrt{\delta} \sigma (q_1(\eta) - \eta) \sum_{n=0}^{N-1}{Z_n} \right ) + O \left ( \delta\sqrt{\delta} i \sum_{n=0}^{N-1}{\phi_1(Z_n)} \right ) + O \left ( \delta^2 \sum_{n=0}^{N-1}{\phi_2(Z_n)} \right ) \right ). \end{align}\] We have \[\begin{align} q_1(\eta) & = \frac{\sin(\eta h_v)}{h_v} = \eta - \frac{1}{6} \eta^3 h_v^2 + O(\eta^5 h_v^4); \\ q_2(\eta) & = 4 \frac{\sin^2 \left ( \frac{\eta h_v}{2} \right)}{h_v^2} = \eta^2 - \frac{1}{12} \eta^4 h_v^2 + O(\eta^6 h_v^4); \\ q_1^2(\eta) & = \frac{\sin^2(\eta h_v)}{h_v^2} = \eta^2 - \frac{1}{3} \eta^4 h_v^2 + O(\eta^6 h_v^4). \end{align}\] Thus \[\begin{align} &\exp{\left ( \frac{1}{2} (a - \sigma^2) \left ( \eta^2 + \xi \eta \delta+ \frac{\xi^2}{3} \delta^2 \right ) T + \frac{1}{2} \left ( \sigma^2 q_1^2(\eta) - a q_2(\eta) \right ) T \right )} \\ &= \exp{\left ( - \frac{1}{6} \left ( \sigma^2 - \frac{a}{4} \right ) \eta^4 h_v^2 T + \frac{1}{2} (a - \sigma^2) \left ( \xi \eta \delta+ \frac{\xi^2}{3} \delta^2 \right ) T + O(\eta^6 h_v^4) \right )} \\ &= 1 - \frac{1}{6} \left ( \sigma^2 - \frac{a}{4} \right ) \eta^4 h_v^2 T + \frac{1}{2} (a - \sigma^2) \left ( \xi \eta \delta+ \frac{\xi^2}{3} \delta^2 \right ) T + O(\delta h_v^2,\delta^2, h_v^4), \\ \end{align}\] and \[\begin{align} &\exp \left ( i \left ( \sqrt{\delta} \sigma (q_1(\eta) - \eta) \sum_{n=0}^{N-1}{Z_n} \right ) + O \left ( \delta\sqrt{\delta} i \sum_{n=0}^{N-1}{\phi_1(Z_n)} \right ) + O \left ( \delta^2 \sum_{n=0}^{N-1}{\phi_2(Z_n)} \right ) \right ) \\ &= 1 - i \frac{1}{6} \sqrt{\delta} \sigma \eta^3 h_v^2 \sum_{n=0}^{N-1}{Z_n} + O \left ( \delta\sqrt{\delta} i \sum_{n=0}^{N-1}{\phi_1(Z_n)} \right ) + O \left ( \delta^2 \sum_{n=0}^{N-1}{\phi_2(Z_n)} \right ). \end{align}\] Hence \[\begin{align} &X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) = X(T,\xi,\eta+T\xi) \left ( \exp{\left ( \sum_{n=0}^{N-1} e_n \right )} - 1 \right ) \\ &= X(T,\xi,\eta+T\xi) \Bigg ( - i \frac{1}{6} \sigma \eta^3 h_v^2 W_T - \frac{1}{6} \left ( \sigma^2 - \frac{a}{4} \right ) \eta^4 h_v^2 T + \frac{1}{2} (a - \sigma^2) \xi \eta \delta T\\ &+ O \left ( \delta\sqrt{\delta} i \sum_{n=0}^{N-1}{\phi_1(Z_n)} \right ) + O \left ( \delta^2 \sum_{n=0}^{N-1}{\phi_2(Z_n)} \right ) + o(\delta, h_v^2) \varphi(T, h_x \xi, h_v \eta) \Bigg ). \end{align}\]

In the case \((\xi, \eta) \in \Omega_{high}\), we have the following result.

We split the proof of the previous Proposition into steps.

First of all, let us calculate a useful bound for \(\mathbb{E}\left [ |X_N(\xi,\eta)|^2 \right ]\).

By 11 we have that \(X_N = X_0 \prod_{n=0}^{N-1} C_n\), where \[C_n = \frac{1 + i \sqrt{\delta} \sigma Z_n q_1(\eta) - \delta\frac{\sigma^2}{2} \left ( Z_n^2 - 1 \right ) q_1^2(\eta)}{1 + \frac{\delta}{2} a q_2(\eta)}.\] Now, we have \[\begin{align} \mathbb{E}\left [ |C_n|^2 \right ] &= 1 - \frac{\delta\left ( a q_2(\eta) - \sigma^2 q_1^2(\eta) \right ) + \delta^2 \left ( \frac{a^2}{4} q_2^2(\eta) - \frac{\sigma^4}{2} q_1^4(\eta) \right )}{1 + \delta a q_2(\eta) + \frac{\delta^2}{4} a^2 q_2^2(\eta)} \\ &= 1 - \frac{\lambda\left ( 4 a \sin^2 \left ( \frac{\eta h_v}{2} \right ) - \sigma^2 \sin^2 (\eta h_v) \right ) + \lambda^2 \left ( 4 a^2 \sin^4 \left ( \frac{\eta h_v}{2} \right ) - \frac{\sigma^4}{2} \sin^4 (\eta h_v) \right )}{1 + \lambda 4 a \sin^2 \left ( \frac{\eta h_v}{2} \right ) + \lambda^2 4 a^2 \sin^4 \left ( \frac{\eta h_v}{2} \right )} \\ &= 1 - \frac{4 \lambda\sin^2 \left ( \frac{\eta h_v}{2} \right ) \left ( a - \sigma^2 \cos^2 \left ( \frac{\eta h_v}{2} \right ) \right ) + 4 \lambda^2 \sin^4 \left ( \frac{\eta h_v}{2} \right ) \left ( a^2 - 2 \sigma^4 \cos^4 \left ( \frac{\eta h_v}{2} \right ) \right )}{1 + \lambda 4 a \sin^2 \left ( \frac{\eta h_v}{2} \right ) + \lambda^2 4 a^2 \sin^4 \left ( \frac{\eta h_v}{2} \right )} \\ &\leqslant 1 - \frac{4 \lambda a \sin^2 \left ( \frac{\eta h_v}{2} \right ) \left ( 1 - \frac{\sigma^2}{a} \right ) + 4 \lambda^2 a^2 \sin^4 \left ( \frac{\eta h_v}{2} \right ) \left ( 1 - \frac{2 \sigma^4}{a^2} \right )}{1 + \lambda 4 a \sin^2 \left ( \frac{\eta h_v}{2} \right ) + \lambda^2 4 a^2 \sin^4 \left ( \frac{\eta h_v}{2} \right )}. \end{align}\] Let us define \(d := a \sin^2 \left ( \frac{\eta h_v}{2} \right )\) and \(\beta:= {\mathord{{\rm min}}}\left \{ 1 - \frac{\sigma^2}{a}, 1 - \frac{2 \sigma^4}{a^2} \right \} \in (0,1)\), since by hypothesis \(\sigma^2 < a\) and \(2 \sigma^4 < a^2\). We get \[\begin{align} \mathbb{E}\left [ |C_n|^2 \right ] \leqslant 1 - \frac{4 \beta(\lambda d + \lambda^2 d^2)}{(1 + 2 \lambda d)^2}, \end{align}\] and thus ?? holds.

Now, we consider the case \(0 < \lambda< 1\), and thus \[\Omega_{high} = \{ (\xi,\eta) \in \mathbb{R}^2 : |\eta| > h_v^{-2p} \} \cap [-\pi h_x^{-1},\pi h_x^{-1}] \times [-\pi h_v^{-1},\pi h_v^{-1}],\] in particular \(|\Omega_{high}| = 4 \pi h_x^{-1}(\pi h_v^{-1} - h_v^{-2p})\).

We have \[X(T,\xi,\eta) = \exp \left ( \frac{1}{2} (\sigma^2-a) \left ( \eta^2 T - \xi \eta T^2 + \frac{\xi^2}{3} T^3 \right ) + i \sigma (\eta - T \xi) W_T \right ).\] Then \[\begin{align} &\mathbb{E}\left [ \left | \int\!\!\!\int_{\Omega_{high}} X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) d\xi d\eta \right |^2 \right ] \\ &\leqslant 4 \pi^2 h_x^{-1} h_v^{-1} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_{high}} \left | X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) \right |^2 d\xi d\eta \right ] \\ &\leqslant 8 \pi^2 h_x^{-1} h_v^{-1} \int\!\!\!\int_{\Omega_{high}} \mathbb{E}\left [ \left | X^N(\xi,\eta + N \delta\xi) \right |^2 + \left | X(T,\xi,\eta+T\xi) \right |^2 \right ] d\xi d\eta \\ &= 8 \pi^2 h_x^{-1} h_v^{-1} \int\!\!\!\int_{\Omega_{high}}{ \mathbb{E}\left [ \left | X^N(\xi,\eta + N \delta\xi) \right |^2 \right ] d\xi d\eta} + f_0(\delta), \end{align}\] where \[\begin{align} f_0(\delta) &= 8 \pi^2 h_x^{-1} h_v^{-1} \int\!\!\!\int_{\Omega_{high}}{ \exp \left ( (\sigma^2-a) \left ( (\eta+T\xi)^2 T - \xi (\eta+T\xi) T^2 + \frac{\xi^2}{3} T^3 \right ) \right ) d\xi d\eta } \\ &= 8 \pi^2 h_x^{-1} h_v^{-1} \int\!\!\!\int_{\Omega_{high}} {\exp \left ( (\sigma^2-a) \left ( \eta^2 T + \xi \eta T^2 + \frac{\xi^2}{3} T^3 \right ) \right ) d\eta d\xi} \\ &\leqslant C h_x^{-2} h_v^{-2} \exp \left ( (\sigma^2-a) C h_v^{-4p} \right ) d\eta d\xi = o(h_v^r), \quad \forall r > 0, \end{align}\] since \({\mathord{{\rm min}}}_{\Omega_{high}}{\left ( \eta^2 T + \xi \eta T^2 + \frac{\xi^2}{3} T^3 \right )} \geqslant C h_v^{-4p}\), with \(C > 0\).

From ?? , we have \[\mathbb{E}\left [ |X^N(\xi,\eta + N \delta\xi)|^2 \right ] \leqslant\left ( 1 - \frac{4 \beta(d + \lambda d^2)}{(1 + 2 \lambda d)^2} \frac{T}{N h_v^2} \right )^N < \exp{\left ( -\frac{4 \beta(d + \lambda d^2) T}{(1 + 2 \lambda d)^2} h_v^{-2} \right )}\] where \[d := a \sin^2 \left ( \frac{\eta h_v}{2} \right ) \geqslant a \sin^2 \left ( \frac{h_v^{1-2p}}{2} \right ) = a \left ( \frac{h_v^{2-4p}}{4} - \frac{h_v^{4-8p}}{48} + O (h_v^{1-2p}) \right ),\] since \(|\eta| \in [h_v^{-2p},\pi h_v^{-1}]\).

Therefore, \[\begin{align} \mathbb{E}\left [ |X^N(\xi,\eta + N \delta\xi)|^2 \right ] &< \exp{\left ( -\frac{4 \beta(d + \lambda d^2) T}{(1 + 2 \lambda d)^2} h_v^{-2} \right )} \\ &< \exp{\left ( -4\beta d T h_v^{-2} \right )} < \exp{\left ( -a\beta T h_v^{-4p} \right )} = o(h_v^r), \quad \forall r > 0. \end{align}\] This concludes the proof.

For \(\lambda\geqslant 1\), we split \(\Omega_{high}\) into two region \[\Omega_1 := \{ ( |\xi|, |\eta|) \in [0, \pi h_x^{-1}] \times [\delta^{-p}, \delta^{-\frac{1}{2}}] \},\] and \[\Omega_2 := \{ ( |\xi|, |\eta|) \in [0, \pi h_x^{-1}] \times [\delta^{-\frac{1}{2}}, \pi h_v^{-1}] \}.\]

\[\begin{align} \mathbb{E}&\left [ \left | \int\!\!\!\int_{\Omega_{high}} X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) d\xi d\eta \right |^2 \right ] \\ &\leqslant 2 \mathbb{E}\left [ \left | \int\!\!\!\int_{\Omega_1} X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) d\xi d\eta \right |^2 \right ] \\ &+ 2 \mathbb{E}\left [ \left | \int\!\!\!\int_{\Omega_2} X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) d\xi d\eta \right |^2 \right ] \\ &\leqslant 8 \pi h_x^{-1} \delta^{-\frac{1}{2}} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_1} \left |X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) \right |^2 d\xi d\eta \right ] \\ &+ 8 \pi^2 h_x^{-1} h_v^{-1} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_2} \left | X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) \right |^2 d\xi d\eta \right ] \\ &\leqslant 16 \pi h_x^{-1} \delta^{-\frac{1}{2}} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_1} \left |X^N(\xi,\eta + N \delta\xi) \right |^2 + \left | X(T,\xi,\eta+T\xi) \right |^2 d\xi d\eta \right ] \\ &+ 16 \pi^2 h_x^{-1} h_v^{-1} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_2} \left | X^N(\xi,\eta + N \delta\xi) \right |^2 + \left | X(T,\xi,\eta+T\xi) \right |^2 d\xi d\eta \right ] \\ &= 16 \pi h_x^{-1} \delta^{-\frac{1}{2}} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_1} \left |X^N(\xi,\eta + N \delta\xi) \right |^2 d\xi d\eta \right ] \\ &+ 16 \pi^2 h_x^{-1} h_v^{-1} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_2} \left | X^N(\xi,\eta + N \delta\xi) \right |^2 d\xi d\eta \right ] + f_1(\delta), \end{align}\] where \[\begin{align} f_1(\delta) &:= 16 \pi h_x^{-1} \delta^{-\frac{1}{2}} \int\!\!\!\int_{\Omega_1} \exp \left ( (\sigma^2-a) \left ( \eta^2 T + \xi \eta T^2 + \frac{\xi^2}{3} T^3 \right ) \right ) d\xi d\eta \\ &+ 16 \pi^2 h_x^{-1} h_v^{-1} \int\!\!\!\int_{\Omega_2} \exp \left ( (\sigma^2-a) \left ( \eta^2 T + \xi \eta T^2 + \frac{\xi^2}{3} T^3 \right ) \right ) d\xi d\eta \\ &\leqslant C \pi h_x^{-2} \delta^{-1} \exp \left ( -a \beta C \delta^{-2p} \right ) + C h_x^{-2} h_v^{-2} \exp \left ( -a \beta C \delta^{-1} \right ), \end{align}\] since \({\mathord{{\rm min}}}_{\Omega_1}{\left ( \eta^2 T + \xi \eta T^2 + \frac{\xi^2}{3} T^3 \right )} \geqslant C \delta^{-1}\), and \({\mathord{{\rm min}}}_{\Omega_2}{\left ( \eta^2 T + \xi \eta T^2 + \frac{\xi^2}{3} T^3 \right )} \geqslant C \delta^{-\frac{1}{2}}\) , \(C > 0\).

From ?? , we have \[\mathbb{E}\left [ |X^N(\xi,\eta + N \delta\xi)|^2 \right ] \leqslant\left ( 1 - \frac{4 \beta(d + \lambda d^2)}{(1 + 2 \lambda d)^2} \frac{T}{N h_v^2} \right )^N < \exp{\left ( -\frac{4 \beta(d + \lambda d^2) T}{(1 + 2 \lambda d)^2} h_v^{-2} \right )}\] where \[d := a \sin^2 \left ( \frac{\eta h_v}{2} \right ) \geqslant a \sin^2 \left ( \frac{\delta^{-p} h_v}{2} \right ) = a \left ( \frac{\delta^{-2p} h_v^2}{4} - \frac{\delta^{-4p} h_v^4}{48} + O (\delta^{-5p} h_v^5) \right ),\] since \(|\eta| \in [\delta^{-p},\delta^{-\frac{1}{2}}]\).

Therefore, \[\begin{align} \mathbb{E}\left [ |X^N(\xi,\eta + N \delta\xi)|^2 \right ] &< \exp{\left ( -\frac{4 \beta(d + \lambda d^2) T}{(1 + 2 \lambda d)^2} h_v^{-2} \right )} \\ &< \exp{\left ( -4\beta d T h_v^{-2} \right )} < \exp{\left ( -a\beta T \delta^{-2p} \right )}, \end{align}\] and \[\begin{align} 16 \pi h_x^{-1} \delta^{-\frac{1}{2}} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_1} \left |X^N(\xi,\eta + N \delta\xi) \right |^2 d\xi d\eta \right ] \leqslant 64 \pi^2 h_x^{-2} \delta^{-1} \exp{\left ( -a\beta T \delta^{-2p} \right )}. \end{align}\] For \((\xi, \eta) \in \Omega_2\), \(d \in [a \sin{\left ( \frac{1}{2 \sqrt{\lambda}} \right )}, a]\), \[\begin{align} \max_{d}{\left ( 1 - \frac{4 \beta(\lambda d + \lambda^2 d^2)}{(1 + 2 \lambda d)^2} \right )} &= 1 - \beta{\mathord{{\rm min}}}_{d}{\left ( 1 - \frac{1}{(1 + 2 \lambda d)^2} \right )} \\ &= 1 - \beta\left ( 1 - \max_{d}{\frac{1}{(1 + 2 \lambda d)^2}} \right ). \end{align}\] Now, we have that \[\max_{d}{\frac{1}{(1 + 2 \lambda d)^2}} = \frac{1}{ \left (1 + 2 \lambda a \sin{\left ( \frac{1}{2 \sqrt{\lambda}} \right )} \right )^2 } =: \alpha\in (0,1),\] and thus \[1 - \frac{4 \beta(\lambda d + \lambda^2 d^2)}{(1 + 2 \lambda d)^2} \leqslant 1 - \alpha\beta=: \theta_0 \in (0,1).\] Then \[\mathbb{E}\left [ |X^N(\xi,\eta + N \delta\xi)|^2 \right ] < \left ( 1 - \frac{4 \beta(\lambda d + \lambda^2 d^2)}{(1 + 2 \lambda d)^2} \right )^N \leqslant\theta_0^N.\] So \[16 \pi^2 h_x^{-1} h_v^{-1} \mathbb{E}\left [ \int\!\!\!\int_{\Omega_2} \left | X^N(\xi,\eta + N \delta\xi) \right |^2 d\xi d\eta \right ] \leqslant 64 \pi^4 h_x^{-2} h_v^{-2} \theta_0^N.\] Hence \[\begin{align} \mathbb{E}&\left [ \left | \int\!\!\!\int_{\Omega_{high}} X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) d\xi d\eta \right |^2 \right ] \\ &\leqslant C \pi h_x^{-2} \delta^{-1} \exp \left ( -a \beta C \delta^{-2p} \right ) + C h_x^{-2} h_v^{-2} \exp \left ( -a \beta C \delta^{-1} \right ) \\ &+ 64 \pi^2 h_x^{-2} \delta^{-1} \exp{\left ( -a\beta T \delta^{-2p} \right )} + 64 \pi^4 h_x^{-2} h_v^{-2} \theta_0^N \\ &= C \pi h_v^{-2} \delta^{-3} \exp \left ( -a \beta C \delta^{-2p} \right ) + C h_v^{-4} \delta^{-2} \exp \left ( -a \beta C \delta^{-1} \right ) \\ &+ 64 \pi^2 h_v^{-2} \delta^{-3} \exp{\left ( -a\beta T \delta^{-2p} \right )} + 64 \pi^4 h_v^{-4} \delta^{-2} \theta_0^N. \end{align}\] Since the first three terms in the previous inequality are terms of order higher than \(h_v^{-4} \delta^{-2} \theta_0^N\), we have, for \(\lambda\geqslant 1\), \[\mathbb{E}\left [ \left | \int\!\!\!\int_{\Omega_{high}} X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) d\xi d\eta \right |^2 \right ] \leqslant C h_v^{-4} \delta^{-2} \theta_0^N.\] Letting \(\theta := \sqrt{\theta_0}\), the thesis follows.

The result is a consequence of Lemma 2 and Lemma 3.

Now we are ready to prove Theorem 2.

The inverse Fourier transform gives \[\begin{align} U_{j,k}^N &= \frac{1}{4\pi^2} \int_{-\frac{\pi}{h_v}}^{\frac{\pi}{h_v}}{\int_{-\frac{\pi}{h_x}}^{\frac{\pi}{h_x}}{\hat{U}^N(\xi,\eta) e^{i \left ( (j-j_0) h_x \xi + (k-k_0) h_v \eta \right )} d\xi d\eta}} \\ &= \frac{1}{4\pi^2} \int_{-\frac{\pi}{h_v}}^{\frac{\pi}{h_v}}{\int_{-\frac{\pi}{h_x}}^{\frac{\pi}{h_x}}{X^N(\xi,\eta + N \delta\xi) e^{i \left ( (j-j_0) h_x \xi + (k-k_0) h_v \eta \right )} d\xi d\eta}} \end{align}\] and \[\begin{align} u_T(x_j,v_k) &= \frac{1}{4 \pi^2}\int_{\mathbb{R}^2}{\hat{u}_T(\xi,\eta) e^{i \langle (x_j,v_k),(\xi,\eta) \rangle} d\xi d\eta} \\ &= \frac{1}{4 \pi^2}\int_{\mathbb{R}^2}{X(T,\xi,\eta+T\xi) e^{i \left ( (j-j_0) h_x \xi + (k-k_0) h_v \eta \right )} d\xi d\eta}. \end{align}\] Hence, we have \[U_{j,k}^N - u_T(x_j,v_k) = \frac{1}{4 \pi^2}\int\!\!\!\int_{\Omega_{low} \cup \Omega_{high}}{\left ( X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi) \right ) e^{i \langle (x_j-x_0,v_k-v_0),(\xi,\eta) \rangle} d\xi d\eta},\] and the thesis follows by Proposition 3 and Proposition 4.

By Parseval’s Theorem, we have \[\sum_{i,j} h_x h_v \left | U^N_{i,j} - u(T,x_i,y_j) \right |^2 = \int\!\!\!\int_{\Omega_{low} \cup \Omega_{high}}{|\hat{u}(0,\xi,\eta + T \xi)|^2 |X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)|^2 d\xi d\eta}.\] By Proposition 3, we have \(|X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)| = |X(T,\xi,\eta+T\xi)| \left ( O(h_v^2) + O(\delta) \right )\), for \((\xi,\eta) \in \Omega_{low}\). By Proposition 4, we have \(|X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)|^2 \leqslant C \theta^{2N}\), for \((\xi,\eta) \in \Omega_{high}\) with \(C>0\).

Now, if the initial datum is a Dirac delta function, then \(|\hat{u}(0,\xi,\eta + T \xi)| = 1\), and we have \[\begin{align} &\int\!\!\!\int_{\Omega_{low} \cup \Omega_{high}}{|\hat{u}(0,\xi,\eta + T \xi)|^2 |X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)|^2 d\xi d\eta} \\ &= \int\!\!\!\int_{\Omega_{low}}{|X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)|^2 d\xi d\eta} + \int\!\!\!\int_{\Omega_{high}}{|X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)|^2 d\xi d\eta} \\ &= O(h_v^4) + O(\delta^2) + O(h_v^{-2} \delta^{-1} \theta^{2N}). \end{align}\] If the initial datum is in \(L_2\), then \(\int \int{|\hat{u}(0,\xi,\eta + T \xi)|^2 d\xi d\eta} < + \infty\) and \[\int\!\!\!\int_{\Omega_{high}}{|\hat{u}(0,\xi,\eta + T \xi)|^2 |X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)|^2 d\xi d\eta} = O(\delta^r), \qquad \text{for any } r > 0,\] so \[\int\!\!\!\int_{\Omega_{low} \cup \Omega_{high}}{|\hat{u}(0,\xi,\eta + T \xi)|^2 |X^N(\xi,\eta + N \delta\xi) - X(T,\xi,\eta+T\xi)|^2 d\xi d\eta} = O(h_v^4) + O(\delta^2).\]

4 Multi-level Monte Carlo↩︎

In this section, we consider the problem of estimating the expectation of a Lipschitz function \(P\) of the solution of the SPDE 1 .

Suppose that we want to approximate \(P\) using a standard Monte Carlo (MC) approach (i.e. averaging over \(N\) independent SPDE simulations). To achieve a \(RMSE\) of order \(\varepsilon\) we would need \(N = O(\varepsilon^{-2})\), \(h_v = O(\varepsilon^{\frac{1}{2}})\), \(\delta= O(\varepsilon)\) and consequently \(h_x = O(\varepsilon^{\frac{3}{2}})\), so the computational cost would be \(O(N h_x^{-1} h_v^{-1} \delta^{-1}) = O(\varepsilon^{-5})\).

Since the complexity of a standard MC is quite high, we aim to improve it by using a Multi-level Monte Carlo (MLMC) strategy. Following a standard approach (see for example [4], [12]), we obtain this result.

Let \(N > 0, h_l := \frac{T}{N^l}\), and \(L = \left \lceil \frac{\log \left ( \sqrt{2} c_1 T^\alpha\varepsilon^{-1} \right )}{\alpha\log(N)} \right \rceil\), (where \(\lceil x \rceil\) denotes the integer part of a real number \(x\)). Then, we have \[N^{-\alpha} \frac{1}{\sqrt{2}} \varepsilon< c_1 h_L^\alpha\leqslant\frac{1}{\sqrt{2}} \varepsilon,\] and, by [i] and [ii] \[\label{MLMC95mean95estimate} \left ( \mathbb{E}[\hat{Y}] - \mathbb{E}[P] \right )^2 \leqslant\frac{1}{2} \varepsilon^2.\tag{13}\] Moreover, we have the following useful inequality \[\sum_{l=0}^L{h_l^{-\gamma}} = h_L^{-\gamma} \sum_{l=0}^L{N^{-\gamma l}} < h_L^{-\gamma} \frac{N^\gamma}{N^\gamma-1},\] and \[\label{hL95estimate} h_L^{-\gamma} < N^\gamma\left ( \frac{\varepsilon}{\sqrt{2} c_1} \right )^{-\frac{\gamma}{\alpha}}.\tag{14}\] Hence \[\label{sum95hl95estimate} \sum_{l=0}^L{h_l^{-\gamma}} < \frac{N^{2\gamma}}{N^\gamma-1} \left ( \sqrt{2} c_1 \right )^{\frac{\gamma}{\alpha}} \varepsilon^{-\frac{\gamma}{\alpha}}.\tag{15}\] Now, let us set \(N_l := \left \lceil 2 \varepsilon^{-2} c_2 h_L^{-\frac{\gamma-\beta}{2}} (1-N^{-\frac{\gamma-\beta}{2}})^{-1} h_l^{\frac{\gamma+\beta}{2}} \right \rceil\). Then, by [iii] \[\sum_{l=0}^L{{\mathrm{var}}(\hat{Y}_l)} < \frac{1}{2} \varepsilon^2 h_L^{\frac{\gamma-\beta}{2}} (1-N^{-\frac{\gamma-\beta}{2}}) \sum_{l=0}^L h_l^{-\frac{\gamma-\beta}{2}} < \frac{1}{2} \varepsilon^2,\] since \[\sum_{l=0}^L h_l^{-\frac{\gamma-\beta}{2}} < h_L^{-\frac{\gamma-\beta}{2}} (1-N^{-\frac{\gamma-\beta}{2}})^{-1}.\] Therefore, by the last inequality and the estimate 13 , we have an \(\varepsilon^2\) upper bound for the \(MSE\).

Finally, we have \[N_l < 2 \varepsilon^{-2} c_2 h_L^{-\frac{\gamma-\beta}{2}} (1-N^{-\frac{\gamma-\beta}{2}})^{-1} h_l^{\frac{\gamma+\beta}{2}} + 1,\] so, by [iv] and 15 , the computational complexity is bounded by \[\begin{align} C &\leqslant c_3 \sum_{l=0}^L N_l h_l^{-\gamma} \\ &< c_3 \left ( 2 \varepsilon^{-2} c_2 h_L^{-\frac{\gamma-\beta}{2}} (1-N^{-\frac{\gamma-\beta}{2}})^{-1} \sum_{l=0}^L h_l^{-\frac{\gamma-\beta}{2}} + \sum_{l=0}^L h_l^{-\gamma} \right ) \\ &< c_3 \left ( 2 \varepsilon^{-2} c_2 (1-N^{-\frac{\gamma-\beta}{2}})^{-2} h_L^{-(\gamma-\beta)} + \frac{N^{2\gamma}}{N^\gamma-1} \left ( \sqrt{2} c_1 \right )^{\frac{\gamma}{\alpha}} \varepsilon^{-\frac{\gamma}{\alpha}} \right ). \end{align}\] Now, by 14 we have \(h_L^{-(\gamma-\beta)} < N^{\gamma-\beta} (\sqrt{2} c_1)^{\frac{\gamma-\beta}{\alpha}} \varepsilon^{-\frac{\gamma-\beta}{\alpha}}\). Thus \[\begin{align} C &< c_3 \left ( 2 c_2 (1-N^{-\frac{\gamma-\beta}{2}})^{-2} N^{\gamma-\beta} (\sqrt{2} c_1)^{\frac{\gamma-\beta}{\alpha}} \varepsilon^{-2 - \frac{\gamma-\beta}{\alpha}} + \frac{N^{2\gamma}}{N^\gamma-1} \left ( \sqrt{2} c_1 \right )^{\frac{\gamma}{\alpha}} \varepsilon^{-\frac{\gamma}{\alpha}} \right ) \\ &\leqslant c_4 \left ( \varepsilon^{-2 - \frac{\gamma-\beta}{\alpha}} + \varepsilon^{-\frac{\gamma}{\alpha}} \right ), \end{align}\] where \(c_4 := c_3 \cdot \max{\{ 2 c_2 (1-N^{-\frac{\gamma-\beta}{2}})^{-2} N^{\gamma-\beta} (\sqrt{2} c_1)^{\frac{\gamma-\beta}{\alpha}}, \frac{N^{2\gamma}}{N^\gamma-1} \left ( \sqrt{2} c_1 \right )^{\frac{\gamma}{\alpha}} \}}\).

Now, let us set \(h_v \approx h^l\), \(\delta\approx h^{2l}\) and \(h_x \approx h^{3l}\) for the level \(l\) in the MLMC, for some \(h<1\). Then, the computational cost at level \(l\) is \(C_l \leqslant c_3 N_l h^{-6l}\), i.e. \(\gamma= 6\).

Moreover, if \(P\) is a Lipschitz function of the SPDE solution, then Corollary 1 implies that the weak and strong errors for the level \(l\) satisfy \[|\mathbb{E}[\hat{P}_l - P]| \leqslant c_1 h^{2 l},\] and \[\begin{align} {\mathrm{var}}{(\hat{P}_l - \hat{P}_{l-1})} &\leqslant{\mathrm{var}}{(\hat{P}_l - P)} + {\mathrm{var}}{(\hat{P}_{l-1} - P)} \\ &\leqslant\mathbb{E}[(\hat{P}_l - P)^2] + \mathbb{E}[(\hat{P}_{l-1} - P)^2] \leqslant c_2 h^{4l}, \end{align}\] so \(\alpha= 2\) and \(\beta= 4 < \gamma\). Notice that we are not in the standard case \(\alpha> \frac{1}{2} \gamma\).

Therefore, following a MLMC strategy, Theorem 5 assures that the complexity for achieving a \(RMSE\) of order \(\varepsilon\) is \(C \leqslant c_4 \left ( \varepsilon^{-2 - \frac{\gamma-\beta}{\alpha}} + \varepsilon^{- \frac{\gamma}{\alpha}} \right ) = 2 c_4 \varepsilon^{-3}\).

5 Numerical tests↩︎

We numerically verify the order of convergence by approximating the \(L_2\) error \[\left ( \sum_{i,j} h_x h_v \mathbb{E}\left [ \left | U^N_{i,j} - u(T,x_i,y_j) \right |^2 \right ] \right )^{\frac{1}{2}}\] by \[\label{L95295err95approx} \left ( \sum_{i,j} h_x h_v E_L \left [ \left | U^N_{i,j} - u(T,x_i,y_j) \right |^2 \right ] \right )^{\frac{1}{2}} := \left ( \sum_{i,j} \frac{h_x h_v}{L} \sum_{l=1}^L \left | U^N_{i,j}(W^{(l)}) - u(T,x_i,y_j;W^{(l)}) \right |^2 \right )^{\frac{1}{2}},\tag{16}\] where \(W^{(l)}\) are independent Brownian motions and \(E_L\) is the empirical mean with \(L\) samples.

Figure [fig1] illustrates the empirical strong error between numerical solutions computed on consecutively refined grids. For each parameter (either \(\delta\) or \(h_v\)), we consider a sequence of discretizations and calculate the differences between successive approximations, resulting in a set of points that quantify the convergence behavior of the numerical method as the grid is refined. The results show that the observed convergence rates are consistent with the theoretical analysis, exhibiting first-order behavior with respect to \(\delta\) and second-order behavior with respect to \(h_v\). This provides a clear visual confirmation that the numerical method achieves the expected accuracy as predicted by the theoretical estimates.

Figure 2 displays the numerical approximation of the solution at the final time \(T = 1\) together with the corresponding pointwise error with respect to the exact solution. In the left panel we report the approximate solution obtained with discretization parameters \(h_v = 10^{-1}\) and \(\delta= 10^{-2}\), while the right panel shows the difference between the approximate and exact solution on the same grid. The initial condition is centered at \(x_0 = 0.1\) and \(v_0 = 5\). The comparison highlights both the structure of the numerical approximation and the regions where the error is most pronounced, providing a clear indication of the accuracy of the scheme under the chosen discretization.

5.1 Comparison of MC and MLMC for the numerical solution of the SPDE↩︎

In this experiment, we investigate the efficiency of the Multilevel Monte Carlo (MLMC) method compared to the standard Monte Carlo (MC) approach for approximating the expected value of a functional of the solution of the SPDE. We performed two tests, corresponding to different spatial resolutions: \(n = 4\) and \(n = 8\), with \(\varepsilon= n^{-2}\). Consequently, the number of independent runs for the standard MC approach was set as \(D = \varepsilon^{-2}\), resulting in \(D = 256\) for \(n = 4\) and \(D = 4096\) for \(n = 8\).

For the MLMC simulations, we first estimated the constants \(c_1\) and \(c_2\) from Theorem 5. Based on preliminary tests:

• For \(n = 4\), we obtained \(L = 2\) levels with \(N_l = [52, 2, 1]\) samples per level.

• For \(n = 8\), we obtained \(L = 3\) levels with \(N_l = [1637, 52, 2, 1]\) samples per level.

We want to point out that once the constants \(c_1\) and \(c_2\) have been estimated, there is no need to repeat the parameter estimation procedure when using different accuracy levels, as the values of \(L\) and \(N_l\) can be scaled accordingly.

The discretization parameters for each level \(l = 0,\dots,L\) were set as \[\tau^{(l)} = 2^{-l}, \qquad h_{v}^{(l)} = 2^{- l / 2},\] ensuring a proper scaling of the temporal and spatial resolutions across levels.

The quantity of interest, or payoff function, was chosen as the \(L_2\)-norm of the SPDE solution, \(P = \|U\|\). The reference exact value \(P_\text{exact}\) was computed as the mean over all \(D\) runs, yielding \[P_\text{exact} = 0.875063 \quad \text{for } n=4, \qquad P_\text{exact} = 0.875063 \quad \text{for } n=8.\]

The results, including the mean estimate, absolute and relative errors, and computational time, are summarized in Table 1. For MLMC, the reported time includes the additional time required to estimate the parameters \(L\) and \(N_l\).

6 Conclusions↩︎

In this work we studied the numerical approximation of degenerate Langevin-type SPDEs in two spatial dimensions, a class of equations that arises naturally in mathematical finance and in the filtering of partially observed diffusions. To handle the mixed deterministic-stochastic structure of the operator \(\mathbf{Y} := \partial_t + v\partial_x\) together with the degeneracy of the diffusion in the \(x\) direction, we introduced a semi-implicit Milstein finite difference scheme, treating the drift term implicitly and the diffusion term with a Milstein-type explicit correction.

Through a Fourier analysis of the scheme, we derived explicit conditions on the coefficients \(a\) and \(\sigma\) (namely \(\sigma^2 < a\) and \(2\sigma^4 < a^2\)) under which the scheme is asymptotically mean-square stable, and we established a full convergence analysis, proving that the scheme achieves order one convergence in the time step \(\delta\) and order two convergence in the velocity mesh size \(h_v\), in the discrete \(L^2\)-in-space, \(L^2\)-in-probability norm. These theoretical rates were confirmed numerically, with the observed convergence orders matching the predicted first-order behavior in \(\delta\) and second-order behavior in \(h_v\).

We then embedded the proposed scheme within a Multilevel Monte Carlo (MLMC) framework in order to reduce the computational cost of estimating expectations of functionals of the SPDE solution. Exploiting the weak and strong error estimates obtained from the Fourier analysis, we showed that the MLMC complexity for achieving a root-mean-square error of order \(\varepsilon\) is \(\mathcal{O}(\varepsilon^{-3})\), a substantial improvement over the \(\mathcal{O}(\varepsilon^{-5})\) complexity of a standard Monte Carlo approach. This theoretical gain was corroborated by numerical experiments, which showed that MLMC achieves comparable or better accuracy than standard Monte Carlo at a fraction of the computational time, with the advantage becoming more pronounced as the target accuracy increases.

Overall, our results demonstrate that combining a semi-implicit Milstein discretization with a multilevel Monte Carlo strategy provides an effective and theoretically grounded approach for the simulation of degenerate Langevin-type SPDEs. Several directions remain open for future research, including the extension of the scheme to the case of non-constant (space- and time-dependent) coefficients \(a_t(x,v)\) and \(\sigma_t(x,v)\), a more refined treatment of boundary effects on bounded domains, and the application of the proposed methodology to the pricing of path-dependent financial derivatives such as Asian options, where the degenerate structure of the SPDE plays a central role.

Acknowledgements↩︎

The first author is member of the INdAM Research Group GNCS.
This research has been mainly developed during a research stay of S. Portaro at University of A Coruña, supported by the Marco Polo Scholarship program of the Alma Mater Studiorum - Università di Bologna. Moroever, S. Portaro gratefully acknowledges the financial support from the Italian Ministry of University and Research (MUR) under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 1, Investment 4.1 Extension of the number of research doctorates and innovative doctorates for public administration and cultural heritage, funded by the European Union - NextGenerationEU. Project Title: Tecniche di Data Science in Infrastrutture HPC per il Monitoraggio Ambientale delle Valli del Reno, Lavino e Samoggia - CUP: J33C23002500002 - Alma Mater Studiorum - Università di Bologna.
C. Vázquez acknowledges the funding from the Galician Government through the grant ED431C 2022/047 (both including FEDER financial support). C. Vázquez also acknowledges the support of CITIC, as a center accredited for excellence within the Galician University System and a member of the CIGUS Network, that receives subsidies from the Department of Education, Science, Universities, and Vocational Training of the Xunta de Galicia. Additionally, CITIC is co-financed by the EU through the FEDER Galicia 2021-27 operational program (Ref. ED431G 2023/01).

References↩︎

---

1. Dipartimento di Matematica, (AM)\(^2\), Alma Mater Studiorum - Università di Bologna, 40126 Bologna, Italy, sascha.portaro@unibo.it↩︎

2. Department of Mathematics and CITIC, University of A Coruña, Campus Elviña, 15071-A Coruña, Spain,carlosv@udc.es↩︎