The implementation comparison between the Euler and trivial coupling schemes for achieving strong convergence

1.   Introduction

In recent decades, there has been an active research movement looking at numerical simulations of solutions to stochastic differential equations (SDEs). The numerical solution of SDEs plays a prominent role in many fields in mathematics, physics, engineering, finance, and other applied sciences. Around a century ago, Bachelier ^[1] employed Brownian motion or the Wiener process to simulate the values of stocks in the Paris Bourse, which is now a well-known technique. Later, Einstein ^[2] employed a comparable mathematical concept in his research on Brownian motion. After that, Wiener ^[3] expanded on the mathematical theory of Brownian motion. Ito ^[4] made a further advancement by laying the groundwork for the stochastic calculus now known as the Ito calculus.

SDEs are commonly used to model financial quantities such as asset prices, interest rates, and derivatives. The solutions to SDEs are continuous-time stochastic processes, unlike ordinary differential equations (ODEs), which have a unique solution for each initial condition if certain conditions on it and its derivatives are met. To accommodate stochastic dynamics, methods for computing SDE solutions are adaptations of those used for solving ODEs.

For SDEs, analytical solutions are frequently unavailable, hence, I must use numerical approaches to solve the problem approximately. The numerical solutions of SDEs were addressed (see, for example, ^[5,6,7]). In ^[8], Farnoosh et al. suggested some analytical and numerical solutions for SDEs. In ^[9], Zhan considered the numerical approximation of random periodic solutions for SDEs. On the other hand, Yin and Gan ^[10] added an error-correcting term to a Milstein scheme to solve for solid SDEs.

Recently, there have been significant advancements in the numerical analysis of SDEs. One issue with Kloeden and Platen's method, which uses a stochastic Taylor series expansion, is that the double stochastic integrals become difficult to express in terms of smaller stochastic integrals when the Wiener process is multidimensional. However, the Fourier series expansion of the Wiener process in multidimensional cases has been used to represent the double integrals ^[11,12,13]. Despite this, many random variable generations are required for each cycle.

To achieve robust convergence for SDEs, Davie ^[14] used coupling and provides order one. Kanagawa ^[15] examined the convergence rate between approximations of solutions with i.i.d random variables using two probability metrics. Fournier ^[16] used the quadratic Vaserstein distance to approximate the Euler scheme, and Rio ^[17] provided a highly precise rate of convergence for the central limit theorem in Vaserstein distance ^[14].

Rio ^[18,19] also produced accurate bound estimates, while Cruzeiro et al. designed a modified Milstein scheme that yields an order one for strong approximation and achieves an order one approach under the non-degeneracy condition. Charbonneau et al. ^[20] analyzed the Vaserstein bound ^[21] using weak convergence and the Strassen-Dudley theorem. Gyongy and Krylov ^[22] used coupling to prove that an approximation to a strong solution on a particular probability space converges. Davie ^[23] applied the Vaserstein constraint to vector SDE solutions and uses the Komlos, Major, and Tusnady theorem to obtain an approximation of order one under the non-degeneracy assumption. It is recommended that these papers be studied thoroughly as they have important implications for understanding this topic.

1.1.   Preliminary definitions

Numerous significant continuous-time Markov processes, such as the Ornstein-Uhlenbeck and Bessel processes can be characterized as solutions to SDEs with drift and diffusion coefficients that are solely dependent on the process's current value. For a one-dimensional process with a one-dimensional Brownian motion, the general form of such an equation is

where $\left(\mathcal{H}_{u}\right) _{u\, \geq \, 0}$ is a standard Wiener process.

Definition 1.1. Assume that $\left(\mathcal{H}_{u}\right) _{u\, \geq \, 0}$ is a standard Brownian motion on a space of probability $(\Omega, \digamma, \mathcal{P})$ with a filtration $\mathbb{F} = (\digamma _{u})_{u\geq 0}$. An adapted process and continuous process $\mathcal{S}_{u} = \mathcal{S}_{u}^{v}$ with the continuous paths satisfying

is called a strong solution of SDE (1.1) with the initial condition $\mathcal{S}_{u_{0}}\in \mathbb{R}$.

For the integrals in Eq (1.2) to be well defined there must be some degree of regularity in $\mathcal{S}_{u}$ and the functions $a$ and $b$. In particular, the following conditions should be achieved

