# Multilevel MC method for weak approximation of stochastic differential equation with the exact coupling scheme

• Yousef Alnafisah
From the journal Open Mathematics

## Abstract

Davie’s exact coupling technique for stochastic differential equations may be used to enhance the convergence of the multilevel Monte Carlo (MC) methodology. Giles developed the multilevel MC technique, which is based on executing the MC method several times with various time increments. It cuts computing costs significantly by executing most simulations at a low cost. The essential concept behind the multilevel MC approach with the exact coupling is discussed in this article. Numerical implementation reveals significant computational savings, which supports the analysis.

MSC 2010: 11-xx; 60-xx; 65-xx; 68-xx

## 1 Introduction

A stochastic differential equation is a mixture of deterministic and probabilistic elements. It is fundamental in many fields of application, and therefore, we find it with great interest in financial and actuarial mathematics, engineering applications, and insurance companies. In recent years, specialists in stochastic differential equations have been trying to find approximate solutions to it, either by using strong approximations or weak approximations. Strong approximation methods have had the most significant use in many scientific papers with a clear difference in applications and conditions. For example, when a stochastic differential equation is one-dimensional, there are often no difficulties in its application. When the Wiener process is multidimensional, the task gets more difficult. The previous articles [1,2,3] explore methods for estimating double integrals in any dimension using Fourier expansion. For the interested reader to learn more about the accomplishment of the simulation of the stochastic differential equation, we refer to [4,5]. Due to its superior properties for many problems, solutions to fractional stochastic differential equations (FSDEs) driven by the Brownian motion have recently received much attention from scientific researchers, see [6,7]. In [8], the authors look at a stochastic viscoelastic wave equation with nonlinear damping and logarithmic nonlinear source terms that established a blow-up result. On the other hand, a strong approximation of an Ito process is not always needed. It is possible that you are only interested in a function of the value of the Ito process at a certain point in time T , one of the first two moments, for example. This method and the use of expectation for some functions give an excellent approximation to the probability distribution of random variables. As a result, the kind of approximation that is needed in this case is significantly weaker than that provided by the pathwise convergence method. In this paper, we will apply the process of weak approximations using the technique of coupling developed by Davie [9].

This paper is structured as follows: Section 2 summarizes the recent stochastic differential equation (SDE) studies and addresses Davie’s methodology [9]. Section 3 discusses the integrated correct relation and the system of Euler. Section 4, offers a numerical example of convergence behavior. In this work, for the exact coupling, we use the multilevel Monte Carlo (MC) technique [10,11] on the two-dimensional SDEs.

## 2 Stochastic differential equations

Assume the SDE is as follows:

(2.1) d Y i ( t ) = O i ( t , Y ( t ) ) d t + k = 1 n B i k ( t , Y ( t ) ) d V k ( t ) , Y i ( 0 ) = Y i ( 0 ) ,

where i = 1 , , n on [ 0 , T ] , Y ( t ) is a n -dimensional vector, and V ( t ) is a n -dimensional path. Further, the coefficients B i k ( t , Y ( t ) ) satisfy a global Lipschitz condition O ( t , Y ) O ( t , y ) A Y y , and B ( t , Y ) B ( t , y ) A Y y , for all t [ t 0 , T ] and Y , y R , with A > 0 is a constant. Now assuming A and B are continuously on t for each Y , then there is a unique solution Y ( t ) to equation (2.1). To find an approximate solution on this interval [ 0 , T ] , we divided this interval into positive integer N equal length of intervals, i.e., S = T / N . Adding the following quadratic terms to the Euler scheme yields the Milstein scheme: k , l = 1 n χ i k l ( j S , y ( j ) ) P k l ( j ) . This will give the following scheme:

(2.2) y i ( j + 1 ) = y i ( j ) + A i ( j S , y ( j ) ) S + k = 1 n B i k ( j S , y ( j ) ) Δ V k ( j ) + k , l = 1 n χ i k l ( j S , y ( j ) ) P k l ( j ) ,

