Stability and error analysis of an implicit Milstein finite difference scheme for a two-dimensional Zakai SPDE

In this article, we propose an implicit finite difference scheme for a two-dimensional parabolic stochastic partial differential equation (SPDE) of Zakai type. The scheme is based on a Milstein approximation to the stochastic integral and an alternating direction implicit (ADI) discretisation of the elliptic term. We prove its mean-square stability and convergence in L2 of first order in time and second order in space, by Fourier analysis, in the presence of Dirac initial data. Numerical tests confirm these findings empirically.


Introduction
The analysis and numerical computation of Zakai equations and other types of stochastic partial differential equation (SPDE) have been extensively studied in recent years. A general form of Zakai equation (see [BC09,GPPP06]) is given by where M is an m-dimensional standard Brownian motion, a is a d × d matrix-valued function, b is a R d -valued function, and γ is a d × m matrix-valued function. This Zakai equation arises from a nonlinear filtering problem: given an m-dimensional observation process M and a d-dimensional signal process Z, the goal is to estimate the conditional distribution of Z given M . If Z satisfies where B is a d-dimensional standard Brownian motion independent of M , σ is a d × d-matrix valued function, and β is a R d -valued function, then the conditional distribution function of Z given M has a density v in L 2 , and it is proved (Theorem 3.1 in [KX99]) that under appropriate conditions, v satisfies (1.1) in a weak sense with Moreover, the solution v to (1.1) can be interpreted as the density -if it exists -of the limit empirical where B i and N are independent Brownian motions, independent of M , and the rest as above.
There are two major approaches to the numerical approximation of the Zakai equation. One is by simulating the particle system (1.3) with a Monte Carlo method, for instance as in [CL99,Cri03,CX10,GPPP06]. The other approach is to directly solve the Zakai SPDE by spatial approximation methods and time stepping schemes, coupled again with Monte Carlo sampling, which is the subject of this paper.
For this second class of methods, several schemes have been developed in earlier works including finite differences [GN97,Gyö98,Gyö99,DG01], finite elements [Wal05,Kru14], and stochastic Taylor schemes [JK09,JK10], but these were restricted to classes of SPDEs not including Zakai equations of the type (1.1).
More recently, methods have been developed and analysed for parabolic SPDEs of the generic form dv = Lv dt + G(v) dM t , (1.4) where L is a second order elliptic differential operator, and G is a functional mapping v onto a linear operator from martingales M into a suitable function space.
Under suitable regularity, for equations of type (1.4), mean-square convergence of order 1/2 is shown for an Euler semi-discretisation in [Lan10] for square-integrable (not necessarily continuous), infinitedimensional martingale drivers. In contrast, [BL13] allow only for continuous martingales but prove convergence of higher order in space and up to 1 in time, in L p and almost surely, for a Milstein scheme and spatial Galerkin approximation of sufficiently high order; this is extended to advection-diffusion equations with possibly discontinuous martingales in [BL12].
Giles and Reisinger [GR12] use an explicit Milstein finite difference approximation to the solution of the following one-dimensional SPDE, a special case of (1.1) for d = 1 and constant coefficients, where T > 0, M is a standard Brownian motion, and µ and 0 ≤ ρ < 1 are real-valued parameters. This is extended in [Rei12] to an approximation of (1.5) with an implicit method on the basis of the σ-θ time-stepping scheme, where the finite variation parts of the double stochastic integral are taken implicit. This is further applied in [RW18] to Multi-index Monte Carlo estimation of expectations of a functional of the solution.
A finite difference scheme for a filtered jump-diffusion process resulting in a stochastic integrodifferential equation is studied in [DL16], where convergence of order 1 in space and 1/2 in time, in L 2 and L ∞ in space, is proven for an Euler time stepping scheme.
The theoretical results in this paper are an extension from those in [GR12] to the multi-dimensional case. They are more specific than those in [DL16] in that we analyse only the case of constant coefficient local SPDEs. In contrast to [BL12,BLS13], we consider only finite-dimensional Brownian motions, as is relevant in our applications. But we specifically include the case of Dirac initial data and extend the results to a practically attractive, semi-implicit alternating direction implicit factorisation in the context of the Milstein scheme.
We want to allow for Dirac initial data because they correspond to the natural situation where all particles in (1.3) start from the same initial position, or a filtering problem with known current state Z 0 in (1.2).
A classical result states that, for a class of SPDEs including (1.6), with initial condition in H 0 , there exists a unique solution v ∈ L 2 (Ω × (0, T ), F, H 0 (R)) [KR81]. This does not include Dirac initial data (1.7), but in fact, the solution to (1.6) and (1.7) is analytically known to be the smooth (in x and y) function (1.8) The availability of a closed-form solution in this case helps us check the validity of our numerical scheme and its convergence rate, although the scheme itself is more widely applicable.
For the SPDE (1.6), we consider both explicit and implicit Milstein schemes. We study the meansquare stability, and the strong convergence of the second moment. This can give us an error bound for the expected error. The advantage over the simpler Euler scheme is that the strong convergence order is improved from 1/2 to 1. As expected, we find that the explicit scheme is stable in the mean-square sense only under a strong CFL-type condition on the timestep, k ≤ Ch 2 for timestep k, mesh size h, and a constant C, while the implicit scheme is mean-square stable under the very mild and somewhat unusual CFL condition k ≤ C| log h| −1 (provided also some constraints on ρ x , ρ y , ρ xy ).
We therefore focus on the implicit scheme, for which we prove first order convergence in the timestep and second order in the spatial mesh size. The analysis is made more difficult by the Dirac initial datum compared to, say, L 2 initial data. We adapt the approach used in [CG07] for the heat equation by studying the convergence for different wave number regions in Fourier space and then assemble the contributions to the error by the inverse transform.
Furthermore, we use an Alternating Direction Implicit (ADI) scheme to approximately factorise the discretisation matrix for the implicit elliptic part in (1.6). This concept is well established for PDEs (see, e.g., [PR55,CS88,HV13]). It is well known that in the multi-dimensional case standard implicit schemes result in sparse banded linear systems, which cannot be solved by direct elimination in a computational cost which scales linearly with the number of unknowns like in the one-dimensional, tridiagonal case. An alternative to advanced iterative linear solvers such as multigrid methods is to reduce the large sparse linear system approximately to a sequence of tridiagonal linear systems, which are computationally easier to handle, by ADI factorisation. To our knowledge, the present work is the first application of ADI to SPDEs. We show that the ADI approximation is also mean-square stable under the same conditions as the original implicit scheme and has the same convergence order.
We note that published analysis of ADI schemes for parabolic PDEs in the presence of mixed spatial derivative terms is currently restricted to constant coefficients (through the use of von Neumann stability analysis; see e.g. [WitH16]). Notwithstanding this, the empirical evidence overwhelmingly suggests that the conclusions drawn there extend to most cases of variable coefficients.
The scheme we propose applies similarly beyond constant coefficients. We give the natural extension to the SPDE (1.1). In that case, additional iterated stochastic integrals (the Lévy area) appear in the Milstein approximation. The efficient, accurate simulation has been studied in the context of SDEs in [KPW92,GL94] and invariably leads to relatively complicated schemes. As the computational effort in the context of the SPDE (1.1) is dominated by the matrix computations from the finite difference scheme, we perform a simple approximation of the stochastic integrals t+k t (W s − W t ) dB s , for correlated Brownian motions W and B, by simple Euler integration with step k 2 . This is sufficiently accurate not to spoil the first order convergence and does not increase the complexity order.
As a specific application, we approximate the equation (1.9) taken from [HK17], with the scheme presented in this paper. Although our analysis (based on Fourier transforms) does not directly apply in this case, the scheme preserves first order convergence in time and second order convergence in space in our numerical tests.
Summarising, the novel contributions of this paper are as follows. We • give a rigorous stability and error analysis in L 2 for a Milstein finite difference scheme for the SPDE (1.6), deriving sharp leading order error terms; • derive pointwise errors for Dirac initial data, which reveal a mild instability for large implicit timesteps and small spatial mesh sizes in this case, not seen in previous studies for L 2 data; • extend the analysis to an alternating direction implicit (ADI) factorisation, which, to our knowledge, is the first application of an ADI scheme to stochastic PDEs; • propose a modification for the more general equation (1.1) through sub-simulation of the Lévy area, which is empirically shown to be of first order.
The rest of this article is structured as follows. We define the approximation schemes in Section 2. Then we analyse the mean-square stability and L 2 -convergence in Sections 3 and 4 in the constant coefficient case of (1.6). Section 5 shows numerical experiments confirming the above findings. Section 6 extends the scheme to variable coefficients as in (1.1) and presents tests for the example (1.9). Section 7 offers conclusions and directions for further research.
2 Approximation and main results