and

Furthermore, for all $u < \infty$, the solution must be exist with probability one.

One of the most important properties of the stochastic integral is that

For more details on stochastic integral, see ^[11].

The vector SDE is defined from the scalar ones, just like the relationship between scalar and vector stochastic differentials. We suppose $\left(\mathcal{H}_{s}\right) _{s\, \geq \, 0}$ is an $m$–dimensional Wiener process with scalar components $\mathcal{H}_{s}^{r},$ $r = 1, 2, ..., m$ are independent w.r.t a common family of $\sigma _{i}$–algebras $(\digamma _{s})_{s\geq 0}$ where $i = 1, 2, ...$. Further, we define a $d$–dimensional vector function $\mathfrak{A}$ and a $d\times m$–matrix function $\mathfrak{B}$ as

and

Then, the $d$–dimensional vector SDE is

The $d-$stochastic integral equation is

1.2.   Existence and uniqueness

The fundamental finding, credited to Ito, is that SDE (1.1) has strong solutions for uniformly Lipschitz functions $a\left(u, \mu \right)$ and $b\left(u, \mu \right)$, and has a unique solution for all initial values $\mathcal{S}_{u_{0}} = \kappa$. In what follows, the initial instant $u_{0}\in \left[0, \mathcal{U}_{0}\right]$ is fixed and arbitrary, and the coefficient functions $a,$ $b:$ $\left[0, \mathcal{U}_{0}\right] \times \mathbb{R}\rightarrow \mathbb{R}$ are given. For the theorem of existence, we need the following conditions

${\rm{E1}}:$ (Measurability): The coefficients functions $a = a\left(u, \mu \right)$ and $b = b\left(u, \mu \right)$ are jointly $\mathcal{L}^{2}$ –measurable in $\left[0, \mathcal{U}_{0}\right] \times \mathbb{R}$.

${\rm{E2}}:$ (Lipschitz condition): There is a $\ell > 0$ (Lipschitz constant) such that

and

for every $u\in \left[0, \mathcal{U}_{0}\right]$ and $\mu, \eta \in \mathbb{ R}$.

${\rm{E3}}:$ (Linear growth bound): There is $k$ (a constant) such that

and

for every $u\in \left[0, \mathcal{U}_{0}\right]$ and $\mu \in \mathbb{R}$.

${\rm{E4}}:$ (Initial value): $\mathcal{S}_{u_{0}}$ is $\digamma _{u_{0}}$ –measurable with $\mathrm{E}\left(\left\vert \mathcal{S} _{u_{0}}\right\vert ^{2}\right) < \infty$.

Theorem 1.1. [11, Theorem 4.5.3] Assume that hypotheses E1–E4 hold, then, SDE (1.1) has a pathwise unique strong solution $\mathcal{S}_{u}$ for $u\in \left[0, \mathcal{U}_{0} \right]$ with

1.3.   Strong order of convergence

By assuming that $\mathrm{s}_{h}$ is a discrete-time approximation with a step size $\hslash$ that converges strongly of order $\delta$ at time $\top = M\hslash$ to the solution $\mathcal{S}_{u}$ of an SDE, we mean that the expected value of the absolute difference between $\mathrm{s}_{\hslash }$ and $\mathcal{S}_{U}$ is less than or equal to $C\hslash ^{\delta }$, where $\hslash$ is in the interval $(0, 1)$. Here, the interval $\left[0, \top \right]$ is evenly divided into the $M = \frac{\top }{ \hslash }$ sub-interval with length $\hslash$. Also, we suppose $C > 0$ is a constant, which is independent of $\hslash$.

2.   Numerical techniques for solving SDEs

Numerous numerical techniques exist for solving SDEs. In this context, I will discuss two significant methods. The first is the Euler-Maruyama approach, which provides a strong convergence rate of $\frac{1}{2}$. The second one is the trivial coupling method, which also yields a strong convergence rate of $\frac{1}{2}$. I will illustrate their convergence behavior through numerical examples. Consider the following SDE:

where $i = 1, 2, ..., q.$

Assuming I have an $m-$dimensional vector $\mathcal{S}(u)$ and a $d$ –dimensional Brownian path $\mathcal{H}(u)$ defined on an interval $[0, \top]$, I can approximate the solution by dividing the interval $[0, \top]$ into $M$ equal sub-intervals of the length $\hslash = \top /M$.