where Δ V k ( j ) = V k ( ( j + 1 ) S ) V k ( j S ) , P k l ( j ) = j S ( j + 1 ) S { V k ( t ) V k ( j S ) } d V l ( t ) , and χ i k l ( t , y ) = m = 1 q B m k ( t , y ) B i l y m ( t , y ) . If the commutativity condition

(2.3) χ i k l ( t , y ) = χ i l k ( t , y )

is satisfied for all y R d , t [ 0 , T ] , and all i , k , l , After that, the Milstein method is reduced to

(2.4) y i ( j + 1 ) = y i ( j ) + A i ( j S , y ( j ) ) S + k = 1 n B i k ( j S , y ( j ) ) Δ V k ( j ) + k , l = 1 n χ i k l ( j S , y ( j ) ) P k l ( j ) .

It should be noted that the previous strategy is solely dependent on the Brownian motion Δ V k ( j ) being implemented . We may use Δ V k ( j ) only and apply the special equations for the Milstein method. This can be implemented from: the observation that P k l ( j ) + P l k ( j ) = 2 F k l ( j ) , where F k l ( j ) = 1 2 Δ V k ( j ) Δ V l ( j ) for k l and F k k ( j ) = 1 2 { ( Δ V k ( j ) ) 2 S } . Also note that scheme (2.4) has the order 1 for n = 1 , but if n > 1 , we obtain the order 1 2 . As delineated in the previous paper [9], we tend to modify to scheme (2.4), which can provide the order 1 with the invertible diffusion. As a result, we will have the exact coupling strategy as follows.

### 2.1 Exact coupling scheme

For the sake of brevity and simplicity, we consider scheme (2.7) with the explicit version and B i k ( y ) depend only on y , and moreover, set the drift term as zero, so

(2.5) y i ( j + 1 ) = y i ( j ) + k = 1 n B i k ( y ( j ) ) Y k ( j ) + k , l = 1 n χ i k l ( y ( j ) ) ( Y k ( j ) Y l ( j ) S δ k l ) .

In [12], we show the order of the exact coupling method, and we obtained the following scheme:

(2.6) Y i ( r , j ) = Y i ( r + 1 , 2 j ) + Y i ( r + 1 , 2 j + 1 ) + k , l = 1 n τ i k l ( Y k ( r + 1 , 2 j + 1 ) Y l ( r + 1 , 2 j ) Y l ( r + 1 , 2 j + 1 ) Y k ( r + 1 , 2 j ) ) + O ( ( S ( r ) ) 3 / 2 ) ,

where τ i k l = 1 2 j c i j χ i k l . and S ( r ) = T 2 r is the step size. Under a condition of nondegeneracy, the interpretation of the normal distribution that was generate in scheme (2.4) will be changed, which in turn will converge to the first order, see [9]. So we apply (2.4) with the increment Y k ( j ) replacing Δ V k ( j ) , obtaining the following approximate solution y i ( j ) ,

(2.7) y i ( j + 1 ) = y i ( j ) + k = 1 n b i k ( y ( j ) ) Y k ( j ) + k , l = 1 n ρ i k l ( y ( j ) ) ( Y k ( j ) Y l ( j ) S δ k l ) ,

when the increment Y k ( j ) are independent normal distribution, i.e., N ( 0 , S ) . Notice that (2.7) is the same as scheme (2.4) except Δ V k ( j ) is substituted by Y k ( j ) and we do not assume Δ V k ( j ) as equal to Y k ( j ) . We now need

Z i k = 1 n b i k ( j S , x ( j ) ) X k ( j ) + k , l = 1 n ρ i k l ( j S , x ( j ) ) ( X k ( j ) X l ( j ) S δ k l ) .

In [12], we demonstrate that the local error of the exact coupling scheme is E y ( r , 1 ) y ( r + 1 , 2 ) 2 D 2 a 2 S 3 , where a and D 2 are functions of y ( r , j ) .