Semi-implicit Milstein finite difference scheme
First, we introduce the numerical scheme to the SPDE (1.6), repeated here for convenience, with Dirac initial v(0, x, y) = δ(x−x 0 )δ(y−y 0 ). We use a spatial grid with uniform spacing h x , h y > 0, and, for T > 0 fixed, N time steps of size k = T /N . Let V n i,j be the approximation to v(nk, ih x , jh y ), n = 1, . . . , N , i, j ∈ Z, where i 0 := [x 0 /h x ], j 0 := [y 0 /h y ], the closest integers to x 0 /h x and y 0 /h y . We approximate v(0, x, y) by To improve the accuracy of the approximation of v in the present case of Dirac initial data, we subsequently choose h x and h y such that x 0 /h x and y 0 /h y are integers and therefore x 0 and y 0 are on the grid.
Extending the implicit Euler scheme in [Rei12] for the 1D case, it is natural to take the drift term implicit, and the terms driven by M x and M y explicit. We will prove later that in this way we obtain better stability (compare Proposition 3.2 to Theorem 3.1). For computational simplicity, in the following, we take the mixed derivative term therein explicit. This is also in preparation for the ADI splitting schemes we will study later.
Using such a semi-implicit Euler scheme, the system of SPDE (1.6) can be approximated by and Z y n = ρ xy Z x n + 1 − ρ 2 xy Z y n , with Z x n , Z y n ∼ N (0, 1) being independent normal random variables. To achieve a higher order of convergence, we introduce the Milstein scheme. Integrating (1.6) over the time interval [nk, (n + 1)k], In the Euler schemes, we approximate all integrands by their value at time nk or (n + 1)k, which is a zero-order expansion in time. By contrast, in the Milstein scheme, we use a first-order expansion for the stochastic integrals, such that we make the approximation v(s, x, y) ≈ v(nk, x, y) for nk < s < (n + 1)k in the first integral and in the second and third. We denote nk as t, and it follows We see that the mixed-derivative terms cancel, and we derive the implicit Milstein scheme as follows, To facilitate its implementation, we combine the scheme with an Alternating Direction Implicit (ADI) factorisation, which has been introduced in [PR55] for parabolic PDEs to approximately factorise the system matrix by matrices which correspond to derivatives in individual directions and which can thus more easily be inverted, while the consistency order is maintained. Applying this principle to the non-Brownian, implicit terms of the SPDE, we obtain Note that there is no substantial benefit in considering second order accurate splitting schemes (such as Craig-Sneyd [CS88] or Hundsdorfer-Verwer [HV13]) as the overall order is limited to 1 by the Milstein approximation to the stochastic integral.
We approximate the second derivative on the right-hand side with D 2 x and D 2 y , but the results for D xx and D yy would be similar.
We can also use the explicit Milstein finite difference scheme to approximate (2.5) but, as we will see, this scheme is stable only under a restrictive condition on the timestep.