2.1.   Euler-Maruyama method

Euler's method is a simple numerical method for approximating solutions to SDEs. The method is based on discretizing the time interval and approximating the SDE by a difference equation. Here is how the method works:

where $\Delta \mathcal{H}_{k}^{(j)} = \mathcal{H}_{k}((j+1)\hslash)- \mathcal{H}_{k}(j\hslash)$.

Where $\mathcal{H}_{k}$ is a Wiener process increment (i.e., a Gaussian random variable with mean zero and variance $\hslash$) and f and g are functions that determine the drift and diffusion coefficients of the SDE.

The approximation $\mathrm{s}_{i}$ obtained by this method is called the Euler-Maruyama approximation. Note that this approximation is a random variable since it depends on the random values of $\mathcal{H}_{1}, ..., \mathcal{H}_{M}$. Euler's method is a $\frac{1}{2}$-order method, which means that the approximation error is proportional to $\hslash$. Many other numerical SDE methods for SDEs, such as the Milstein or the Runge-Kutta method, have higher-order accuracy. However, Euler's method is simple, easy to implement, and often used as a starting point for more sophisticated methods.

2.2.   Derivation of trivial coupling from the Milstein method

The Milstein scheme is a numerical method for approximating the solutions to SDEs. It is an extension of the Euler-Maruyama that scheme can provide a higher order of convergence.

The Milstein scheme is obtained by adding the quadratic terms.

to the Euler scheme, which will give the following scheme

where

and

The Milstein scheme is a first-order method, meaning the approximation error is proportional to $\hslash ^{2}$. This is an improvement over the $\frac{1}{2}$-order Euler-Maruyama scheme. The Milstein scheme is also more accurate when the diffusion coefficient is highly nonlinear. I will not use the Milstein method to find approximate solutions to SDEs; I want to use it to find the trivial coupling method. I can obtain the trivial coupling from Milstein's process by using the fact mentioned later in this paper from the Brownian motion symmetry.

Both methods are easy to apply to the MATLAB program, with some properties of the Brownian motion in the trivial coupling scheme.

2.3.   Two-level approximation

To approximate the solution to SDE (1.3) using Euler or other numerical schemes, I must generate the increments $\Delta \mathcal{H}_{k}^{(j)}$. To achieve this, I can use Levy's construction of the Brownian motion to simulate a sequence of approximations that converge to the solution.

That is,

where $r\in \mathbb{N}$ and

with $\hslash ^{(r)} = \frac{\top }{2^{r}}.$

The two-level approximation in Eq (2.4) is referred to as the trivial coupling. To generate the normal distribution in Eq (2.4) for a given level $r$, I can first generate the increments in the left hand side $\Delta \mathcal{H}_{k}^{(r, j)}$, and then conditionally generate the increments in the right hand side. I repeat this process for each level $r+2$, $r+3$, and so on until i obtain the Brownian path $\mathcal{H}(u)$.

It is difficult to extend the Milstein scheme to $d\geq 2$, but I can implement a special class of equations using only the $\Delta \mathcal{H} _{k}^{(j)}$. This can be achieved by observing that $A_{kl}^{(j)}+A_{lk}^{(j)} = 2B_{kl}^{(j)}$, where

if $k\neq l$ and

3.   The implementation of Euler and trivial coupling schemes

This section will focus on examining two-dimensional SDEs and testing the convergence of the Euler and trivial coupling schemes. The equations being considered are given by:

There are two separate standard Brownian motions, $G_{1}(t)$ and $G_{2}(t)$. The accuracy of the final solutions obtained by the Euler and trivial coupling methods will be compared. Generally, obtaining explicit solutions to SDEs is impossible, so simulation methods are used to approximate the solutions. To estimate the error of an approximation, I can compare the results obtained with two different step sizes. Specifically, I can take the difference between the results obtained using a step size of $\hslash$ and those obtained using a step size of $\hslash /2$.