### 2.2 Weak order of convergence

We thus set the step size S , which splits the interval [ 0 , T ] into equal steps. Y ( t ) is the solution for the SDE, where C is a positive constant that is independent of S . If any smooth function f is used

E ( f ( y S ) ) E ( f ( Y ( T ) ) ) C S γ , S ( 0 , 1 ) ,

then, at time T = N S , a discrete time approximations y S , with the step size S , converges to the solution Y ( t ) with order γ .

## 3 Multilevel method for weak order of SDEs

Giles is developing the multilevel MC technique and based on using the MC technique many times with various time increments. So, let us assume that we have the SDE.

(3.1) d x i ( t ) = k = 1 n b i k ( x ( t ) ) d W k ( t ) , x i ( 0 ) = x i ( 0 )

where i = 1 , , q on [ 0 , T ] , x ( t ) is a n -dimensional vector, and W ( t ) is a n -dimensional Brownian path. We will estimate E f ( x ( T ) ) , by the mean of N successive estimates, i.e., 1 N i = 1 N f ( x S ( i ) ) , where x S ( i ) is a solution of (3.1) and f is globally Lipschitz. The multilevel approach has the benefit of reducing the estimation’s computational load. In the multilevel MC simulations, we examine MC simulations with various time steps S ( r ) = T 2 r , while the standard method depends on a fixed time step. Now, we assume that μ r = E f ( x S ( r ) ) . Then, we will estimate μ n , where n is large enough.

From Giles [10], we have

(3.2) μ n = μ 0 + μ 1 μ 0 + μ 2 μ 1 ± + μ n 1 μ n 2 + μ n μ n 1 = μ 0 + k = 1 n ( μ k μ k 1 ) .

The multilevel MC approach depends on approximating all terms in the right-hand side of (3.2) independently.

First, for μ 0 , we have

(3.3) μ 0 = E f ( x S 0 )

and

(3.4) Y 0 = 1 N 0 i = 1 N f ( x S 0 ( i ) ) .

Next, for μ k μ k 1 , we have

(3.5) μ k μ k 1 = 1 N k i = 1 N k ( f ( x S ( k ) ( i ) ) f ( x S ( k 1 ) ( i ) ) ) Y ˆ k .

Thus, the right-hand side of (3.2) becomes

(3.6) 1 N 0 i = 1 N f ( x S 0 ( i ) ) + k = 1 n 1 N k i = 1 N k ( f ( x S ( k ) ( i ) ) f ( x S ( k 1 ) ( i ) ) ) .

It is obvious that the quantity ( f ( x S ( k ) ) f ( x S ( k 1 ) ) ) is the outcome of two approximations with successive time steps S ( k ) and S ( k 1 ) , but the same Brownian path. The variance may be reduced by using the same Brownian path, so that a smaller value for N k can be obtained.

Now, V k represents the variance of a single sample ( f ( x S ( k ) ) f ( x S ( k 1 ) ) ) . Y ˆ k simple estimator’s variance is expressed as follows:

(3.7) Var ( Y ˆ k ) = Var 1 N k i = 1 N k ( f ( x S ( k ) ( i ) ) f ( x S ( k 1 ) ( i ) ) ) = V k N k .

As a result, if the variance of the starting level is

Var ( Y 0 ) = 1 N 0 i = 1 N Var ( f ( x S 0 ( i ) ) ) = V 0 N 0 ,

then the variance of Y ˆ = k = 1 n Y ˆ k is expressed as follows:

(3.8) Var ( Y ˆ ) = Var k = 0 n Y ˆ k = k = 0 n Var ( Y ˆ k ) = k = 0 n V k N k .

Furthermore, the computational load being proportional to

k = 0 n N k S k .

Therefore, by choosing N k to be proportional to V k S k , the variance, for a fixed computational load, is minimized.

The following theorem presents the general application of the multilevel MC approach, see [10]. By applying this theorem to the exact coupling, we obtain the computational load results.