Main convergence results
The following theorems describe the mean-square stability and convergence of the implicit finite difference scheme (2.3) and the ADI scheme (2.4). We make the following assumption: Section 3 shows that Assumption (2.1) is a sufficient condition for stability of the schemes (2.3) and (2.4). 1 If ρ xy = 0, these conditions reduce to 2ρ 2 x ≤ 1 and 2ρ 2 y ≤ 1, which is analogous to the condition for mean-square stability in the 1-dimensional case in [Rei12]. In the worst case, |ρ xy | = 1, sufficient conditions are ρ x , ρ y ≤ 1/ √ 6, and ρ x ρ y ≤ 1/12.
First, we recall the setting of our numerical schemes. For T > 0 fixed, we discretise [0, T ] by N steps, and let k = T /N be the timestep. We discretise space R 2 with mesh sizes h x and h y .
The following theorem shows the convergence of the implicit Milstein difference scheme (2.3). The constant θ ∈ (0, 1) therein is determined by the parameters ρ x , ρ y and ρ xy . The proof of Lemma 4.5 gives an explicit value; though not sharp, this is sufficient to highlight the divergence for h x , h y → 0 when k is fixed.
and h y > 0 be mesh sizes. Then, under Assumption 2.1, there exists θ ∈ (0, 1), independent of h x , h y and k, such that the implicit Milstein finite difference scheme (2.3) has the error expansion . , E 5 , and R are random variables with bounded first and second moments, all independent of h x , h y and k.
Proof See Section 4. ✷ for some β, C 0 > 0 independent of h x , h y and k, or, equivalently, then the implicit Milstein scheme (2.3) has the error expansion For the ADI discretisation scheme (2.4), a similar convergence result holds.
Theorem 2.2 Under the conditions of Remark 2.1, the error of the ADI scheme (2.4) has the same order as for the implicit Milstein scheme, where E 1 , E 2 , E 3 and R are random variables with bounded first and second moments.
Proof See Section 4.
✷ Theorem 2.1 and Theorem 2.2 state the convergence pointwise in space and L 2 in probability. If we consider L 2 convergence in space, then by applying Parseval's theorem, we can get further results.
Corollary 2.2 Under the conditions of Theorem 2.1, the L 2 error in space and in probability of the implicit Milstein scheme (2.3) at time T satisfies, If the initial condition lies in L 2 , then Proof See Section 4. ✷

