Analysis of sparse grid multilevel estimators for multi-dimensional Zakai equations

TL;DR

This paper evaluates the accuracy and computational efficiency of multilevel estimators for multi-dimensional Zakai equations, demonstrating that MLMC on sparse grids achieves optimal complexity, especially in higher dimensions.

Contribution

The paper introduces a detailed analysis of sparse grid multilevel Monte Carlo methods for Zakai equations, showing their optimal complexity and advantages over traditional methods.

Findings

MLMC on sparse grids has optimal $O( ext{error}^{-2})$ complexity in 2D.

MLMC on regular grids has $O( ext{error}^{-2}( ext{log error})^2)$ complexity.

Numerical tests confirm the theoretical complexity results.

Abstract

In this article, we analyse the accuracy and computational complexity of estimators for expected functionals of the solution to multi-dimensional parabolic stochastic partial differential equations (SPDE) of Zakai-type. Here, we use the Milstein scheme for time integration and an alternating direction implicit (ADI) splitting of the spatial finite difference discretisation, coupled with the sparse grid combination technique and multilevel Monte Carlo sampling (MLMC). In the two-dimensional case, we find by detailed Fourier analysis that for a root-mean-square error (RMSE) $\varepsilon$, MLMC on sparse grids has the optimal complexity $O(\varepsilon^{-2})$, whereas MLMC on regular grids has $O(\varepsilon^{-2}(\log\varepsilon)^2)$, standard MC on sparse grids $O(\varepsilon^{-7/2}(|\log\varepsilon|)^{5/2})$, and MC on regular grids $O(\varepsilon^{-4})$. Numerical tests confirm these findings empirically. We give a discussion of the higher-dimensional setting without detailed proofs, which suggests that MLMC on sparse grids always leads to the optimal complexity, standard MC on sparse grids has a fixed complexity order independent of the dimension (up to a logarithmic term), whereas the cost of MLMC and MC on regular grids increases exponentially with the dimension.

Tables

Table 1. Table 1: log 2 ⁡ | 𝔼 ​ [ Δ ​ P ( l 1 , l 2 ) N ] | subscript 2 𝔼 delimited-[] Δ subscript superscript 𝑃 𝑁 subscript 𝑙 1 subscript 𝑙 2 \log_{2}\big{|}\mathbb{E}[\Delta P^{N}_{(l_{1},\,l_{2})}]\big{|} and log 2 ⁡ Var ​ [ Δ ​ P ( l 1 ​ l 2 ) N ] subscript 2 Var delimited-[] Δ subscript superscript 𝑃 𝑁 subscript 𝑙 1 subscript 𝑙 2 \log_{2}\mathrm{Var}[\Delta P^{N}_{(l_{1}\,l_{2})}] with h x = 2 − l 1 , h y = 2 − l 2 formulae-sequence subscript ℎ 𝑥 superscript 2 subscript 𝑙 1 subscript ℎ 𝑦 superscript 2 subscript 𝑙 2 h_{x}=2^{-l_{1}},\,h_{y}=2^{-l_{2}} , k = 4 − 3 𝑘 superscript 4 3 k=4^{-3} .

	0	1	2	3	4
$0$	-0.0819	-7.28	-9.04	-11.00	-12.98
$1$	-7.30	-13.73	-15.49	-17.44	-19.43
$2$	-9.05	-15.48	-17.24	-19.19	-21.17
$3$	-10.99	-17.43	-19.19	-21.13	-23.12
$4$	-12.98	-19.42	-21.17	-23.13	-25.11

Table 2. Table 2: Log values of mean, variance, cost for δ ​ L l 𝛿 subscript 𝐿 𝑙 \delta L_{l} , h 0 = 1 / 2 , k 0 = 1 / 8 formulae-sequence subscript ℎ 0 1 2 subscript 𝑘 0 1 8 h_{0}=1/2,\ k_{0}=1/8 .

	$l=0$	$l=1$	$l=2$	$l=3$	$l=4$	$l=5$
${\log }_2\left|𝔼\left[\delta P_l\right]\right|$	-0.0647	-7.62	-9.44	-11.42	-13.30	-15.22
${\log }_2\left(\mathrm{Var}\left[\delta P_l\right]\right)$	-9.15	-16.93	-20.71	-24.66	-28.20	-31.83
${\log }_2\left(\mathrm{Cost}\left[\delta P_l\right]\right)$	13.64	18.50	22.02	25.22	28.44	31.65

Equations

$$\mathrm{d}v$$

$$v(0,x,y)=\delta (x-x_0)\otimes \delta (y-y_0)$$

$$Lv$$

$$v(T,x,y)=\frac{\exp\Big{(}-\frac{\big{(}x-x_{0}-\mu_{x}T-\sqrt{\rho_{x}}W_{T}^{x}\big{)}^{2}}{2(1-\rho_{x})T}-\frac{\big{(}y-y_{0}-\mu_{y}T-\sqrt{\rho_{y}}W_{T}^{y}\big{)}^{2}}{2(1-\rho_{y})T}\Big{)}}{2\pi\sqrt{(1-\rho_{x})(1-\rho_{y})}\,T}\,.$$

$$\mathrm{d}v=-\mu \frac{\partial v}{\partial x}\mathrm{d}t+\frac{1}{2}\frac{{\partial }^{2}v}{\partial x^2}\mathrm{d}t-\sqrt{\rho }\frac{\partial v}{\partial x}\mathrm{d}{W}_t,(t,x)\in (0,T)\times \mathbb{R},$$

$$P_t={\int }_0^{\infty }{\int }_0^{\infty }v(t,x,y)\mathrm{d}x\mathrm{d}y,$$

$$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)v(t,x,y)\,\mathrm{d}x\,\mathrm{d}y,\qquad g\Big{(}\int_{0}^{\infty}\int_{0}^{\infty}v(t,x,y)\,\mathrm{d}x\,\mathrm{d}y\Big{)},$$