## Theorem 3.1

Let f ( x ( T ) ) be a functional of the solution of SDE (3.1) with a given Brownian path W ( t ) , and f ( x S ( k ) ) represents the corresponding approximation to x ( T ) , utilizing a time step numerical discretization S k = T M k .

If there exist independent estimators Y ˆ k using N k MC simulation, and the positive constants α 1 2 , β , b 1 , b 2 , b 3 such that

1. E [ f ( x S ( k ) ) f ( x ( T ) ) ] b 1 S k α

2. E [ Y ˆ k ] = E [ f ( x S 0 ) ] , k = 0 E [ f ( x S ( k ) ) f ( x S ( k ) 1 ) ] , k > 0

3. V [ Y ˆ k ] b 2 S k β N k

4. B k , the computational complexity of Y ˆ k , is bounded by

C k b 3 N k S k .

Then there is a b 4 positive constant., s.t. for any ε < e 1 , there are certain values n and N k for which the multi-level estimator:

Y ˆ = k = 0 n Y ˆ k

has a limit on the mean-square-error (MSE)

MSE = E [ ( Y ˆ E [ f ( x ( T ) ) ] ) 2 ] < ε 2

having a computational cost B with bound

B b 4 ε 2 , β > 1 , b 4 ε 2 ( log ε ) 2 , β = 1 , b 4 ε 2 ( 1 β ) / α , 0 < β < 1 .

## Proof

Look at Giles [10].□

The previous theorem may be used to obtain β , which is twice the scheme’s strong order.

We’ll examine the order of the variance in the following applications, i.e., V k = O ( S k β ) , and use Theorem 3.1 to calculate the computational load.

## 4 The multilevel MC method for the exact coupling method using scheme (2.7)

In [13], Giles and Szpruch have worked on the multilevel approach and obtained O ( ε 2 ) computational load although they employ a more stringent demand for the function’s regularity f than Lipschitz. So, we showed in [14] that under nondegeneracy, the order of method (2.7) with the exact coupling is O ( S ) , and a single sample’s variance is

(4.1) V k = V ( f ( x S ( k ) ) f ( x S ( k 1 ) ) ) = O ( S k 2 ) .

As a result of (4.1) and because f is a Lipschitz function, the vector and matrix functions fulfill the Lipschitz conditions with uniformly bounded derivatives, and we have the weak order O ( S ) for (2.7) with the invertibility condition for matrix b i k ( t , X ( t ) ) . So that