Fourier analysis of mean-square stability
Recall the SPDE (1.6), The Fourier transform of (3.1) yields subject to the initial data v(0) = e −iξx 0 −iηy 0 . For the remainder of the analysis, we take µ x = µ y = 0. This does not alter the results (see Remark 2.3 in [Rei12] for the 1d case).
For the numerical solution, we can use a discrete-continuous Fourier decomposition (note that we approximate v(nk, ih x , jh y ) by V n i,j ) (3.4) In the last step, we integrate by substitution, By analogy with the theoretical solution v(t) = X(t) v(0), we make the ansatz . We can regard X n (ξ, η) as the numerical approximation to X(nk) in (3.3).
We say that the scheme is asymptotically mean-square stable, provided for any (ξ, η) This concept has been defined in the context of systems of SDEs in Definition 2.2, 3., in [BK10], and we apply it here to a fixed wave number in the Fourier domain. A generalisation to SPDEs is analysed in [LPT17] (see Definition 2.1 therein). A link between (3.6) and mean-square stability of the SPDE discretisation can be established using Parseval's equality, if it can be shown that the 2-norm of X n diminishes. We will not do this here but show convergence in L 2 (for fixed T) directly under the same conditions; see [Rei12] for mean-square stability and convergence of a 1-d parabolic SPDE. If (3.6) holds without any restriction between h x , h y and k, we call it unconditionally stable. This leads to three conditions summarised in Assumption 2.1, as shown by the following Theorem 3.1.
Theorem 3.1 The implicit Milstein finite difference scheme (2.3) is unconditionally stable in the mean-square sense of (3.6) provided Assumption 2.1 holds.
Proof By inserting (3.4) and (3.5) in (2.3), we have where To ensure mean-square stability, it is necessary and sufficient (given the multiplicative form and time-homogeneity of (3.7)) that for any (ξ, η) i.e., This is equivalent to (3.10) Note that Z y n = ρ xy Z x n + 1 − ρ 2 xy Z y n , with Z x n , Z y n ∼ N (0, 1) being independent normal random variables, hence E Z 2 n,x Z 2 n,y = 1 + 2ρ 2 xy , E Z n,x Z n,y = ρ xy , E Z 3 n,x Z n,y = E Z n,x Z 3 n,y = 3ρ xy , and E c x ρ x kZ n,x + c y ρ y k Z n,y 2 = c 2 x ρ x k + c 2 y ρ y k + 2c x c y √ ρ x ρ y ρ xy k.
✷ Now we prove the stability of the ADI scheme.
Proof By insertion in (2.4), we have for a x , a y ≤ 0, the stability also holds for the ADI scheme. ✷ Proposition 3.2 The explicit Milstein (finite difference) scheme (2.5) is stable in the mean-square sense provided Proof To ensure L 2 stability in this case, we need E 1+(a x +a y )k−ic x ρ x kZ n,x −ic y ρ y k Z n,y +b x ρ x k(Z 2 n,x −1)+b y ρ y k( Z 2 n,y −1)+d √ ρ x ρ y kZ n,x Z n,y 2 < 1.
For simplicity, we denote u = | sin ξhx 2 |, v = | sin ηhy 2 |, then we have This leads to the two sufficient conditions in (3.12). ✷ It follows that if ρ xy = 0, the stability conditions are So it is sufficient that k/h 2 x ≤ 1/6, and k/h 2 y ≤ 1/6. If |ρ xy | = 1, which is the worst case in (3.12), the stability conditions are So it is sufficient to ensure k/h 2 x ≤ 1/24, and k/h 2 y ≤ 1/24.