$$V_{i,j}^0=h_x^{-1}{h}_y^{-1}{\delta }_{(i_0,j_0)}=\{\begin{array}{ll}{h}_x^{-1}{h}_y^{-1},& i=i_0,j=j_0,\\ 0,& \text{otherwise}.\end{array}$$

$$\displaystyle\!\!\!\!\bigg{(}I+\frac{\mu_{x}k}{2h_{x}}D_{x}+\frac{\mu_{y}k}{2h_{y}}D_{y}-\frac{k}{2h_{x}^{2}}D_{xx}-\frac{k}{2h_{y}^{2}}D_{yy}\!\bigg{)}V^{n+1}\!=\bigg{(}\!I\!-\!\frac{\sqrt{\rho_{x}k}Z_{n,x}}{2h_{x}}D_{x}\!-\!\frac{\sqrt{\rho_{y}k}\widetilde{Z}_{n,y}}{2h_{y}}D_{y}$$

$$\displaystyle+\frac{\rho_{x}k(Z_{n,x}^{2}\!\!-\!1)}{8h_{x}^{2}}D_{x}^{2}\!+\!\frac{\rho_{y}k(\widetilde{Z}_{n,y}^{2}\!\!-\!1)}{8h_{y}^{2}}D_{y}^{2}\!+\!\frac{\sqrt{\rho_{x}\rho_{y}}kZ_{n,x}\widetilde{Z}_{n,y}}{4h_{x}h_{y}}D_{xy}\!\bigg{)}\!V^{n}$$

$$=:M_n{V}^n,$$

$$(D_{x}V{)}_{i,j}=V_{i+1,j}-V_{i-1,j},(D_{y}V{)}_{i,j}=V_{i,j+1}-V_{i,j-1},$$

$$(D_{xx}V{)}_{i,j}=V_{i+1,j}-2V_{i,j}+V_{i-1,j},(D_{yy}V{)}_{i,j}=V_{i,j+1}-2V_{i,j}+V_{i,j-1},$$

$$(D_{xy}V{)}_{i,j}=V_{i+1,j+1}-V_{i-1,j+1}-V_{i+1,j-1}+V_{i-1,j-1},$$

$$\displaystyle\bigg{(}I+\frac{\mu_{x}k}{2h_{x}}D_{x}-\frac{k}{2h_{x}^{2}}D_{xx}\bigg{)}\bigg{(}I+\frac{\mu_{y}k}{2h_{y}}D_{y}-\frac{k}{2h_{y}^{2}}D_{yy}\bigg{)}V^{n+1}=M_{n}V^{n}\!.$$

$$2{\rho }_x^2(1+2|{\rho }_{xy}|)$$

$$k/\min \{h_x^2,h_y^2\}\le \lambda ,$$

$$\sqrt{\mathbb{E}\big{[}|V_{i,j}^{N}-v(T,ih_{x},jh_{y})|^{2}\big{]}}\leq C(h_{x}^{2}+h_{y}^{2}).$$

$$P_T={\int }_0^{\infty }{\int }_0^{\infty }v(T,x,y)\mathrm{d}x\mathrm{d}y,$$

$$P_{(l_1,l_2)}^N$$

$$\sqrt{\mathbb{E}\left[|P_{(l_1,l_2)}^N-P{|}^2\right]}\le C(h_x^2+h_y^2),$$

$${\Delta }_i{P}_{\mathbf{l}}^N=\{\begin{array}{ll}{P}_{\mathbf{l}}^N-P_{\mathbf{l}-{\mathbf{e}}_i}^N,& \text{if}{l}_i>0,\\ P_{\mathbf{l}}^N& \text{if}{l}_i=0,\end{array}$$

$$\Delta P_{(l_1,l_2)}^N=P_{(l_1,l_2)}^N-P_{(l_1,l_2-1)}^N-P_{(l_1-1,l_2)}^N+P_{(l_1-1,l_2-1)}^N.$$

$$\Big{|}\mathbb{E}\left[\Delta P^{N}_{(l_{1},\,l_{2})}\right]\Big{|}=O(h_{x}^{2}h_{y}^{2}),\qquad\mathbb{E}\Big{[}\left|\Delta P^{N}_{(l_{1},\,l_{2})}\right|^{2}\Big{]}=O(h_{x}^{4}h_{y}^{4}).$$

$$\displaystyle\Big{|}\mathbb{E}\left[\Delta_{1}P^{N}_{(l_{1},\,l_{2})}\right]\Big{|}=O(h_{x}^{2}),\qquad\mathbb{E}\Big{[}\left|\Delta_{1}P^{N}_{(l_{1},\,l_{2})}\right|^{2}\Big{]}=O(h_{x}^{4}),$$

$$\displaystyle\Big{|}\mathbb{E}\left[\Delta_{2}P^{N}_{(l_{1},\,l_{2})}\right]\Big{|}=O(h_{y}^{2}),\qquad\mathbb{E}\Big{[}\left|\Delta_{2}P^{N}_{(l_{1},\,l_{2})}\right|^{2}\Big{]}=O(h_{y}^{4}).$$

$$P_l^N=\sum _{(l_1,l_2)\in {\mathcal{I}}_l}\Delta P_{(l_1,l_2)}^N.$$

$$\sqrt{\mathbb{E}\left[|P_l^N-P{|}^2\right]}=O(l{2}^{-2l}).$$

$$\widehat{Y}=\frac{1}{M}\sum_{m=1}^{M}\sum_{(l_{1},\,l_{2})\in\mathcal{I}_{l}}\!\!\!\Big{(}\Delta\widehat{P}_{(l_{1},\,l_{2})}^{N}\Big{)}^{(m)},$$

$$W=\sum_{l_{1}+l_{2}\leq l+1}2^{l_{1}+l_{2}}\cdot O(k^{-1})=O(l\,2^{3l})=O\big{(}\varepsilon^{-\frac{3}{2}}|\log\varepsilon|^{\frac{5}{2}}\big{)}.$$

$$W=O\big{(}\varepsilon^{-\frac{7}{2}}|\log\varepsilon|^{\frac{5}{2}}\big{)}.$$

$${\mathcal{I}}_l=\{(l_1,l_2)\in {\mathbb{N}}_0^2:l_1+l_2\le l+1\},k_l=k_0\cdot 2^{-2l},N_l=\frac{T}{k_l}.$$

$$\mathbb{E}[P]=\mathbb{E}[P_0]+\sum _{l=1}^{\infty }\mathbb{E}[P_l-P_{l-1}].$$

Full text

1 Introduction

The focus of this paper is the efficient simulation of the two-dimensional SPDE

[TABLE]

for $x,y\in\mathbb{R},\ 0<t\leq T$, subject to the Dirac initial datum

[TABLE]

for given $x_{0}$ and $y_{0}$, where $W=(W^{x},\,W^{y})$ is a two-dimensional standard Brownian motion with correlation $\rho_{xy}$ on a probability space $(\Omega,\mathcal{F},\mathbb{P})$,

[TABLE]

$\mu_{x},\mu_{y}$ and $0\leq\rho_{x},\rho_{y}<1$, $-1\leq\rho_{xy}\leq 1$ are real-valued parameters.

This is a special case of the Zakai equation from stochastic filtering where $v$ describes the distribution of the filter given a signal process $W$ (see [1, 9]).

A classical result states that, for a class of SPDEs including (1.1), with initial condition in $L_{2}$, there exists a unique solution $v\in L_{2}(\Omega\times(0,T),\mathcal{F},L_{2}(\mathbb{R}^{2}))$ [13]. This does not include Dirac initial data (1.2), but in fact, the solution to (1.1) and (1.2) can be analytically derived as the smooth (in $x$ and $y$) function

[TABLE]

More commonly, however, such a closed-form solution is not available, for instance in the case of variable coefficients. We will focus on (1.1) for the analysis, but the numerical methods we investigate apply similarly to such a wider class.

There is a large body of recent literature on the numerical solution of SPDEs. Most closely related to the present work, Giles and Reisinger in [8] used an explicit Milstein finite difference approximation to the solution of the one-dimensional SPDE

[TABLE]

where $T>0$, $W$ is a standard Brownian motion, and $\mu$ and $0\leq\rho<1$ are real-valued parameters; [17] extended the discretisation to an implicit method on the basis of the $\sigma$-$\theta$ time-stepping scheme. This 1-d SPDE (1.4) has also been used to model default risk in large credit portfolios (see [4]).

For the 2-d SPDE (1.1), we will use an implicit method such that under some constraints on $\rho_{x},\rho_{y},\rho_{xy}$, the scheme is unconditionally mean-square stable. Furthermore, we use an Alternating Direction Implicit (ADI) factorisation which is more convenient computationally than a purely implicit scheme, and is also unconditionally mean-square stable (see [19]).

We consider the following functional of the solution,

[TABLE]

which is a two-dimensional version of the loss in [4, 8] which represents the proportion of the defaulted firms there. In the context of filtering, $P_{t}$ is related to the cumulative distribution function of the filter given the signal.

The functional in (1.5) is a special case of more general linear and nonlinear functionals of the form

[TABLE]

with $f$ being the Heaviside function and $g$ the identity in the case of $P_{t}$. Preliminary derivations indicate that our analysis may be extended from $P_{t}$ to these cases for sufficiently smooth $f$ and $g$, by judicious multivariate Taylor expansion, but this involves exceedingly lengthy calculus and is beyond the scope of this work.

A classical approach to approximating $\mathbb{E}[P_{t}]$ is the standard Monte Carlo method, using a suitable approximation scheme for (1.1) and sampling of $W$ on a discrete time mesh. For the SPDE (1.1) and standard schemes, to achieve a root mean square error (RMSE) $\varepsilon$, this requires an overall computational cost $O(\varepsilon^{-4})$, as we require $O(\varepsilon^{-2})$ samples, $O(\varepsilon^{-1})$ time steps (e.g., for the Euler-Maruyama scheme with weak order 1), and $O(\varepsilon^{-1/2})$ mesh points in each direction (e.g., for central difference schemes with order 2). One way to reduce the cost is the MLMC method (see [6]) by using the SPDE solution on paths with a coarse timestep and spatial mesh as a control variate of solutions on paths with a fine timestep and mesh. As a result, for a fixed accuracy $\varepsilon$, the cost can be reduced significantly to $O(\varepsilon^{-2}(\log\varepsilon)^{2})$ by standard MLMC methodology as in [8]. However, this complexity of the MLMC method is not optimal as in the one-dimensional case in [8]. The reason is that the cost of each sample on higher levels increases with the same order as the variance decays. Moreover, the total cost of MLMC increases exponentially in the dimension.

The approach taken here is to approximate the SPDE (1.1) by the sparse grid combination technique. Sparse grids were first introduced to solve high-dimensional PDEs on a tensor product space in [20]. The error bounds in [2] show that sparse finite element approximations can alleviate the curse of dimensionality in the numerical implementation of certain elliptic PDEs with sufficiently smooth solutions. In contrast to the finite element method, the combination technique, first introduced in [10], decomposes the solution into contributions from simple tensor product grids with different resolutions in each dimension. For a survey of methods and early results see [3]. The analysis in [15] and [16] shows that the computational cost for the combination technique applied to the Poisson problem with sufficient regularity is independent of the dimension, up to a logarithmic term, for finite elements and finite differences, respectively. Hendricks, Ehrhardt and Günther combine the sparse grid combination technique with the ADI scheme for diffusion equations in [12].

This sparse grid method has been extended to multi-index Monte Carlo (MIMC) in the context of SPDEs in [11]. MIMC can be viewed as the sparse grid combination technique applied to equations with stochasticity, with optimised number of samples for the individual terms in the combination formula (the “hierarchical surpluses”), akin MLMC. Giles, Kuo and Sloan summarised these ideas applied to elliptic PDEs with finite-dimensional uncertainty in the coefficients in [7]. The MIMC method was applied in [18] to a 1-d SPDE (1.4), where the timestep and space mesh are coupled for stability. However, optimal complexity is not achieved in this space-time method.

In this paper, given a fixed timestep and Brownian path, we solve the SPDE using the sparse grid combination technique in space. Then, to evaluate $\mathbb{E}[P_{t}]$, we use $M$ independent samples of the hierarchical surpluses and calculate the average. The benefit here is that with a RMSE $\varepsilon$, the total cost is fixed with $O(\varepsilon^{-7/2})$ up to a logarithmic term, as the cost for one sample is $O(\varepsilon^{-3/2}|\log\varepsilon|^{5(d-1)/2})$, and the number of samples needed is $M=O(\varepsilon^{-2})$. Hence, this will improve on the complexity of the MLMC method, whose total cost is $O(\varepsilon^{-1-d/2})$, when the dimension $d>5$.

To recover the optimal complexity $O(\varepsilon^{-2})$, we further combine the sparse combination technique and MLMC, in a different way from standard MIMC. In this way, the total cost is $O(\varepsilon^{-2})$ independent of dimension.

The rest of this article is structured as follows. We define the approximation schemes in Section 2. Section 3 gives a Fourier analysis of the sparse combination estimators. Section 4 shows numerical experiments confirming the above findings, and Section 5 generalises the problem to higher dimensions. Section 6 offers conclusions and directions for further research.

2 Approximation and main results

For simplicity of presentation, we initially restrict ourselves with the description of the schemes and their analysis to the two-dimensional case. The extension to higher dimensions is discussed in Section 5.

Moreover, we focus on the case of constant coefficients as in (1.1) and the functional (1.5). While the numerical method itself is directly applicable to the variable coefficient case and more general functionals, the Fourier analysis we perform is tailored to the present setting.

2.1 Semi-implicit Milstein finite difference scheme

We use a spatial grid with uniform spacing $h_{x},\,h_{y}>0$, and let $V_{i,j}^{n}$ be the approximation to $v(nk,ih_{x},jh_{y})$, $n=1,\ldots,N$, $i,j\in\mathbb{Z}$. We assume for simplicity that $i_{0}:=x_{0}/h_{x}$ and $j_{0}:=y_{0}/h_{y}$ are integers. Then we approximate $v(0,x,y)=\delta(x-x_{0})\otimes\delta(y-y_{0})$ by

[TABLE]

We use the semi-implicit Milstein scheme to approximate (1.1), proposed in [19]:

[TABLE]

where $M$ is a random operator and

[TABLE]

and $\widetilde{Z}_{n}^{y}=\rho_{xy}Z_{n}^{x}+\sqrt{1-\rho_{xy}^{2}}Z_{n}^{y}$ with $Z_{n}^{x},Z_{n}^{y}\sim N(0,1)$ independent normals. Briefly, the terms on the left-hand side of (2.2) correspond to the implicit approximation of the operator $L$ in (1.1); the second and third terms on the right-hand side are the Euler-Maruyama approximation of the stochastic integral; and the last three terms the Milstein correction for strong first order 1. Note that the cross-derivative term on the left-hand side cancels out with a Milstein term, as detailed in [19].

To save computational cost, we combine the scheme with an Alternating Direction Implicit (ADI) factorisation [14],

[TABLE]

A detailed analysis of the $L_{2}$ stability and convergence of these schemes can be found in [19]. Here, we state the main result which is relevant here. We make the following technical assumptions throughout333They are used in the proof of Lemma 3.3.:

[TABLE]

Theorem 2.1** **(Corollary 2.1 in [19])

Let $\lambda>0$ be fixed and (2.4) be satisfied. Then there exists $C>0$ such that for any $k,h_{x},h_{y}$ with

[TABLE]

for the solution to the implicit Milstein scheme (2.2),

[TABLE]

This convergence result also holds for the ADI scheme (2.3).

2.2 Sparse combination estimators

Now we consider the specific functional

[TABLE]

as discussed in the introduction, where $v$ is the solution to (1.1) and (1.2).

By introducing integer multi-indices $\bm{l}=(l_{1},l_{2})\in\mathbb{N}^{2}$ as the refinement levels of the spatial mesh in dimensions $x$ and $y$, respectively, we denote by $P_{(l_{1},\,l_{2})}^{N}$ the discrete approximations to $P$ with mesh sizes $h_{x}=h_{0}\cdot 2^{-l_{1}}$, $h_{y}=h_{0}\cdot 2^{-l_{2}}$, and fixed timestep $k$ with $k\leq\lambda\min\{h_{x}^{2},h_{y}^{2}\}$. Let $N=T/k$. We then use the trapezoidal approximation

[TABLE]

where $V_{i,j}^{N}$ is the solution to (2.2).

Proposition 2.1

Assume (2.4) holds. Let $P$ be given by (2.6) and $P^{N}_{(l_{1},\,l_{2})}$ by (2.7). Then for fixed $\lambda>0$ there exists $C>0$ such that for any $(l_{1},l_{2})\in\mathbb{N}_{0}^{2}$, $k\leq\lambda\min\{h_{x}^{2},h_{y}^{2}\}$,

[TABLE]

where $h_{x}=h_{0}\cdot 2^{-l_{1}},\,h_{y}=h_{0}\cdot 2^{-l_{2}}$, $N=T/k$.

**Proof ** Similar to Proposition 2.2 in [18]. $\Box$

Let $\Delta_{i}$ be the first-order difference operator along directions $i=1,2$, defined as

[TABLE]

with $\bm{e}_{i}$ the canonical unit vectors in $\mathbb{R}^{2}$, i.e., $(\bm{e}_{i})_{j}=\delta_{ij}$, and $\bm{l}=(l_{1},l_{2})$. We also define the first-order mixed difference operator $\Delta=\Delta_{1}\otimes\Delta_{2}$. Hence, for $l_{1}>0,l_{2}>0$,

[TABLE]

We will prove in Section 3.3 the following theorem.

Theorem 2.2

Assume (2.4) holds. Let $h_{x}=h_{0}\cdot 2^{-l_{1}},\ h_{y}=h_{0}\cdot 2^{-l_{2}}$, and for any fixed $\lambda>0$, $k/\min\{h_{x}^{2},h_{y}^{2}\}\leq\lambda$, $N=T/k$. Then for $\Delta P^{N}_{(l_{1},\,l_{2})}$ from (2.10),

[TABLE]

**Proof ** See Section 3.3. $\Box$

Remark 2.1

We emphasise that (2.11) does not follow from (2.8), but ascertains essentially that the difference operator cancels out any terms in the error expansion which depend on $h_{x}$ or $h_{y}$ alone. Establishing this property, which depends on the regularity of the problem, is the crucial step for the application of the sparse combination technique in the deterministic case of PDEs (see [10, 3, 16]). It is precisely the condition required for the Multi-Index Monte Carlo method of [11]. Indeed, the approach here is interpretable as a specific multi-index decomposition applied to the spatial variables, treating time separately.

Corollary 2.1

With the same setting as in Theorem 2.2, the first order differences of $P^{N}_{(l_{1},l_{2})}$ derived from (2.7) have the first and second moments

[TABLE]

Following the ideas in [3, 11], given a sequence of index sets $\mathcal{I}_{l}=\{(l_{1},\,l_{2})\in\mathbb{N}_{0}^{2}:\,l_{1}+l_{2}\leq l+1\}$, the approximation on level $l$ is defined as

[TABLE]

Note that we use the same $k$ for all $\Delta P_{(l_{1},\,l_{2})}^{N}$, $(l_{1},\,l_{2})\in\mathcal{I}_{l}$. We have the following.

Theorem 2.3

Assume (2.4) holds. Let $k\leq\lambda\,h_{0}^{2}\,2^{-2l}$ for a fixed $\lambda>0$, $N=T/k$, $P$ given by (2.6) and $P_{l}^{N}$ given by (2.12). Then

[TABLE]

**Proof ** See Section 3.4. $\Box$

The standard Monte Carlo estimator for \mathbb{E}\big{[}P_{l}^{N}\big{]} using $M$ samples is defined as

[TABLE]

where \Big{(}\Delta\widehat{P}_{(l_{1},\,l_{2})}^{N}\Big{)}^{(m)} is the $m$-th sample for the difference on spatial grid level $(l_{1},\,l_{2})$ of the SPDE approximation using $N$ time steps.

To reduce the bias below $\varepsilon$, we can choose l=\big{[}\frac{1}{2}(-\log_{2}\varepsilon+\log_{2}|\log\varepsilon|)\big{]}, where $[\,\cdot\,]$ is the closest integer. Since $k\leq\lambda h_{0}^{2}2^{-2l}=O(2^{-2l})$, the computational cost for one sample of $P_{(l_{1},l_{2})}^{N}$ is

[TABLE]

Using standard Monte Carlo sampling, we need $M=O(\varepsilon^{-2})$ samples to reduce the variance below $\varepsilon^{2}$, hence the total computational cost to achieve a RMSE $\varepsilon$ is

[TABLE]

2.3 Sparse combination MLMC estimators

Next, we combine the sparse grid combination technique with the MLMC method. We will show that this combination leads to the optimal complexity $O(\varepsilon^{-2})$.

Let $P_{l}\coloneqq P_{l}^{N_{l}}$ be an approximation to $P$ as in (2.12) using a numerical discretisation with timestep $k_{l}$ and index set $\mathcal{I}_{l}$, where

[TABLE]

For $(l_{1},\,l_{2})\in\mathcal{I}_{l}$, the mesh sizes are $h_{x}^{(l_{1})}=h_{0}\cdot 2^{-l_{1}}$, $h_{y}^{(l_{2})}=h_{0}\cdot 2^{-l_{2}}$. By using different levels of refinement, we have the following identity:

[TABLE]

Then we define the difference operators (acting on the level)

[TABLE]

where we denote $P_{-1}\coloneqq 0$. Thus, the approximation $\mathbb{E}[P_{l^{\ast}}]$ to $\mathbb{E}[P]$ at level $l^{\ast}$ is

[TABLE]

Therefore, instead of directly simulating $P_{l^{\ast}}$, we simulate $\delta P_{l},\,l=0,1,\ldots,l^{\ast}$ separately. The key point is that we use the same Brownian path for $P_{l}$ and $P_{l-1}$ to calculate $\delta P_{l}=P_{l}-P_{l-1}$ such that the variance is considerably reduced. We have:

[TABLE]

where the first two inequalities follow from Theorem 2.3 and the third is immediate.

Let $\widehat{Y}_{l}$ be an estimator for \mathbb{E}\big{[}\delta P_{l}\big{]} using $M_{l}$ independent samples $\delta\widehat{P}_{l}^{(m)}$ of $\delta P_{l}$,

[TABLE]

The MLMC estimator is then defined as $\widehat{P}_{l^{\ast}}\coloneqq\sum_{l=0}^{l^{\ast}}\widehat{Y}_{l}$ or

[TABLE]

where \Big{(}\Delta\widehat{P}_{(l_{1},\,l_{2})}^{N_{l}}\Big{)}^{(m)} is the $m$-th sample for the difference operator on spatial grid level $(l_{1},\,l_{2})$ using $N_{l}$ time steps. Following [6, Theorem 3.1], choosing $M_{l}$ to minimise the computational cost for a fixed variance, using (2.16) and noticing that the variance decreases strictly faster than the cost increases (the polynomial terms in $l$ being negligible), we achieve the optimal complexity $O(\varepsilon^{-2})$.

3 Fourier analysis of the sparse combination error expansion

We will prove Theorem 2.2 and Theorem 2.3 in this section. We employ a Fourier transform and then analyse the different wave number regions separately, a technique that has been successfully used to derive error expansions for numerical approximations to PDEs (see [5]) and SPDEs (see [18]).

3.1 Fourier transform of the solution

Define the Fourier transform pair

[TABLE]

The Fourier transform of (1.1) yields

[TABLE]

subject to the initial data $\widetilde{v}(0)=\mathrm{e}^{-\mathrm{i}\xi x_{0}-\mathrm{i}\eta y_{0}}.$ For the remainder of the analysis, we take $\mu_{x}=\mu_{y}=0$. This does not alter the results (see Remark 2.3 in [17] for 1d).

The solution to (3.1) is

[TABLE]

where

[TABLE]

For the numerical solution, we can use a discrete-continuous Fourier pair

[TABLE]

where $i_{0}=x_{0}/h_{x}$, $j_{0}=y_{0}/h_{y}$. It follows from (2.1) that $\widetilde{V}^{0}(\xi,\eta)=1.$ Then we have

[TABLE]

where we make the ansatz

[TABLE]

It follows from (3.2) that the exact solution of $X(t_{n+1})$ given $X(t_{n})$ is

[TABLE]

where $W_{t_{n+1}}^{x}-W_{t_{n}}^{x}=:\sqrt{k}Z_{n,x},\ W_{t_{n+1}}^{y}-W_{t_{n}}^{y}=:\sqrt{k}\widetilde{Z}_{n,y}$ are the Brownian increments.

Now we consider the numerical approximation of $X_{n}$. Let

[TABLE]

and $e_{n}$ is the logarithmic error between the numerical solution and the exact solution introduced during $[nk,(n+1)k]$. Aggregating over $N$ time steps, at $t_{N}=kN=T$,

[TABLE]

where

[TABLE]

is the exact solution at time $T$. Moreover, inserting $V_{i,j}^{n}$ from (2.2) to (3.3), we have

[TABLE]

where

[TABLE]

Hence, $C_{n}$ can also be expressed as

[TABLE]

3.2 Fourier transform of the sparse combination estimator

With the notation from Section 2.2, omitting $N$ in $P$ to simplify the notation,

[TABLE]

where $\chi(h_{x},\xi),\ \chi(h_{y},\eta)$ is defined in a distributional sense as

[TABLE]

Note that $\chi$ only appears multiplied by the smooth, fast decaying function $X_{N}$ and in integral form, such that this is well-defined.

We recall from (2.10) the sparse combination estimator

[TABLE]

We assume $h_{x},\ h_{y}<1$. Even though $P_{(l_{1},\,l_{2})},\ P_{(l_{1},\,l_{2}-1)},\ P_{(l_{1}-1,\,l_{2})},\ P_{(l_{1}-1,\,l_{2}-1)}$ have different Fourier domains, we define $\Omega_{\text{low}}$ shared by all of them, for $0<p<1/2$,

[TABLE]

Then we define $I_{(l_{1},\,l_{2})}$ as

[TABLE]

where $\Omega_{(l_{1},\,l_{2})}=[-\pi/h_{x},\pi/h_{x}]\times[-\pi/h_{y},\pi/h_{y}]$. Then

[TABLE]

and

[TABLE]

Denote

[TABLE]

Then we have

[TABLE]

and

[TABLE]

3.3 Proof of Theorem 2.2

In this section, we give a proof of Theorem 2.2. The splitting into different wave number regions is motivated by the analysis of the heat equation in [5]. We further separate $\Omega_{(l_{1},\,l_{2})}\setminus\Omega_{\text{low}}=\Omega_{\text{mid}}\cup\Omega_{\text{high}}$. Hence $\Omega_{(l_{1},\,l_{2})}$ is divided into three regions,

[TABLE]

We state three lemmas without proof, before giving the proof of Theorem 2.2.

Lemma 3.1** **(Low wave region)

For $G_{(l_{1},\,l_{2})}(\xi,\eta)$ introduced in (3.10), there exists a constant $C>0$, such that for any $l_{1},l_{2}\in\mathbb{N}$,

[TABLE]

where $R(T)$ is a random variable satisfying $|\mathbb{E}[R(T)]|\leq C,\;\mathbb{E}[|R(T)|^{2}]\leq C.$

**Proof ** See Appendix A.1. $\Box$

Lemma 3.2** **(Middle wave region)

For the middle wave region, we have,

[TABLE]

**Proof ** See Appendix A.2. $\Box$

Lemma 3.3** **(High wave region)

For the high wave region, we have

[TABLE]

**Proof ** See Appendix A.3. $\Box$

**Proof **[Theorem 2.2] We have derived in Section 3.2 that

[TABLE]

Lemma 3.2 and Lemma 3.3 give that

[TABLE]

for all $r>0$. Similar estimates hold for $I_{(l_{1}-1,\,l_{2})}$, $I_{(l_{1},\,l_{2}-1)}$, and $I_{(l_{1}-1,\,l_{2}-1)}$. Combining this with Lemma 3.1, the two inequalities in (2.11) follow. $\Box$

3.4 Proof of Theorem 2.3

**Proof **[Theorem 2.3] Since from Proposition 2.1,

[TABLE]

we have

[TABLE]

$\Box$

4 Numerical tests

In this section, we test the theoretical convergence results empirically.

4.1 Mean and variance of hierarchical increments

First, we verify numerically for $l_{1},\,l_{2}>0$, $h_{x}=h_{0}\cdot 2^{-l_{1}},\,h_{y}=h_{0}\cdot 2^{-l_{2}}$, the result from Theorem 2.2 that

[TABLE]

We choose parameters $h_{0}=1,\ T=1,\ x_{0}=y_{0}=2,\ \mu_{x}=\mu_{y}=0.0809,\ \rho_{x}=\rho_{y}=0.2,\ \rho_{xy}=0.45$, and truncate the domain to $(x,y)\in[-8,12]\times[-8,12]$.

Table 1 shows \log_{2}\big{|}\mathbb{E}[\Delta P^{N}_{(l_{1},\,l_{2})}]\big{|} and $\log_{2}\mathrm{Var}[\Delta P^{N}_{(l_{1}\,l_{2})}]$ with fixed timestep $k=4^{-3}$, and different levels of mesh refinement. We can see from the table that the mean decreases by around 2 going from level $l$ to $l+1$, and the variance decreases by approximately 4, consistent with the theoretical prediction. Figure 1 depicts the corresponding contour plots.

4.2 Mean and variance of sparse grid increments

Next, we estimate \mathbb{E}\big{[}\delta P_{l}\big{]} and \mathrm{Var}\big{[}\delta P_{l}\big{]}, where $\delta P_{l}$ is given by (2.15). From Table 1 we deduce that for the index set $I_{l}=\{(l_{1},l_{2}):l_{1}+l_{2}\leq l+1\}$, the terms on the ‘boundary’ (i.e., $l_{1}=0$ or $l_{2}=0$) will dominate. Although this does not affect the total order of complexity, we can further optimise the cost by modifying the index set such that the contribution from the boundary and ‘interior’ are similar.

Define therefore $\widehat{I}_{l}$ as the indices for interior meshes, and $\widetilde{I}_{l}$ for boundary meshes,

[TABLE]

We balance the contributions from these two sets by finding, for some fixed $N$,

[TABLE]

and using the index set

[TABLE]

For example, from Table 1, $l^{\ast}=2$. So we have

[TABLE]

We use $k_{l}=k_{0}\cdot 2^{-2l}$ and $\mathcal{I}_{l}$ for the construction of $\delta P_{l}$. Table 2 verifies

[TABLE]

Figure 2 are the corresponding plots. The fitted slopes in Figure 2 are $-1.91,\ -3.73$, and $3.27$, respectively, compared to the theoretical asymptotic values of $-2$, $-4$, and $3$ (neglecting logarithmic terms, which will play a role for low levels).

The cost is counted here as the total number of operations. Specifically, under the ADI scheme (2.3), the cost of one numerical realisation of (1.1) with mesh size $h_{x},\,h_{y}$ and timestep $k$ is

[TABLE]

Since the mean square error for the estimator can be expressed as sum of the variance and the square of the weak error, we split the accuracy “budget” as

[TABLE]

and optimize over $\alpha$. Since $V_{l}C_{l}\sim l^{3}\,2^{-l},$ the variance decays with levels more rapidly than the cost increases with levels, and thus the dominant cost is on level [math]. Therefore, we choose $\alpha$ relatively small to reduce the cost on level [math], and hence we reduce the total cost. To find the optimal $\alpha$, one approach is to approximate the total cost given different $\alpha$, and choose the one which minimise the complexity.

Figure 3a is the loglog plot of total cost among all the methods mentioned before: standard Monte Carlo, sparse combination, multilevel Monte Carlo, sparse combination with MLMC. All the schemes use the optimal $\alpha$ for each accuracy. We can see from the graph that standard MC gives the cost $O(\varepsilon^{-4})$, and sparse combination gives $O(\varepsilon^{-7/2})$, as expected. As for MLMC and sparse combination with MLMC methods, both yield approximately $O(\varepsilon^{-2})$, which verifies our proof as the log term in MLMC cost is negligible in this plot.

Figure 3b compares $\varepsilon^{2}W$ between sparse combination with MLMC and MLMC alone (without sparse combination). The total computational cost of sparse combination with MLMC is approximately proportional to $\varepsilon^{-2}$, hence $\varepsilon^{2}W$ does not vary significantly for different accuracy $\varepsilon$. However, as the total cost of the multilevel scheme is proportional to $\varepsilon^{-2}(\log\varepsilon)^{2}$, and we can see from the figure that $\varepsilon^{2}W$ increases as $\varepsilon$ goes to zero. For the black line, we use $\alpha=0.1$. For the blue dotted line, each scheme uses the optimal $\alpha$ for different accuracies. As a result, sparse combination with MLMC achieves the optimal order of complexity.

5 Generalisation to higher dimensions

We can generalise the SPDE to $d$ dimensions, where $d\geq 3$. Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, where there is given a $d$-dimensional standard Brownian motion $W$ with correlation matrix $\Sigma=(\rho_{k,m})$. The natural extension of the SPDE (1.1) is

[TABLE]

for $x\in\mathbb{R}^{d},\ 0<t\leq T$, where $\mu\in\mathbb{R}^{d}$ and $\rho\in(0,1)^{d}$ are parameters, subject to the Dirac initial datum $v(0,x)=\delta(x_{1}-x_{1,0})\otimes\cdots\otimes\delta(x_{d}-x_{d,0})$, where $x_{0}\in\mathbb{R}_{+}^{d}$ is given.

We use a spatial grid with uniform spacing $h_{x_{i}}>0$, $i=1,\ldots,d$, and let $V_{i_{1},\ldots,i_{d}}^{m}$ be the approximation to $v(mk,i_{1}h_{x_{1}},\ldots,i_{d}h_{x_{d}})$, $m=1,\ldots,N$, $i_{j}\in\mathbb{Z}$, where $i_{0,j}\coloneqq[x_{0,j}/h_{x_{j}}]$, the closest integers to $x_{0,j}/h_{x_{j}}$. We approximate $v(0,x_{1},\ldots,x_{d})$ by

[TABLE]

In analogy with (1.5), we define a linear functional by

[TABLE]

Similar as before, we use an implicit finite difference scheme to approximate (5.1), and we use the trapezoidal rule for $P$ with a truncation of the domain. We have the following conjectures and results.

Conjecture 5.1

Assume the implicit finite difference scheme is stable 444We expect this to be true for ‘small enough’ correlations as in the two-dimensional case for (2.4), but it is not obvious what the conditions will be for specific $d>2$ without performing the analysis., and the timestep $k$ and mesh size $h_{x_{i}}$ satisfy

[TABLE]

for arbitrary fixed $\lambda>0$. Then we have the error expansion

[TABLE]

where $N=T/k$ is the number of time steps.

Next we apply the sparse combination method to (5.2) with fixed $k$ satisfying (5.3). Similar to Section 2.2, let $\Delta_{i}$ be the first-order difference operator along directions $i=1,\ldots,d$, defined as in (2.9), and $\Delta=\Delta_{1}\otimes\cdots\otimes\Delta_{d}$.

Conjecture 5.2

Assume the implicit finite difference scheme is stable. Let $h_{x_{i}}=h_{0}\cdot 2^{-l_{i}}$, $i=1,\cdots,d$, and $k$ be the timestep such that for an arbitrary fixed $\lambda>0$, $k\leq\lambda\min\{h_{x_{1}}^{2},h_{x_{2}}^{2},\cdots,h_{x_{d}}^{2}\}$, where $N=T/k$ the number of time steps. Then the first and second moments of $\Delta P^{N}_{\bm{l}}$ satisfy

[TABLE]

Given a sequence of index sets

[TABLE]

the approximation on level $l$ is defined (similar to (2.12)) as

[TABLE]

Note that we use the same $k$ for all $\Delta P^{N}_{\bm{l}}$, $(l_{1},\cdots,l_{d})\in\mathcal{I}_{l}$. Then we have:

Proposition 5.1

Assume Conjectures 5.1 and 5.2 to be true. Then, for $P$ given by (5.2) and $P_{l}^{N}$ by (5.5), we have

[TABLE]

Moreover, by choosing l=\big{[}\frac{1}{2}(-\log_{2}\varepsilon+(d-1)\log_{2}|\log\varepsilon|)\big{]}, the computational cost $W$ to achieve a RMSE $\varepsilon$ using standard Monte Carlo estimation is

[TABLE]

Next, we combine MLMC with the sparse combination scheme. Similar to Section 2.3, let $P_{l}\coloneqq P_{l}^{N_{l}}$ be an approximation to $P$ as in (2.12) using a discretisation with timestep $k_{l}$ and index set $\mathcal{I}_{l}$ defined above. Then we define

[TABLE]

where we denote $P_{-1}\coloneqq 0$. Thus the approximation to $P$ at level $l^{\ast}$ has the form

[TABLE]

In this way, we simulate $\delta P_{l},\,l=0,1,\ldots,l^{\ast}$ instead of directly simulating $P_{l^{\ast}}$, such that the variance of $\delta P_{l}=P_{l}-P_{l-1}$ is considerably reduced by using the same Brownian path for $P_{l}$ and $P_{l-1}$.

Let $\widehat{Y}_{l}$ be an estimator for \mathbb{E}\big{[}\delta P_{l}\big{]} using $M_{l}$ samples. Each estimator is an average of $M_{l}$ independent samples, where $\delta\widehat{P}_{l}^{(m)}$ is the $m$-th sample arising from a single SPDE approximation,

[TABLE]

The MLMC estimator is defined as

[TABLE]

where \Big{(}\Delta\widehat{P}_{\bm{l}}^{N_{l}}\Big{)}^{(m)} is the $m$-th sample for the difference on spatial grid level $\bm{l}=(l_{1},\ldots,l_{d})$ of the SPDE approximation using $N_{l}$ time steps. Following [6], through optimising $M_{l}$ to minimise the computational cost for a fixed variance, we can achieve the optimal complexity $O(\varepsilon^{-2})$ in this case.

Proposition 5.2

Assume Conjectures 5.1 and 5.2 to be true. Given $\delta P_{l}$ from (5.7), there exist $C_{1},\,C_{2},\,C_{3}>0$, such that

[TABLE]

Then, given a RMSE $\varepsilon$, the MLMC estimator (5.8) leads to a total complexity $O(\varepsilon^{-2})$.

**Proof ** The first inequalities follow directly from (5.4). The complexity then is a small modification of [6, Theorem 3.1]. $\Box$

If we use MLMC without sparse combination to estimate $\mathbb{E}[P_{t}]$, by letting $h_{x_{j}}=h_{0}\cdot 2^{-l}$, $k=k_{0}\cdot 2^{-2l}$, then for constants $C_{1},C_{2},C_{3}>0$, independent of $h$ and $k$,

[TABLE]

Similarly, we have the following result.

Proposition 5.3

Given a RMSE $\varepsilon$, the total computational cost $W$ for $P$ in (5.2) using MLMC satisfies

[TABLE]

Comparing (5.9) with (5.6), we can see that when $d>5$, the sparse combination scheme with standard MC performs better than MLMC on regular grids.

6 Conclusion

We considered a two-dimensional parabolic SPDE and a functional of the solution. We analysed the accuracy and complexity of a sparse combination multilevel Monte Carlo estimator, and showed that, by using a semi-implicit Milstein finite difference discretisation (2.2), we achieved the order of complexity $O(\varepsilon^{-2})$ for a RMSE $\varepsilon$, whereas the cost using standard Monte Carlo is $O(\varepsilon^{-4})$, and that using MLMC is $O(\varepsilon^{-2}(\log\varepsilon)^{2})$. When generalising to higher-dimensional problems, sparse combination with MLMC maintains the optimal complexity, whereas MLMC has an increasing total cost as the dimension increases.

Further research will apply this method to a Zakai type SPDE with non-constant coefficients. Another open question is a complete analysis of the numerical approximation of initial-boundary value problems for the considered SPDE.

Appendix A Proofs of some auxiliary stability results

A.1 Proof of Lemma 3.1

**Proof ** For $(\xi,\eta)\in\Omega_{\text{low}}$, we have $G_{(l_{1},\,l_{2})}(\xi,\eta)$ as in (3.11). From (3.5), we have

[TABLE]

where

[TABLE]

By a similar proof to Lemma 4.1 in [19], we have

[TABLE]

Here, $R_{0},\ R_{1}$, $R_{2}$, $R_{3}$ are random variables satisfying

[TABLE]

where $f_{i}(\cdot)$ are polynomial functions. So in the low wave region, we get

[TABLE]

Here,

[TABLE]

is the numerical approximation to

[TABLE]

by the trapezoidal rule with error $O(h_{x}^{2})+O(h_{y}^{2})$. Therefore,

[TABLE]

By the same reason,

[TABLE]

Therefore, we get

[TABLE]

Similarly,

[TABLE]

$\Box$

A.2 Proof of Lemma 3.2

**Proof ** As

[TABLE]

it is enough to show for $\Omega_{\text{mid}}^{1}=\{0<|\xi|<h_{x}^{-p},\ h_{y}^{-p}<|\eta|<\pi/h_{y}\}$,

[TABLE]

We have

[TABLE]

where $X(T)$ is introduced in (3.2), and $X_{N}$ is given by (3.6). Specifically,

[TABLE]

As $X(T)$ is given in closed form, a direct calculation gives

[TABLE]

Hence, in the following we focus on

[TABLE]

Since $h_{x}h_{y}\chi(h_{x},\xi)\chi(h_{y},\eta)$ has finite 2-norm, it is justifiable to consider

[TABLE]

First we denote

[TABLE]

so that

[TABLE]

Then we introduce $\widetilde{X}_{n}$ such that

[TABLE]

where

[TABLE]

It has been proved in [18] that

[TABLE]

[TABLE]

Thus there exists a constant $C$ such that

[TABLE]

and as a result,

[TABLE]

Since

[TABLE]

it follows that

[TABLE]

Therefore,

[TABLE]

$\Box$

A.3 Proof of Lemma 3.3

**Proof ** For the same reason as in the proof of Lemma 3.2, it is sufficient to prove

[TABLE]

Note that

[TABLE]

By (3.6), we have

[TABLE]

We denote

[TABLE]

In this case,

[TABLE]

where

[TABLE]

As we have $k<\lambda\min\{h_{x}^{2},h_{y}^{2}\}$,

[TABLE]

Also,

[TABLE]

We define a constant

[TABLE]

then

[TABLE]

Therefore we have

[TABLE]

$\Box$