E ( f ( x S ( k ) ) f ( x ( T ) ) = O ( S k ) .

Thus, the single sample variance V k is O ( S k 2 ) for the standard estimator (3.5). In addition, the best option for N k is proportional to ( S k 3 / 2 ) asymptotically. As a result, since the variance is β = 2 > 1 , the computational load will be O ( ε 2 ) , according to Theorem 3.1.

As a consequence, the number of simulations N k will drop rapidly as the variance of the multilevel MC for scheme (2.7) with the exact coupling decreases. So the main order of the convergence of this method will be obtained from the first level of the simulation, which will be the dominant term with the computational load O ( ε 2 ) .

### 4.1 Numerical implementation

In this section, we will look at invertible two-dimensional SDE equations and demonstrate numerically how well the multilevel MC method performs for the exact coupling method (2.7). We will also look at how reducing the number of simulations affects the computing load. To see the outcomes, the following two-dimensional invertible SDEs will be used for the Matlab implementation.

(4.2) d X 1 ( t ) = ( sin ( X 2 ( t ) ) ) 2 d W 1 ( t ) 1 1 + X 1 2 ( t ) d W 2 ( t ) , d X 2 ( t ) = 1 1 + X 2 4 ( t ) d W 1 ( t ) + ( cos ( X 1 ( t ) ) ) 2 d W 2 ( t ) , for 0 t 1 , with X 1 ( 0 ) = 2 and X 2 ( 0 ) = 0 .

W 1 ( t ) and W 2 ( t ) are two independent standard Brownian motions.

At each level of simulation, the Lipschitz function f ( x S ( k ) ) = sin ( x S ( k ) ) will be used to estimate the following amount ( f ( x S ( k ) ) f ( x S ( k 1 ) ) ) , where x S ( k ) and x S ( k 1 ) are two discrete approximations with distinct time steps S ( k ) and S ( k 1 ) . However, they both follow the Brownian path. The number of simulations N k , the step-size S , and the result of the simple estimator Y ˆ k with its confidence interval are presented in Table 1 for each level.

Table 1

The multilevel MC with the exact coupling two-dimensional SDE

Level k ( N k ) The number of simulations The step size S = T 2 ( n 1 ) ( Y ˆ k ) The simple estimator Confidence interval of Y ˆ k
1 4,194,304 1 0.0973 (0.0977, 0.0969)
2 1,482,910 0.5 0.0248 (0.0242, 0.249)
3 524,288 0.25 0.0074 (0.0072, 0.0077)
4 185,364 0.125 0.0033 (0.0031, 0.0037)
5 65,536 0.0625 0.0015 (0.0013, 0.0017)
6 23,170 0.03125 0.00075 (0.00058, 0.00090)
7 8,192 0.015625 0.00052 (0.00040, 0.00063)
8 2,896 0.0078125 0.00015 (0.000040, 0.00027)
9 1,024 0.00390625 0.000033 ( 0.00006 , 0.00013)

Now we will implement a Matlab code to obtain the basic estimator’s findings at each level, i.e., Y ˆ k in (3.5) and the number of levels will be M = 9 . For the estimate, the starting value of the number of simulations N k will be 2 22 , and also the same Brownian path will be used for each simulation. We calculate the 90% confidence interval of Y ˆ k , ( μ ψ , μ + ψ ) , given its mean value μ in the numerical simulations, where

(4.3) ψ = Z α / 2 × σ N k .

The sample size is N k , the confidence coefficient is Z α / 2 , and the standard deviation is σ . We compute the standard deviation as follows, and we obtain a N k estimate for the amount ( f ( x h ( k ) ) f ( x h ( k 1 ) ) ) , for each level M . We might then express their values as a vector and compute their standard deviation σ using the Matlab function (Std).

We obtain two terms for the multilevel MC. The first comes from the estimate of the top level μ 0 , which includes a step-size S 0 , i.e.,

E f ( x S ) = 1 N i = 1 N f ( x S ( i ) ) 0.6840 .

The remaining terms are derived from the sum of the simple estimator Y ˆ k , which will be

k = 1 M 1 N k i = 1 N k ( f ( x S ( k ) ( i ) ) f ( x S ( k 1 ) ( i ) ) ) 0.1364 .

As a result, the estimate for the exact coupling multilevel MC technique is expressed as follows:

(4.4) μ M = 1 N i = 1 N f ( x S ( i ) ) + k = 1 M 1 N k i = 1 N k ( f ( x S ( k ) ( i ) ) f ( x S ( k 1 ) ( i ) ) ) 0.6840 + 0.1364 = 0.8204 .

## 5 Conclusion

We established a weak convergence method to numerically solve stochastic differential equations with a possibly degenerate diffusion coefficient. We used an exact coupling in a two-dimensional case, which has a condition that the SDE should be invertible.

1. Conflict of interest: The author states no conflict of interest.

## References

[1] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag, Berlin, 1995. Search in Google Scholar

[2] T. Rydén and M. Wiktrosson, On the simulation of iteraled Itô integrals, Stochastic Process Appl. 91 (2001), 151–168. 10.1016/S0304-4149(00)00053-3Search in Google Scholar

[3] M. Wiktorsson, Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent Brownian motions, Ann. Appl. Probab. 11 (2001), 470–487. 10.1214/aoap/1015345301Search in Google Scholar

[4] Y. Alhojilan, Explicit order 3/2 Runge-Kutta method for numerical solutions of stochastic differential equations by using Itô-Taylor expansion, Open Math. 17 (2019), 1515–1525. 10.1515/math-2019-0124Search in Google Scholar

[5] Y. Alnafisah, The exact coupling with trivial coupling (Combined Method) in two-dimensional SDE with non-invertiblity matrix, Dyn. Syst. Appl. 28 (2019), 111–142. Search in Google Scholar

[6] Y. Alnafisah and H. M. Ahmed, Neutral delay Hilfer fractional integrodifferential equations with fractional brownian motion, Evol. Equ. Control Theory 11 (2022), no. 3, 925–937, http://doi.org/10.3934/eect.2021031. Search in Google Scholar

[7] H. M. Ahmed, Noninstantaneous impulsive conformable fractional stochastic delay integro-differential system with Rosenblatt process and control function, Qual. Theory Dyn. Syst. 21 (2022), 15, http://doi.org/10.1007/s12346-021-00544-z. Search in Google Scholar

[8] A. Benramdane, N. Mezouar, M. Alqawba, S. Boulaaras, and B. Cherif, Blow-up for a stochastic viscoelastic lamé equation with logarithmic nonlinearity, J. Funct. Spaces 2021 (2021), 9943969, http://doi.org/10.1155/2021/9943969. Search in Google Scholar

[9] A. M. Davie, Pathwise approximation of stochastic differential equations using coupling, www.maths.ed.ac.uk/sandy/coum.pdf. Search in Google Scholar

[10] M. B. Giles, Multilevel Monte Carlo path simulation, Oper. Res. 56 (2008), 607–617. 10.1287/opre.1070.0496Search in Google Scholar

[11] M. B. Giles, Improved multilevel Monte Carlo convergence using the Milstein scheme, in: Monte Carlo and Quasi-Monte Carlo Methods, A. Keller, S. Heinrich, H. Niederreiter, eds Springer, Berlin, Heidelberg, 2008, pp. 343–358, http://doi.org/10.1007/978-3-540-74496-2_20. Search in Google Scholar

[12] Y. Alnafisah, Comparison between Milstein and exact coupling methods using MATLAB for a particular two-dimensional stochastic differential equation, J. Inf. Sci. Eng. 36 (2020), 1223–1232. Search in Google Scholar

[13] M. B. Giles and L. Szpruch, Antithetic multilevel Monte Carlo estimation for multidimensional SDEs without Lévy area simulation, arXiv:1202.6283. Search in Google Scholar

[14] Y. Alnafisah, Order-one convergence for exact coupling using derivative coefficients in the implementation, Dyn. Syst. Appl. 28 (2019), 573–585. Search in Google Scholar

[15] N. Fournier, Simulation and approximation of Lévy-driven stochastic differential equations, ESIAM Probab. Stat. 15 (2011), 233–248. 10.1051/ps/2009017Search in Google Scholar

[16] I. Gyöngy and N. Krylov, Existence of strong solutions for Itôas stochastic equations via approximations, Probab. Theory Related Fields 105 (1996), 143–158. 10.1007/BF01203833Search in Google Scholar

[17] P. E. Kloeden, E. Platen, and I. Wright, The approximation of multiple stochastic integrals, Stoch. Anal. Appl. 10 (1992), 431–441. 10.1080/07362999208809281Search in Google Scholar

[18] A. Davie, KMT theory applied to approximations of SDE, in: Stochastic Analysis and Applications, D. Crisan, B. Hambly, T. Zariphopoulou (eds), Springer, 2014, pp. 185–201. 10.1007/978-3-319-11292-3_7Search in Google Scholar

[19] Y. Alnafisah, Exact coupling method for Stratonovich stochastic differential equation using non-Degeneracy for the diffusion, IEEE Access 7 (2019), 7442–7447. 10.1109/ACCESS.2018.2888945Search in Google Scholar