Fourier analysis of L 2 -convergence
Extending the analysis in [CG07] for the standard 1D (deterministic) heat equation to our 2D SPDE, we compare the numerical solution to the exact solution in Fourier space first by splitting the Fourier domain into two wave number regions. Assume p is a constant satisfying 0 < p < 1 4 . Then we define the low wave number region by and the high wave number region by Note that both X n and X(nk) are functions of ξ and η. The idea of the convergence proof is that X n is a good approximation to X(nk) in the low wave region, and they both damp exponentially in the high wave region.
Proof See Section 4.1. ✷ Lemma 4.2 Under Assumption 2.1, there exists C > 0 independent of h x , h y , and k, such that where 0 < θ < 1 is independent of h x , h y , and k.
Proof See Section 4.2.

✷
The following theorem shows mean square convergence of the implicit finite difference scheme (2.3).
Proof [Theorem 2.1] By Lemma 4.1 and Lemma 4.2, the inverse Fourier transform gives where x i = ih x , y j = jh y , E 1 , . . . , E 5 , and R are random variables with bounded first and second moments, N = T /k, 0 < θ < 1, and θ is independent of h x , h y and k. ✷ Next we give a proof of Corollary 2.2, the L 2 convergence in space and probability of the implicit finite difference scheme (2.3).

Low wave number region (proof of Lemma 4.1)
For the low wave region, we consider the case where both ξ, η are small. It follows from (3.3) that the exact solution of X(t n+1 ) given X(t n ) is where M x t n+1 − M x tn ≡ √ kZ n,x , M y t n+1 − M y tn ≡ √ k Z n,y are the Brownian increments.
Now we consider X n , the numerical approximation of X(nk). Let 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 , where is the exact solution at time T .
Remark 4.1 We can derive the exact leading order term by taking the inverse Fourier transform. For instance, the leading order error in h x is and similar for h y (replacing 'x' by 'y'); the leading order error in k can be found by the same technique but is significantly lengthier and hence omitted.

High wave number region (proof of Lemma 4.2)
Now we consider the case when either ξ or η is large.

Convergence of the ADI scheme (proof of Theorem 2.2)
Theorem 2.2 states that the error of the ADI method (2.4) has the same order as the implicit Milstein scheme (2.3). We now give a proof as follows.
and e n is the logarithmic error between the numerical solution and the exact solution introduced during [nk, (n + 1)k]. From (3.11), C n has the form In the low wave region, the numerical solutions are close to the exact solutions. We get from Taylor expansion that In the high wave region, we have where u = sin 2 hxξ 2 /( hxξ 2 ) 2 , v = sin 2 hyη 2 /( hyη 2 ) 2 . By the same reasoning as for the implicit scheme (2.3), the integration over the high wave region is of higher order than h 2 x and h 2 y given condition (2.9). Then the inverse Fourier transform gives the result. ✷