For the SDE given by Eq (3.1), the MATLAB code is used to estimate the absolute error $\epsilon = \frac{1}{R}E(\left\vert \mathrm{s}_{(r)}-\widehat{\mathrm{s} }_{(r)}\right\vert),$ where $\mathrm{s}_{(r)}$ is the Euler and trivial coupling approximation to (3.1) using a step size of $\hslash$, $\widehat{\mathrm{s}}_{(r)}$ is another approximation for (2.1) using a step size of $\hslash /2$ and $r$ is a range of values of $\hslash$. The same Brownian path will be used for both approximations. The MATLAB code will be run with different steps (10, 20, 40, 80,160 and 320) over many paths (e.g., $R = 10,000$).

The Euler and trivial coupling schemes have achieved strong convergence with a value of $\delta = \frac{1}{2}$. A previous study (^[24,25,26]) demonstrated that the proposed algorithm converges faster than existing algorithms. Since analytical solutions for SDEs are typically unknown, I use approximate solutions to compare the two methods. To determine the approximation error for each method, I calculate the absolute error for varying numbers of steps, using the same number of simulations for both methods. I generate 10,000 different Brownian paths with varying step sizes over the interval $[0, 1]$, and the results for the experimental error and elapsed time for the trivial method are shown in Table 1. I use THE MATLAB code with different step sizes and many paths $R$ to obtain the results in Table 1.

Figure 1 depicts a logarithmic plot of the absolute error for six different time steps, showing that the trivial scheme achieves strong convergence with an order of one-half. Both methods use the same six-step sizes of $0.1, 0.05, 0.025, 0.0125, 0.00625,$ and $0.003125.$ Table 1 presents the trivial scheme's experimental error and elapsed time. Table 1 and Figure 1 indicate that the strong approximation error decreases as the step size decreases. The expected strong convergence for the trivial scheme is one-half, as stated in Alnafisah's paper (^[24,25,26]).

On the other hand, I have the Euler method, a simple and widely used numerical method for solving SDEs. It involves approximating the solution of the SDE at each step by using the derivative of the solution evaluated in the previous step. This approximation is then used to compute the solution at the current step. The Euler method is computationally inexpensive and straightforward to implement, but it has some limitations. For example, it is only accurate for small time steps and may introduce significant errors when the SDE has nonlinearities or the solution changes rapidly. Nonetheless, the Euler method remains useful and popular for solving SDEs, especially as a starting point for more sophisticated numerical methods. The Euler method is an alternative to the previous trivial coupling method.

So, command [Error(1, i) = log(Eulererrorq2('b', [2; 0], 1, S(1, i)))] calculates the absolute value of the difference between the approximate solution $\mathrm{s}_{(r)}$ and the solution $\widehat{\mathrm{s}}_{(r)}$ of the SDE with different step sizes. Table 2 provides the experimental error for each of the six-time steps and the elapsed time for the Euler method.

Figure 2 is a plot of $\log \epsilon$ against $\log(h)$. It is easy to see that the slope is $0.0.49037$, which again shows the strong convergence of $\frac{1}{2}$-order method.

Comparing the results in Tables 1 and 2, I observe in both methods that as step size decreases, the estimate of the absolute error also decreases. I can also observe in the previous tables and plots that the trivial coupling and Euler methods strongly converge with order one-half. I emphasize that I applied these methods over the same Brownian paths ($R$ = 10,000) for the same step sizes. It can also be seen that using the Euler scheme can reduce the total computational time. I see from the tables that there is a significant difference between the elapsed time. The trivial method code takes much time to obtain the result, especially when the Brownian path is very big, but the Euler method code takes a few seconds.

4.   Conclusions

Generally, an explicit solution for an SDE cannot be obtained. Therefore, I used simulation to approximate the solution and determine the convergence behavior. In this study, I used MATLAB to simulate the trivial coupling and Euler methods and obtained an approximate solution for the SDE. Both methods exhibited an order of one-half convergence. I then applied these methods to a two-dimensional SDE to compare their efficiency. Furthermore, I calculated the error values for the trivial coupling and Euler methods to compare their strong order and computation time. Our findings suggest that the Euler method was faster than the trivial method for solving invertible SDEs. However, the trivial method required assuming the nondegeneracy condition for the diffusion term and was computationally expensive. In contrast, the Euler method does not require this condition. Therefore, I can conclude that the Euler method is more effective than the trivial method for solving invertible SDEs.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The researcher would like to thank the Deanship of Scientific Research, Qassim University, for funding the publication of this project.

Conflict of interest

This work does not have any conflict of interest.