Numerical tests
In this section, we illustrate the stability and convergence results from the previous section by way of empirical tests.
where M (l) are independent Brownian motions and E L is the empirical mean with L = 100 samples.
Here, Figure 2(a) shows the convergence in h = h x = h y with k = 2 −12 small enough, which demonstrates second order convergence in h. Figure 2(b) shows the convergence in k with h = h x = h y = 2 −6 small enough to ensure sufficient accuracy of the spatial approximation. One can clearly observe first order convergence in k.
In Figure 3, we illustrate the dependence of the approximation error in the L 2 -norm on the correlation parameters. The error increases as a function of ρ x and ρ y . The error for ρ x = ρ y ≤ 0.3 (see Figure 3(a)) varies between roughly 10 −3 and 3 · 10 −3 , the error being smallest for ρ x = ρ y = 0 (the PDE case), and largest for large ρ x = ρ y and ρ xy between 0.1 and 0.4. For larger ρ x and ρ y , the error increases sharply. The stability region from Assumption 2.1 is marked in dark blue, which shows that stable results are obtained even outside the region where mean-square stability is proven. We found problems only for ρ x = ρ y ≥ 0.8. This discrepancy is partly due to the fact that Assumption 2.1 is sufficient, but not necessary, as some of the estimates are not sharp. Figure 3(b) shows a similar behaviour when varying ρ x and ρ y independently for fixed ρ xy .    Figure 4 shows the singular behaviour of the solution for large k and small h, as predicted by Theorem 2.1. Figure 5 investigates the behaviour of the error in this regime further, with h x = h y = h in Figure 5(a), and h y = 2 −1 fixed, h x = h in Figure 5(b). We calculate the L 2 error in space, and compare different scenarios. The top black line shows the error with k = 2 −2 fixed, and ρ x = ρ y = 0.6, ρ xy = 0.1. We can see that as h goes to zero, the error diverges with rate h −1/2 (choosing h y = 2 −1 fixed enables more refinements in h x to show the asymptotic behaviour better). The blue line second from top plots the error with k = 2 −4 fixed instead. Note that in this case the error will eventually diverge for h going zero, but this is not visible yet for this level of h. The next red line plots the error for ρ x = ρ y = 0, with k = 2 −2 fixed. Then the SPDE (1.6) becomes a PDE and divergence does not appear. Finally, the bottom blue dotted line plots the error for the SPDE (1.6) with initial condition v(0, x, y) For this smooth initial condition, the solution does not diverge for large k and small h. Hence, this verifies that the divergence is a result of the interplay of singular data and stochastic terms only, as shown in Corollary 2.2. We emphasise that the instability is so mild that it is only visible in artificial numerical tests, while for reasonably small k, in particular for k ∼ h 2 as would be chosen in practice, no instabilities occur.  6 An extended scheme and tests for a more general SPDE To approximate the general Zakai SPDE (1.1), by the Milstein scheme, we approximate in the last term The corresponding ADI implicit Milstein scheme is Here, {D i } 1≤i≤d are first order difference operators, and {D ij } 1≤i,j≤d are second order difference operators, X is the vector of mesh points ordered the same way as V , and, by slight abuse of notation, we denote by a(X), b(X), γ(X) the diagonal matrices such that each element of the diagonal corresponds to the function evaluated at the corresponding mesh point.
Notice the presence of an iterated Itô integral t+k t (M p (s) − M s (t)) dM l (s), called Lévy area, as is common in multi-dimensional Milstein schemes. It has been proved in [CC80,MG02] that there is no way to achieve a better order of strong convergence than for the Euler scheme by using solely the discrete increments of the driving Brownian motions. An efficient algorithm for the approximate simulation of the Lévy area has been proposed in [Wik01], building on earlier work in [KPW92,GL94] and based on an approximation of the distribution of the tail-sum in a truncated infinite series representation derived from the characteristic functions of these integrals. The best complexity of sampling a single path to obtain strong error ε is ε −3/2 , and the algorithm fairly complex.
Instead, to estimate this term in the time interval [t, t + k], we further divide the interval into O(k −1 ) steps and perform a simple Euler approximation. We find numerically that this still leads to first order convergence in time, and second order convergence in space. To balance the leading order error, the optimal choice is O(k) = O(h 2 ) = O(ε). Therefore, the estimate of the Lévy area in each time-step increases the computation time by O(k −1 ) = O(ε −1 ) for one step, whereas the matrix calculation for each time step is also O(ε −1 ), and hence the order of total complexity does not change.
Moreover, the path simulation including the Lévy areas can be performed separately beforehand using vectorisation, leading to further speed-up. Now we apply this method to an SPDE from [HK17], with Dirac initial u(0, x, y) = δ(x − x 0 )δ(y − y 0 ). This SPDE models the limit empirical measure of a large portfolio of defaultable assets in which the asset value processes are modelled by Heston-type stochastic volatility models with common and idiosyncratic factors in both the asset values and the variances, and default is triggered by hitting a lower boundary.
The notation for Y follows the same principle as above for X.
Figure 6(a) shows the density for a single Brownian path, with k = 2 −8 , h x = 5/16, and h y = 1/80. Figure 6(b) compares the computational cost under different time-stepping schemes: the Milstein scheme, a "modified" Milstein scheme, and the Euler scheme. Here, for the Milstein scheme we approximate the Lévy area by sub-timestepping as explained above, while in the "modified" Milstein scheme we drop t+k t (W s − W t ) dB s but keep the one-dimensional iterated integrals as they are known analytically. We expect that the latter will lead to a worse convergence in time (for non-zero ξ 1 , ρ 2,1 , ρ 1,1 ), which is verified in Figure 7(b).
In Figure 6(b), from a coarsest mesh with h x = 0.625, h y = 0.025, and k = 0.25, we keep decreasing the time-step k by a factor of 4, and the spatial mesh width by a factor of 2. This shows that the cost, measured by time elapsed in simulating one path, increases by a factor of 16, demonstrating that these three schemes result in the same order of complexity.
where V N 2i,2j (k, h x /2, h y /2) is the numerical solution to v(T, ih x , jh y ) with mesh size h x /2 and h y /2, and V N i,j (k, h x , h y ) uses a coarse mesh h x and h y . Both share the same Brownian path and same time step k, thus the univariate error in k cancels and we should see the correct convergence order in h x and h y . Here, Figure 7(a) shows the convergence in h = h x = h y with k = 2 −4 fixed, which demonstrates second order convergence in h.
Similarly, we study the error in k in terms of i,j h x h y E |V N i,j (k, h x , h y ) − V 2N i,j (k/2, h x , h y )| 2 1/2 , using now the difference between two solutions with same mesh size, same Brownian path, but different time-steps. Figure 7(b) shows the convergence in k with fixed h x = 5/8, h y = 1/40, under the Milstein scheme, "modified" Milstein scheme, and Euler scheme. The timestep k decreases by a factor of 4 from one level to the next. We also plot two blue dashed lines with slope 1/2 and 1 as reference. One can clearly observe first order convergence in k for the Milstein scheme, and half order convergence in k for the Euler scheme. As for the "modified" Milstein scheme, although it appears to converge with first order on coarse levels (due to dominance of the terms which converge with first order for this level of accuracy), the asymptotic order is seen to be lower.

Conclusions
We studied a two-dimensional parabolic SPDE arising from a filtering problem. We proved meansquare stability and pointwise as well as L 2 -convergence for a semi-implicit Milstein discretisation scheme. To reduce the complexity, we also implemented an ADI version of the scheme, and provided corresponding convergence results.
Further research is needed to analyse almost sure convergence, which is of interest for filtering applications and does not follow directly from our analysis.
Another open question is a complete analysis of the numerical approximation of initial-boundary value problems (as opposed to problems posed on R d ) for the considered SPDE, when the regularity at the boundary is lost. For example, for the 1-d SPDE with constant coefficients on the half-line, dv = −µ ∂v ∂x dt + 1 2 ∂ 2 v ∂x 2 dt − √ ρ ∂v ∂x dM t , (t, x) ∈ (0, T ) × R + , v(t, 0) = 0, with initial condition v(0, ·) ∈ H 1 , the second derivative can be unbounded, i.e., v(t, ·) / ∈ H 2 . This and more general forms have been studied in [KR81, BHH + 11]. In such cases, the assumptions on Galerkin approximations in papers previously mentioned such as in [BL12,BLS13] are not established in the literature, hence a new approach for the numerical analysis is to be developed.