A symmetry-adapted numerical scheme for SDEs

TL;DR

This paper introduces a geometric numerical scheme for stochastic differential equations that leverages Lie symmetries to preserve invariant properties, improving the accuracy and structure of numerical solutions.

Contribution

The paper presents a symmetry-adapted numerical scheme for SDEs, utilizing Lie symmetries to maintain invariant properties during approximation.

Findings

Applied to linear SDEs with successful symmetry preservation

Established theoretical error estimates for the proposed scheme

Demonstrated improved geometric properties in numerical solutions

Abstract

We propose a geometric numerical analysis of SDEs admitting Lie symmetries which allows us to individuate a symmetry adapted coordinates system where the given SDE has notable invariant properties. An approximation scheme preserving the symmetry properties of the equation is introduced. Our algorithmic procedure is applied to the family of general linear SDEs for which two theoretical estimates of the numerical forward error are established.

Equations

$$\nabla (F)=({\partial }_{x^i}(F^j)){|}_{\stackrel{j=1,...,k}{i=1,...,n}}.$$

$${\Phi }_{\ast }(Y)=(\nabla (\Phi )\cdot Y)\circ {\Phi }^{-1}.$$

$$X_t^i=X_0^i+{\int }_0^t{\mu }^i(X_s)ds+{\int }_0^t\sum _{\alpha =1}^m{\sigma }_{\alpha }^i(X_s)d{W}_s^{\alpha }+.$$

$${\int }_0^t\sum _{i,\alpha }({\sigma }_{\alpha }^i{)}^2(X_s)ds<+\infty ,{\int }_0^t\sum _i|{\mu }^i(X_s)|ds<+\infty$$

$$L=\sum _{\alpha =1}^m\sum _{i,j=1}^n\frac{1}{2}{\sigma }_{\alpha }^i{\sigma }_{\alpha }^j{\partial }_{x^i{x}^j}+\sum _{i=1}^n{\mu }^i{\partial }_{x^i}.$$

$$d{F}_t=L(f)(X_t)dt+\nabla (f)(X_t)\cdot \sigma (X_t)\cdot d{W}_t.$$

$$L(\Phi )$$

$$\nabla (\Phi )\cdot \sigma$$

$$Y(\mu )-L(Y)$$

$$\left[Y,{\sigma }_{\alpha }\right]$$

$$X_n^i=X_{n-1}^i+{\mu }^i(X_{n-1})\Delta t_n+\sum _{\alpha =1}^m{\sigma }_{\alpha }^i(X_{n-1})\Delta W_n^{\alpha },$$

$${\bar{X}}_n^i$$

$$\Delta {\mathbb{W}}_n^{1,1}=\frac{1}{2}((\Delta W_n{)}^2-\Delta t_n).$$

$${ϵ}_N={\left(\mathbb{E}[\Vert X_T-X_N{\Vert }^2]\right)}^{1/2}\le C(T,\mu ,\sigma )h^{1/2}.$$

$${\overline{ϵ}}_N={\left(\mathbb{E}[\Vert X_T-{\bar{X}}_N{\Vert }^2]\right)}^{1/2}\le \bar{C}(T,\mu ,\sigma )h.$$

$$x_n=F(x_{n-1},\Delta t_n).$$

$$\Phi (x_n)=F(\Phi (x_{n-1}),\Delta t_n).$$

$${\Phi }^{-1}(F(\Phi (x),\Delta t))=F(x,\Delta t)$$

$$Y^i(F(x,\Delta t))-Y^k{\partial }_{x^k}(F)(x,\Delta t)=0$$

$$x_n=F(x_{n-1},\Delta t,\Delta W^1,...,\Delta W^m).$$

$${\Phi }^{-1}(F(\Phi (x),\Delta t,\Delta W^{\alpha }))$$

$$Y^i(F(x,\Delta t,\Delta W^{\alpha }))-Y^k{\partial }_{x^k}(F)(x,\Delta t,\Delta W^{\alpha })$$

$$F^i(x)=x^i+{\mu }^i(x)\Delta t+{\sigma }_{\alpha }^i(x)\Delta W^{\alpha }.$$

$$0=Y_j(F^i)(x)-Y_j^i(F(x))$$

$$\begin{array}[]{c}Y^{i}_{j}(x)+\mu^{k}(x)\partial_{x^{k}}(Y^{i}_{j})(x)\Delta t+\frac{1}{2}\sum_{\alpha}\sigma^{k}_{\alpha}\sigma^{h}_{\alpha}\partial_{x^{k}x^{h}}(Y^{i}_{j})(x)\Delta t+\sigma^{k}_{\alpha}(x)\partial_{x^{k}}(Y^{i}_{j})(x)\Delta W^{\alpha}=\\ =Y^{i}_{j}(x+\mu\Delta t+\sigma_{\alpha}\Delta W^{\alpha}).\end{array}$$

$$Y_j^i(x)+{\mu }^k(x){\partial }_{x^k}(Y_j^i)(x)\Delta t+\frac{1}{2}\sum _{\alpha }{\sigma }_{\alpha }^k{\sigma }_{\alpha }^h{\partial }_{x^k{x}^h}(Y_j^i)(x)\Delta t+{\sigma }_{\alpha }^k(x){\partial }_{x^k}(Y_j^i)(x)\Delta W^{\alpha }=$$

$$=B_j^i+C_{j,k}^i{x}^k+C_{j,k}^i{\mu }^k(x)\Delta t+C_{j,k}^i{\sigma }_{\alpha }^k(x)\Delta W^{\alpha }$$

$$=B_j^i+C_{j,k}^i(x^k+{\mu }^k(x)\Delta t+{\sigma }_{\alpha }^k(x)\Delta W^{\alpha })$$

$$=Y_j^i(x+\mu \Delta t+{\sigma }_{\alpha }\Delta W^{\alpha }).$$

$$Y_j^i(x_0+{\sigma }_{\alpha }\Delta W^{\alpha })=Y_j^i(x_0)+{\sigma }_{\alpha }^k(x_0){\partial }_{x^k}(Y_j^i)(x_0)\Delta W^{\alpha }.$$

$$d{X}_t=\left(a\tanh (X_t)-\frac{b^2}{2}{\tanh }^3(X_t)\right)dt+b\tanh (X_t)d{W}_t,$$

$$Y=\tanh (x){\partial }_x$$

$${\Phi }_1(x)$$

Full text

1 Introduction

The exploitation of special geometric structures in numerical integration of both ordinary and partial differential equations (ODEs and PDEs) is nowadays a mature subject of the numerical analysis often called geometric numerical integration (see e.g. [17, 24, 28, 38]). The importance of this research topic is due to the fact that many differential equations in mathematical applications have some particular geometrical features such as for example a conservation law, a variational origin, an Hamiltonian or symplectic structure, a symmetry structure etc.. The development of geometrically adapted numerical algorithms permits to obtain suitable integration methods which both preserve the qualitative properties of the integrated equations and have a more efficient numerical behaviour with respect to the corresponding standard discretization schemes.

In comparison the study of geometric numerical integration of stochastic differential equations (SDEs) is not so well developed. In the current literature the principal aims consist in producing numerical stochastic integrators which are able to preserve the symplectic structure (see e.g. [1, 35, 40]), some conserved quantities (see e.g. [5, 23, 32]) or the variational structure (see e.g. [2, 3, 22, 42]) of the considered SDEs. For the study of the algebraic structure of stochastic expansions in order to achive an optimal efficient stochastic integrators at all orders see [11].

Although the exploitation of Lie symmetries of ODEs and PDEs (see e.g. [36]) to obtain better numerical integrators is an active research topic (see e.g. [4, 10, 31, 30] and references therein), the application of the same techniques in the stochastic setting to the best of our knowledge is not yet pursued, probably because the concept of symmetry of a SDE has been quite recently developed (see e.g. [6, 9, 8, 7, 16, 27, 29, 33]).

In this paper we introduce two different numerical methods taking advantage of the presence of Lie symmetries in a given SDE in order to provide a more efficient numerical integration of it.

We first propose the definition of an invariant numerical integrator for a symmetric SDE as a natural generalization of the corresponding concept for an ODE. When one tries to construct general invariant numerical methods in the stochastic framework, in fact, a non trivial problem arises. Since both the SDE solution as well as the Brownian motion driving it are continuous but not differentiable processes, the finite differences discretization could not converge to the SDE solution. We give some necessary and sufficient conditions in order that the two standard numerical methods for SDEs (the Euler and the Milstein schemes) are also invariant numerical methods. By using this result, in particular, we are able to identify a class of privileged coordinates systems where it is convenient to make the discretization procedure.

Our second numerical method, based on a well-defined change of the coordinates system, is inspired by the standard techniques of reduction and reconstruction of a SDE with a solvable Lie algebra of symmetries (see [8, 26]). Indeed a SDE with a solvable Lie algebra of symmetries can be reduced to a triangular system and, when the number of symmetries is sufficiently high, the latter can be explicitly integrated. In the stochastic setting the explicit integration concept is of course a quite different notion with respect to the deterministic one. Indeed the evaluation of an Ito integral, a necessary step in the reconstruction of a reduced SDE, can only be numerically implemented.

We apply our two proposed numerical techniques to the general linear SDEs, being the first non-trivial class of symmetric equations. In this case the two algorithmic methods can be harmonized in such a way to produce the same simple family of best coordinates systems for the discretization procedure. Interestingly, the identified coordinate changes are closely related to the explicit solution formula of linear SDEs. Although the integration formula of linear SDEs is widely known, it is certainly original the recognition of the proposed numerical scheme for linear SDEs as a particular implementation of a general procedure for SDEs with Lie symmetries. We finally point out that the simple case of affine drift and diffusion coefficients plays an important role since any SDE with real analytic drift and diffusion coefficients can be seen locally like it.

Moreover we theoretically investigate the numerical advantages of the new numerical scheme for linear SDEs. More precisely we obtain two estimates for the forward numerical error which, in presence of an equilibrium distribution, guarantee that the constructed method is numerically stable for any size of the time step $h$. This means that for any $h>0$ the error does not grow exponentially with the maximum-integration-time $T$, but it remains finite for $T\rightarrow+\infty$. This property is not shared by standard explicit or implicit Euler and Milstein methods. The obtained estimates can be considered original results mainly because the coordinate changes involved in the formulation of our numerical scheme are strongly not-Lipschitz, and so the standard convergence theorems can not be applied. Our theoretical results are also numerically illustrated.

It is interesting to note that the main part of the theory, in particular the definitions of strong symmetry of a SDE and of numerical scheme, can be easily extended to Stratonovich type SDE driven by general noises ([6]), for example by exploiting rough paths theory. Unfortunately, since the proofs of Theorem 5.1 and Theorem 5.2 use in essential way the (forward and bachward) Ito formula, the long terms estimates obtained here cannot be straightforwardly generalized to the rough paths driven SDEs framework. At the same time we think that some ideas in the proof of Theorem 5.2 can be suitable exploited to obtain pathwise estimates of long term error in the rough paths setting. The article is structured as follows: in Section 2 we recall the notion of strong symmetry of a SDE and we describe the two standard discretization schemes used in the rest of the paper. In Section 3 we propose two numerical procedures adapted with respect to the Lie symmetries of a SDE. We apply the proposed integration methods to general one and two-dimensional linear SDEs in Section 4 . In Section 5 some theoretical estimates showing the stability and efficiency of our adapted-to-symmetries numerical schemes in linear-SDEs are proved. In the last section we expose some numerical experiments confirming the theoretical estimates obtained in the previous section.

2 Preliminaries

2.1 Strong symmetries of SDEs

In the following $M$ will be an open subset of $\mathbb{R}^{n}$, and $x^{1},...,x^{n}$ will be the standard coordinates system of $M$. If $F:M\rightarrow\mathbb{R}^{k}$ we denote by $\nabla(F)$ the Jacobian of $F$ i.e. the matrix-valued function

[TABLE]

Furthermore we can identify the vector fields $Y\in TM$ with the functions $Y:M\rightarrow\mathbb{R}^{n}$, and if $\Phi_{M}\rightarrow M^{\prime}$ is a diffeomorphism (with $M^{\prime}\subset\mathbb{R}^{n}$) we introduce the pushforward

[TABLE]

Definition 2.1

Let $(\Omega,\mathcal{F},\mathcal{F}_{t},\mathbb{P})$ be a filtered probability space. Let $\mu$ and $\sigma$ be two smooth functions defined on $M$ and valued in n vectors and $n\times m$ matrices respectively. A solution to a SDE($\mu,\sigma$) is a pair $(X,W)$ of adapted processes such that

i) W is a $\mathcal{F}_{t}$-Brownian motion in $\mathbb{R}^{m}$;

ii) For $i=1,2,...,n$

[TABLE]

Remark 2.2

In particular all the integrals are meaningful if a.s.:

[TABLE]

Definition 2.3

A solution $(X,W)$ to a SDE($\mu,\sigma$) on $(\Omega,\mathcal{F},\mathcal{F}_{t},\mathbb{P})$ is said a strong solution if $X$ is adapted to the filtration $\mathcal{F}^{W}_{t}$ generated by the BM $W$ and completed with respect to $\mathbb{P}$.

Of course a solution $(X,W)$ is called a weak solution when it is not strong.

In this paper we fix a Brownian motion $W$, that is we consider only strong solutions of a SDE($\mu,\sigma$) and, consequently, we denote them simply by $X$. For a symmetry analysis via weak solutions of SDEs see [9].

A solution $X$ to a SDE($\mu,\sigma$) is a diffusion process admitting as infinitesimal generator:

[TABLE]

It is particularly useful for obtaining stochastic differentials the following celebrated formula.

Theorem 2.4** (Ito formula)**

Let $X$ be a solution of the SDE $(\mu,\sigma)$ and let $f:M\rightarrow\mathbb{R}$ be a smooth function. Then $F=f(X)$ has the following stochastic differential

[TABLE]

We recall important definitions of symmetries of a SDE.

Definition 2.5** (strong finite symmetry)**

We say that a diffeomorphism $\Phi$ is a (strong finite) symmetry of the SDE $(\mu,\sigma)$ if for any solution $X$ to the SDE $(\mu,\sigma)$ also $\Phi(X)$ is a solution to the SDE $(\mu,\sigma)$.

By using Ito formula it is immediate to prove the following result.

Theorem 2.6

A diffeomorphism $\Phi$ is a symmetry of the SDE $(\mu,\sigma)$ if and only if

[TABLE]

It is well-known that vector fields acting as infinitesimal generators of one parameter transformation groups are the most important tools in Lie group theory.

Definition 2.7** (strong infinitesimal symmetry)**

A vector field $Y$ is said a (strong infinitesimal) symmetry of the SDE $(\mu,\sigma)$ if the group of the local diffeomorphism $\Phi_{a}$ generated by $Y$ is a symmetry of the SDE $(\mu,\sigma)$ for any $a\in\mathbb{R}$.

The following determining equations for (any) infinitesimal symmetries are well-known (see, e.g., [16]). For their generalization to a weak solution case see [9].

Theorem 2.8** (Determining equations)**

A vector field $Y$ is an infinitesimal symmetry of the SDE $(\mu,\sigma)$ if and only if

[TABLE]

where $\sigma_{\alpha}$ is the $\alpha$-column of $\sigma$ ($\alpha=1,...,m$) and $[\cdot.\cdot]$ are the standard Lie brackets between vector fields.

2.2 Numerical integration of SDEs

For convenience of the reader, we recall the two main numerical methods for simulating a SDE and a theorem on the strong convergence of these methods (for a detailed description see e.g. [25]).

Consider the SDE having coefficients $(\mu,\sigma)$, driven by the Brownian motion $W$, and let $\{t_{n}\}_{n}$ be a partition of $[0,T]$. The Euler scheme for the equation $(\mu,\sigma)$ with respect to the given partition is provided by the following sequence of random variables $X_{n}\in M$

[TABLE]

where $\Delta t_{n}=t_{n}-t_{n-1}$ and $\Delta W^{\alpha}_{n}=W^{\alpha}_{t_{n}}-W^{\alpha}_{t_{n-1}}$. The Milstein scheme for the same equation $(\mu,\sigma)$ is instead constituted by the sequence of random variables $\bar{X}_{n}\in M$ such that

[TABLE]

where $\Delta\mathbb{W}^{\alpha,\beta}_{n}=\int_{t_{n-1}}^{t_{n}}{(W^{\beta}_{s}-W^{\beta}_{t_{n-1}})dW^{\alpha}_{s}}$. We recall that when $m=1$ we have that

[TABLE]

Theorem 2.9

Let us denote by $X_{t}$ the exact solution of a SDE $(\mu,\sigma)$ and by $X_{N}$ and $\bar{X}_{N}$ the N-step approximations according with Euler and Milstein scheme respectively. Suppose that the coefficients $(\mu,\sigma)$ are $C^{2}$ with bounded derivatives and put $t_{n}=\frac{nT}{N}$ and $h=\frac{T}{N}$. Then there exists a constant $C(T,\mu,\sigma)$ such that

[TABLE]

Furthermore when the coefficients $(\mu,\sigma)$ are $C^{3}$ with bounded derivatives then there exists a constant $\bar{C}(T,\mu,\sigma)$ such that

[TABLE]

Proof. See Theorem 10.2.2 and Theorem 10.3.5 in [25].

Theorem 2.9 states that $X_{N}$ and $\bar{X}_{N}$ strongly converge in $L^{2}(\Omega)$ to the exact solution $X_{T}$ of the SDE $(\mu,\sigma)$, where the order of the convergence with respect to the step size variation $h=\frac{T}{N}$ is $\frac{1}{2}$ in the Euler case and $1$ in the Milstein one.

Nevertheless the theorem gives no information on the behaviour of the numerical approximations when we fix the step size $h$ and we vary the final time $T$. In the standard proof of Theorem 2.9 one estimates the constants $C(T,\mu,\sigma)$ and $\bar{C}(T,\mu,\sigma)$ by proving that there exist two positive constants $K(\mu,\sigma),K^{\prime}(\mu,\sigma)$ such that $C(T,\mu,\sigma)=\exp(T\cdot K(\mu,\sigma))$ and $\bar{C}(T,\mu,\sigma)=\exp(T\cdot K^{\prime}(\mu,\sigma))$, by using Gronwall Lemma. In some situations the exponential growth of the error is a correct prediction (see for example [34]).

Of course this fact does not mean that in any case the errors $\epsilon_{n}$ and $\bar{\epsilon}^{\prime}_{n}$ exponentially diverge with the time $T$. Indeed if the SDE $(\mu,\sigma)$ admits an equilibrium distribution it could happen that the two errors remain bounded with respect to the time $T$. Unfortunately this favorable situation does not happen for any values of the step size $h$, but only for values within a certain region. The phenomenon just described is known as the stability problem for a discretization method of a SDE. This problem, and the corresponding definition, is usually stated and tested for some specific SDEs (see e.g. [19, 41] for the geometric Brownian motion, see e.g. [18, 39] the Ornstein-Uhlenbeck process, see e.g. [20, 21] for non-linear equations with a Dirac delta equilibrium distribution, and see e.g. [43] for more general situation). We will show some numerical examples of the stability phenomenon for general linear SDEs in Section 6.

3 Numerical integration via symmetries

3.1 Invariant numerical algorithms

When a system of ODEs admits Lie-point symmetries then invariant numerical algorithms can be constructed (see e.g. [31, 30, 10, 4]). By completeness we recall the definition of an invariant numerical scheme for a system of ODEs, in the simple case of one-step algorithms. The obvious extension for multi-step numerical schemes is immediate. The discretization of an ODEs system is a function $F:M\times\mathbb{R}\rightarrow M$ such that if $x_{n},x_{n-1}\in M$ are the $n,n-1$ steps respectively and $\Delta t_{n}$ is the step size of our discretization we have that

[TABLE]

If $\Phi:M\rightarrow M$ is a diffeomorphism we say that the discretization defined by the map $F$ is invariant with respect to the map $\Phi$ if it happens that

[TABLE]

If we require that the previous property holds for any $x_{n}\in\mathbb{R}^{n}$ and for any $\Delta t_{n}\in\mathbb{R}_{+}$ we get

[TABLE]

for any $x\in M$ and $\Delta t\in\mathbb{R}$. If $\Phi_{a}$ is an one-parameter group generated by the vector field $Y=Y^{i}(x)\partial_{x^{i}}$, by deriving the relation $\Phi_{-a}(F(\Phi_{a}(x),\Delta t))=F(x,\Delta t)$ with respect to $a$, we obtain the relation

[TABLE]

which guarantees that the discretization $F$ is invariant with respect to $\Phi_{a}$, generated by $Y$.

We can extend the previous definition to the case of a SDE in the following way. Let us discuss an integration scheme which depends only on the time $\Delta t$ and on the Brownian motion $\Delta W^{\alpha}_{n},\alpha=1,\dots,m$ (as for example the Euler method). The same discussion for integration methods depending also on $\Delta\mathbb{W}^{\alpha,\beta}_{n}$ or other random variables (as the Milstein method) is immediate. In the stochastic case the discretization is a map $F:M\times\mathbb{R}\times\mathbb{R}^{m}\rightarrow M$ and we have

[TABLE]

Equations (4) and (5) become

[TABLE]

Since Ito integral strongly depends on the fact that the approximation is backward (and not forward), we stress again that it is not easy to prove that a given discretization $X_{n}$ converges to the real solution of the SDE $(\mu,\sigma)$. For this reason we give a theorem which provides a sufficient (and necessary) condition in order that Euler and Milstein discretizations of a SDE are invariant with respect to any strong symmetries $Y_{1},...,Y_{r}$.

Theorem 3.1

Let $Y_{1},...,Y_{r}$ be strong symmetries of a SDE $(\mu,\sigma)$. When $Y^{i}_{j}=Y_{j}(x^{i})$ are polynomials of first degree in $x^{1},...,x^{n}$, then the Euler discretization (or the Milstein discretization) of the SDE $(\mu,\sigma)$ is invariant with respect to $Y_{1},...,Y_{r}$. If for a given $x_{0}\in M$, $\operatorname{span}\{\sigma_{1}(x_{0}),\dots,\sigma_{m}(x_{0})\}=\mathbb{R}^{n}$, also the converse holds.

Proof. We give the proof for the Euler discretization because for the Milstein discretization the proof is very similar. In the case of Euler discretization we have that

[TABLE]

The discretization is invariant if and only if

[TABLE]

Recalling that $Y_{j}$ is a symmetry for the SDE $(\mu,\sigma)$ and therefore it has to satisfy the determining equations (2) and (3), we have that the Euler discretization is invariant if and only if

[TABLE]

Suppose that $Y_{j}^{i}=B_{j}^{i}+C^{i}_{j,k}x^{k}$, then

[TABLE]

Conversely, suppose that the Euler discretization is invariant and so equality (8) holds. Let $x_{0}$ be as in the hypotheses of the theorem and choose $\Delta t=0$. Then

[TABLE]

Since $\Delta W^{\alpha}$ are arbitrary and $\operatorname{span}\{\sigma_{1}(x_{0}),...\sigma_{m}(x_{0})\}=\mathbb{R}^{n}$, $Y^{i}_{j}$ must be of first degree in $x^{1},...,x^{n}$.

Remark 3.2

The affinity of the coefficients $Y_{j}^{i}$ in Theorem 3.1 strongly depends on the fact that Euler and Milstein numerical approximations depend in an affine way from the noise $\Delta t,\Delta W^{\alpha},\Delta\mathbb{W}^{\alpha,\beta}$. If we coonsider a non affine numerical approximation we can have non affine symmetries $Y_{1},...,Y_{r}$ (see the discussion below).

Theorem 3.1 can be fruitfully applied in the following way. If $Y_{1},...,Y_{r}$ are strong symmetries of a SDE we search a diffeomorphism $\Phi:M\rightarrow M^{\prime}\subset\mathbb{R}^{n}$ (i.e. a coordinate change) such that $\Phi_{*}(Y_{1}),...,\Phi_{*}(Y_{r})$ have coefficients of first degree in the new coordinates system $x^{\prime 1},...,x^{\prime n}$. We discretize the transformed SDE $\Phi(\mu,\sigma)$ using the Euler discretization, obtaining a discretization $\tilde{F}(x^{\prime},\Delta t,\Delta W^{\alpha})$ which is invariant with respect to $\Phi_{*}(Y_{1}),...,\Phi_{*}(Y_{r})$. As a consequence the discretization $F=\Phi(\tilde{F}(\Phi^{-1}(x),\Delta t,\Delta W^{\alpha})$ is invariant with respect to $Y_{1},...,Y_{r}$. It is easy to prove that if the map $\Phi$ is Lipschitz we have that the constructed discretization converges in $L^{1}$ to the solution, while if the map $\Phi$ is only locally Lipschitz, the weaker convergence in probability can be established.

The existence of the diffeomorphism $\Phi$ allowing the application of Theorem 3.1 for general $Y_{1},...,Y_{r}$ is not guaranteed. Furthermore, even when the map $\Phi$ exists, unfortunately in general it is not unique. Consider for example the following one-dimensional SDE

[TABLE]

which has

[TABLE]

as a strong symmetry. There are many transformations $\Phi$ which are able to put $Y$ with coefficients of first degree, for example the following two transformations:

[TABLE]

Indeed we have that

[TABLE]

While the map $\Phi_{1}$ transforms equation (9) into a geometrical Brownian motion, the transformation $\Phi_{2}$ reduces equation (9) to a Brownian motion with drift. By applying Euler method by means of $\Phi_{1}$ we obtain a poor numerical result ( in fact $\Phi_{1}$ is not a Lipschitz function and in this circumstance errors are amplified). By exploiting $\Phi_{2}$ to make the discretization we obtain instead an exact simulation. The example shows that this first approach strongly depends on the choice of the diffeomorphism $\Phi$ (which has to be invertible in terms of elementary functions). So it is better to have another procedure able to individuate the best coordinate system for performing the SDE discretization.

3.2 Adapted coordinates and triangular systems

We introduce a further possible use of Lie’s symmetries in the numerical simulation of a SDE which turns out to be relevant only in the stochastic framework. Indeed in the deterministic setting one can obtain a completely explicit result.

Suppose that $M=M_{1}\times M_{2}$, with standard cartesian coordinates $x_{1}^{1},...,x_{1}^{r},x_{2}^{1},...,x^{n-r}_{2}$, and consider the following triangular SDE

[TABLE]

where $\mu^{i}_{1},\sigma^{i}_{1,\alpha}$ do not depend on $x^{i}_{1},...,x^{r}_{1}$. The above SDE is triangular in the variables $(x^{1}_{1},...,x^{r}_{1})$. By discretizing a triangular SDE $(\mu,\sigma)$ we reasonable aspect a better behavior than in the general case. Furthermore if $X^{1}_{2,t},...,X^{n-r}_{2,t}$ can be exactly simulated with $\sigma^{i}_{2,\alpha},\mu^{i}_{2}$ growing at most polynomially, we can conjecture that the error grows polynomially with respect to the maximal integration time $T$.

We recall that the triangular property of stochastic systems is closely related with their symmetries and in particular to SDEs with a solvable Lie algebra of symmetries. In order to briefly explain the connection between symmetries and the triangular form of SDEs, we introduce the following definitions (for more details see [8]).

Definition 3.3

A set of vector fields $Y_{1},...Y_{r}$ on $M$ is called regular on $M$ if, for any $x\in M$, the vectors $Y_{1}(x),...,Y_{r}(x)$ are linearly independent.

Definition 3.4

Let $Y_{1},...,Y_{r}$ be a set of regular vector fields on $M$ which are generators of a solvable Lie algebra $\mathcal{G}$. We say that $Y_{1},...,Y_{r}$ are in canonical form if there are $i_{1},...,i_{l}$ such that $i_{1}+...+i_{l}=r$ and

[TABLE]

where $G^{h}_{k}:M\rightarrow\operatorname{Mat}(i_{h},i_{k})$ are smooth functions.

Theorem 3.5

Let a SDE $(\mu,\sigma)$ admit $Y_{1},...,Y_{r}$ as strong symmetries and let us suppose that $Y_{1},...,Y_{r}$ constitute a solvable Lie algebra in canonical form. Then the SDE $(\mu,\sigma)$ assumes a triangular form with respect to $x^{1},....,x^{r}$.

Proof. The proof is an application of the determing equations and Definition 3.4 (see [8]).

As a notable consequence when we have a SDE $(\mu,\sigma)$ admitting a solvable regular Lie algebra $Y_{1},...,Y_{r}$ we can apply a methodology similar to the one proposed in the previous subsection. Indeed we can start by searching a map $\Phi:M\rightarrow M^{\prime}$ such that $\Phi(Y_{1}),...,\Phi(Y_{r})$ constitute a solvable Lie algebra in canonical form so implying that $\Phi(\mu,\sigma)$ is a triangular SDE. We can discretize $\Phi(\mu,\sigma)$ according with one of standard methods obtaining a discretization $\tilde{F}$. By composing $\tilde{F}$ with $\Phi$ we obtain a discretization $F(x,\Delta t,\Delta W^{\alpha})=\Phi^{-1}(\tilde{F}(\Phi(x),\Delta t,\Delta W^{\alpha})$ which, when $\Phi$ is Lipschitz, has the property of being a more simple triangular discretization scheme. Differently from Theorem 3.1, in the present situation we can always construct the diffeomorphism $\Phi$, as the following proposition states.

Proposition 3.6

Let $\mathcal{G}$ be an $r$-dimensional solvable Lie algebra on $M$ such that $\mathcal{G}$ has constant dimension $r$ as a distribution of $TM$. Then, for any $x_{0}\in M$, there exist a set of generators $Y_{1},...,Y_{r}$ of $\mathcal{G}$ and a local diffeomorphism $\Phi:U(x_{0})\rightarrow M^{\prime}$, such that $\Phi_{*}(Y_{1}),...,\Phi_{*}(Y_{r})$ are generators in canonical form for $\Phi_{*}(\mathcal{G})$.

Proof. See [8].

We conclude by pointing out that for a general solvable Lie algebra $Y_{1},...,Y_{r}$, the map $\Phi$, whose existence is guaranteed by Proposition 3.6, does not transform $\Phi_{*}(Y_{1}),...,\Phi_{*}(Y_{r})$ into a set of vector fields with coefficients of first degree in $x^{\prime 1},...,x^{\prime n}$. For this reason and by Theorem 3.1, the discretization $F$ constructed by using the diffeomorphism $\Phi$ and the usual Euler discretization algorithm is not invariant with respect to $Y_{1},...,Y_{r}$.

Nevertheless if we consider solvable Lie algebras satisfying a special relation, then $\Phi_{*}(Y_{1}),...,\Phi_{*}(Y_{r})$ will have coefficients of first degree in $x^{\prime 1},...,x^{\prime r}$.

Proposition 3.7

Suppose that the Lie algebra $\mathcal{G}=\operatorname{span}\{Y_{1},...,Y_{r}\}$ is such that $[[\mathcal{G},\mathcal{G}],[\mathcal{G},\mathcal{G}]]=0$. Then the coefficients of $\Phi_{*}(Y_{1}),...,\Phi_{*}(Y_{r})$ are of first degree in $x^{\prime 1},...,x^{\prime r}$. Moreover one can choose $\Phi$ such that the coefficients of $\Phi_{*}(Y_{1}),\dots,\Phi_{*}(Y_{r})$ are of first degree in all the variables $x^{\prime 1},...,x^{\prime n}$.

Proof. Let us suppose that $Y_{1},...,Y_{k}$ generates $\mathcal{G}^{(1)}=[\mathcal{G},\mathcal{G}]$. Then $\Phi^{*}(Y_{i})=(\delta^{l}_{i})$ for $i=1,...,k$. Using the fact that $[Y_{i},\mathcal{G}^{(1)}]\subset\mathcal{G}^{(1)}$ and the fact that $\Phi_{*}(Y_{1}),...,\Phi_{*}(Y_{r})$ are in canonical form, we must have that $\Phi_{*}(Y_{k+1}),...,\Phi_{*}(Y_{r})$ do not depend on $x^{\prime k+1},...,x^{\prime r}$ and their coefficients must be of first degree in $x^{\prime 1},...,x^{\prime r}$.

The second part of the proposition follows from the well known fact that when the vector fields $Z_{1},...,Z_{r}$ generate an integrable distribution, it is possible to choose a local coordinate system such that the coefficients of $Z_{1},...,Z_{r}$ do not depend on $x^{\prime r+1},...,x^{\prime n}$.

4 General linear SDEs

We first consider the one-dimensional linear SDE

[TABLE]

where $a,b,c,d\in\mathbb{R}$ and we apply the procedure previously presented in order to obtain a symmetry adapted discretization scheme.

Although it is possible to prove that equation (10) for $ad-bc\not=0$ does not admit strong symmetries (see [9]), we can look at equation (10) as a part of a two dimensional system admitting Lie symmetries.

Let us consider the system

[TABLE]

on $\mathbb{R}\times\mathbb{R}_{+}=M$, consisting of the original linear equation and the associated homogeneous one. It is simple to prove, by solving the determining equations (2) and (3), that the system (11) admits the following two strong symmetries:

[TABLE]

The more general adapted coordinate system system for the symmetries $Y_{1},Y_{2}$ is given by

[TABLE]

where $l\in\mathbb{R}$ and $f:\mathbb{R}_{+}\rightarrow\mathbb{R}$ is a smooth function. Indeed in the coordinate system $(x^{\prime},z^{\prime})^{T}=\Phi(x,z)$ we have that

[TABLE]

In order to guarantee that the Euler and Milstein discretization schemes are invariant, by Theorem 3.1 it is sufficient to choose $f(z)=-\frac{k}{z}$ for some constant $k$.

In the new coordinates the original two dimensional SDE becomes

[TABLE]

In the following, for simplicity, we consider the discretization scheme only for $l=0$. The Euler integration scheme becomes:

[TABLE]

and the Milstein scheme:

[TABLE]

We note that when $k=-\frac{d}{c}$ the two discretization schemes coincide.

Coming back to the original problem, in the Euler case we get:

[TABLE]

and in the Milstein case we obtain:

[TABLE]

Remark 4.1

There is a deep connection between equations (22) and (23) and the well-known integration formula for scalar linear SDEs. Indeed the equation (10) admits as solution

[TABLE]

where

[TABLE]

Equation (22) and (23) can be viewed as the equations obtained by expanding the integrals in formula (24) according with stochastic Taylor’s Theorem (see [25]). This fact should not surprise since the adapted coordinates obtained in Subsection 3.2 were introduced exactly to obtain formula (24) from equation (11). Since the discretizations schemes (22) and (23) are closely linked with the exact solution formula of linear SDEs we call them exact methods (or exact discretizations) for the numerical simulation of linear SDEs.

Let us now consider the following two dimensional SDE

[TABLE]

The previous equation can be solved explicitly. In particular the homogeneous linear part has solution given by (see, e.g. [12])

[TABLE]

where $\mu=\alpha+\frac{\sigma^{\prime 2}}{2}$. Thus the solution of the initial equation is

[TABLE]

The Euler discretization of the previous equation becomes:

[TABLE]

where $\Delta t_{n}=t_{n}-t_{n-1}$ and $\Delta W^{i}_{n}=W^{i}_{t_{n}}-W^{i}_{t_{n-1}}$.

5 Theoretical estimation of the numerical forward error for linear

SDEs

We provide an explicit estimation of the forward error associated with the exact numerical schemes proposed in the previous section for simulating a general linear SDE. The explicit solution of a linear SDE is well known and the use of the resolutive formula for its simulation is extensively used, but in the literature, to the best of our knowledge, there is no explicit estimation of the forward error.

5.1 Enunciates of the Theorems

Dividing $[0,T]$ in $N$ parts we obtain $N+1$ instants $t_{0}=0,t_{n}=nh,t_{N}=T$, with $h=\frac{T}{N}$. We denote by $X^{N,T}_{t}$ the approximate solution given by exact Euler method, $\bar{X}^{N,T}_{t}$ the approximate solution with respect to exact Milstein method and by $X_{t}$ the exact solution to the linear SDE. In the following we will omit $T$ where it is possible.

Theorem 5.1

For all $t,T\in\mathbb{R},t\in[0,T]$, we have

[TABLE]

where $h=\frac{T}{N}$, $g$ is a continuous function and $f$ is a strictly positive continuous function such that for $x\rightarrow+\infty$

[TABLE]

with $C(a,c)\in\mathbb{R}_{+}$.

Theorem 5.2

For all $t,T\in\mathbb{R},t\in[0,T]$, we have that

[TABLE]

where $h=\frac{T}{N}$, $\bar{g}$ is a continuous function and $f$ is a strictly positive continuous function such that for $x\rightarrow+\infty$

[TABLE]

with $C^{\prime}(a,c)\in\mathbb{R}_{+}$.

Before giving the proof of the two previous theorems we propose some remarks. We recall that a linear SDE with $ad-bc\not=0$ has an equilibrium distribution if and only if $a-\frac{c^{2}}{2}<0$. Furthermore the equilibrium distribution admits a finite first moment if and only if $a<0$ and a finite second moment if and only if $a+\frac{c^{2}}{2}<0$. Since we approximate the Ito integral up to the order $h^{1/2}$, the three cases in Theorem 5.1 follow from the fact that for giving an estimate of the error in Euler discretization we need a second moment control. More precisely we can expect a bounded error with respect to $T$ only when the second moment is finite as $T\rightarrow+\infty$.

Since in the Milstein case a finite first moment sufficies, in the second theorem we obtain that the error does not grow with $T$ when $a<0$. We can obtain an analogous estimate for the Euler method when $d=0$, i.e. in the case in which the Milstein and Euler discretizations coincide (situation similar to the additive-noise-SDEs setting). The use of only the first moment finitess for estimating the error has a price: indeed we obtain an $h^{1/2}$ dependence of the error. We remark that the techniques used in the proof of Theorem 5.2 exploit some ideas from the recent rough path integration theory (see e.g. [15]), and in particular this circumstance explains the $\frac{1}{2}$ order of convergence. This fact induces us to conjecture that our proof probably works also in the general rough path framework (for example for fractional Brownian motion by following [14]). If in Theorem 5.2 we do not require an uniform-in-time estimate, we can apply the methods used in the proof of Theorem 5.1 for obtaining an error convergence of order $1$.

Essentially the above theorems prove that for $a+\frac{c^{2}}{2}<0$ and for $a<0$ respectively, our symmetry adapted discretization methods are stable for any value of $h$. In Section 6 we give a comparison between the stability of the adapted-coordinates schemes with respect to the standard Euler and Milstein ones, via numerical simulations.

We conclude by noting that Theorem 5.1 and Theorem 5.2 cannot be deduced in a trivial way from the standard theorems about the convergence of Euler and Milstein methods (such as Theorem 2.9). Indeed the Euler and Milstein discretizations of equations (16) and (17) do not have Lipschitz coefficients. Furthermore even if a given discretization $(X^{\prime}_{n},Z^{\prime}_{n})$ of the system composed by (16) and (17) should converge to the exact solution in $L^{2}(\Omega)$, being the coordinate change $\Phi$ ( introduced in Section 4) not globally Lipschitz, it does not imply that the transformed discretization $(X_{n},Z_{n})$ converges to the exact solution $(X,Z)$ of the equation (11) in $L^{2}(\Omega)$. Finally, as pointed out in Subsection 2.2, Theorem 2.9 does not guarantee an uniform-in-time convergence as Theorem 5.1 and Theorem 5.2 instead state.

For proving the theorems we need the following two lemmas. The second allows to avoid very long calculations (see Appendix A).

Lemma 5.3

Let $W_{t}$ be a Brownian motion, $\alpha,\beta\in\mathbb{R}$ and $n\in\mathbb{N}$ then for any $t\in\mathbb{R}_{+}$

[TABLE]

is a continuous function of $t$ and in particular it is locally bounded. Moreover we have that

[TABLE]

Proof. The proof is based on the fact that $W_{t}$ is a normal random variable with zero mean and variance equal to $t$.

Lemma 5.4

Let $F:\mathbb{R}^{2}\rightarrow\mathbb{R}$ be a smooth function such that $F(0,0)=0$ and such that

[TABLE]

for some $\alpha\in 2\mathbb{N}$, for any $h$ and for some continuous function $L:\mathbb{R}\rightarrow\mathbb{R}_{+}$. Then there exists an increasing function $C:\mathbb{R}\rightarrow\mathbb{R}$ such that

[TABLE]

If furthermore $\partial_{w}(F)(0,0)=0$ and

[TABLE]

there exists an increasing function $C^{\prime}:\mathbb{R}\rightarrow\mathbb{R}$ such that

[TABLE]

Proof. The statements of the lemma derive as special cases from Lemma 5.6.4 and Lemma 5.6.5 in [25].

5.2 Proof of Theorem 5.1

We consider the case $t=T$. In fact we will find that our estimate is uniform for $t\leq T$. Using the notations in Remark 4.1 we can write $X_{T}=I_{1}+I_{2}$ where

[TABLE]

Also the approximation $X^{N}_{T}$ can be written as the sums of two integrals of the form $X^{N}_{T}=I^{N}_{1}+I^{N}_{2}$ where

[TABLE]

Obviously the strong error $\epsilon_{N}$ can be estimated by $\|I_{1}-I^{N}_{1}\|_{2}+\|I_{2}-I^{N}_{2}\|_{2}$, where hereafter $\|\cdot\|_{\alpha}=(\mathbb{E}[|\cdot|^{\alpha}])^{1/\alpha}$.

5.2.1 Estimate of $\|I_{1}-I^{N}_{1}\|_{2}$

Setting $\Psi_{s,t}=\Phi_{t}(\Phi_{s})^{-1}$ for any $s<t$, we obtain (with $\Delta t_{i}=h$)

[TABLE]

[TABLE]

[TABLE]

By Jensen’s inequality

[TABLE]

and by Fubini theorem we have to calculate $\mathbb{E}[(\Psi_{t,T}-\Psi_{t_{i-1},T})^{2}]$. Since

[TABLE]

and $\Psi_{s,t}=\Psi_{s,u}\Psi_{u,t}$ for any $s\leq u\leq t$ we obtain that

[TABLE]

because $\Psi_{t,T}$ and $\Psi_{t_{i-1},t}$ are independent as a consequence of the Brownian increments independence.

It is simple to note that the function

[TABLE]

satisfies $F_{1}(0,0)=0$ and, by Lemma 5.3,

[TABLE]

Thus, by Lemma 5.4, there exists an increasing function $C_{1}(h)$

[TABLE]

Using Lemma 5.3 we get

[TABLE]

obtaining

[TABLE]

where

[TABLE]

5.2.2 Estimate of $\|I_{2}-I^{N}_{2}\|_{2}$

We first consider $I_{2}=(d)\Phi_{T}\int_{0}^{T}{(\Phi_{t})^{-1}dW_{t}}$. Since Ito integral involves adapted processes we cannot bring $\Phi_{T}$ under the integral sign. However it is possible to take advantage of the backward integral formulation which allows to integrate processes that are measurable with respect to the (future) filtration $\mathcal{F}^{t}=\sigma\{W_{s}|s\in[t,T]\}$. In particular when $X_{s}$ is $\mathcal{F}^{t}$-measurable then

[TABLE]

where $\{t^{n}_{i}\}|_{i}$ is a sequence of $n$ points partitions of the interval $[0,T]$, having amplitude decreasing to [math] and the limit is understood in probability.

When $F$ is a regular function, $F(W_{t},t)$ is a process which is measurable with respect to both the filtrations $\mathcal{F}_{t}$ and $\mathcal{F}^{t}$; therefore one can calculate either $\int_{0}^{T}{F(W_{t},t)dW_{t}}$ and $\int_{0}^{T}{F(W_{t},t)d^{+}W_{t}}$.

The next well-known lemma says that we can write $I_{2}$ in terms of a backward integral, which allows to bring $\Phi_{T}$ under the integral sign.

Lemma 5.5

Let $F:\mathbb{R}^{2}\rightarrow\mathbb{R}$ be a $C^{2}$-function such that

[TABLE]

Then

[TABLE]

Proof. We report the proof for convenience of the reader (see, e.g., [37]). Setting

[TABLE]

since $F$ is $C^{2}$ then also $\tilde{F}$ is $C^{2}$. From this fact one deduces that

[TABLE]

By equating the two expressions one obtains the final formula.

Remark 5.6

Since in the proof of Lemma 5.5 the backward Ito integral plays a fundamental role, the same result cannot be obtained within a pathwise approach such as the rough paths one.

Since

[TABLE]

and $\partial_{w}(F)(w,t)=-cF(w,t)$, by Lemma 5.5, we can write

[TABLE]

Introducing $\tilde{I}_{2}=(d)\int_{0}^{T}{\Psi_{t,T}d^{+}W_{t}}$ and

[TABLE]

we have that

[TABLE]

We first consider the term $\|\tilde{I}_{2}-\tilde{I}^{N}_{2}\|_{2}$. The process $\tilde{I}^{N}_{2}$ can be written as $\int_{0}^{T}{(d)H_{t}dW^{+}_{t}}$ where $H_{t}$ is the $\mathcal{F}^{t}-$ measurable process given by

[TABLE]

where $1_{(t_{i-1},t_{i}]}$ is the characteristic function of the interval $(t_{i-1},t_{i}]$. By Ito’s isometry and Fubini’s Theorem we obtain

[TABLE]

Since Brownian motion has independent increments, we have that

[TABLE]

Introducing the function:

[TABLE]

which satisfies $H(0,0)=0$, by Lemma 5.4 and Lemma 5.3 we obtain

[TABLE]

where $C_{2}(h)$ is an increasing function and, finally,

[TABLE]

where

[TABLE]

In order to estimate the other term in the right-hand side of (36) we note that by introducing

[TABLE]

we have

[TABLE]

and

[TABLE]

By applying Lemma 5.5 to $K_{i}(t_{i},W_{t_{i}})$ we can write

[TABLE]

From the previous equality, by Ito isometry and Minkowski’s integral inequality we get

[TABLE]

where

[TABLE]

When $t_{i-1}<t\leq t_{i}$, by independence

[TABLE]

Introducing

[TABLE]

we have that $F_{2}(0,0)=F_{3}(0,0)=0$ and $\mathbb{E}[|\partial_{w}(F_{i})(t,W_{t_{i}}-W_{t})|^{2}]$, $\mathbb{E}[|\partial_{ww}(F_{i})(t,W_{t_{i}}-W_{t})|^{2}]$, $\mathbb{E}[|\partial_{t}(F_{i})(t,W_{t_{i}}-W_{t})|^{2}]\leq L(t_{i}-t)$ and so, by Lemma 5.4, there exist two continuous increasing functions $C_{3}(t),C_{4}(t)$ such that

[TABLE]

Since by independence

[TABLE]

analogously we can prove that there exists an increasing function $C_{5}$ such that

[TABLE]

For the second term in the right-hand side of (36), we have finally the following estimate

[TABLE]

where $G_{1}(T)$ and $G_{2}(T)$ are given by (35) and (39) respectively.

5.3 Proof of Theorem 5.2

We make the proof only for $a<0$, since in the other case the estimate are equal to the Euler case and can be addressed with the same proof. We introduce the two integrals

[TABLE]

5.3.1 Estimate of $\|I_{1}-\bar{I}^{N}_{1}\|_{1}$

First we note that (with $\delta t_{i}=h$)

[TABLE]

where we have taken, $n\in\mathbb{N}$, $\frac{1}{2n}+\frac{1}{\alpha}=1$ and $1<\alpha<2$ such that $\alpha a+\alpha(\alpha-1)\frac{c^{2}}{2}\leq 0$ (the last condition guarantees that when $T\rightarrow\infty$ we have $\mathbb{E}[\Psi_{t_{i},T}^{\alpha}]\rightarrow 0$). By Jensen’s inequality and Lemma 5.4 we can derive the following estimate:

[TABLE]

where $C_{5}(h)$ is an increasing function and in the last inequality we have used the fact that the function $F_{4}(t,W_{t})=\Psi_{t,h}-\Psi_{0,h}$ is such that $F_{4}(0,0)=0$. By Lemma 5.3, we have that

[TABLE]

and so

[TABLE]

where

[TABLE]

5.3.2 Estimate of $\|I_{2}-\bar{I}^{N}_{2}\|_{1}$

First we note that

[TABLE]

where $\alpha,n$ are as in the previous subsection. We introduce the following notation

[TABLE]

where we have used Lemma 5.5 and the fact that $\Psi_{s,t}=\Phi_{t}(\Phi_{s})^{-1}$. By introducing also $\hat{I}_{2,t_{i}}=\int_{t_{i-1}}^{t_{i}}{\Psi_{t,t_{i}}d^{+}W_{t}}$ and

[TABLE]

we have that

[TABLE]

It is simple to see that the two norms on the right-hand side of the previous expression do not depend on $t_{i}$ but only on the difference $h=t_{i}-t_{i-1}$, so we study the functions (with $\Psi_{t_{i},t_{i}}=1$):

[TABLE]

By a well-known consequence of Ito isometry (see, e.g., [13]) we can estimate the function $Z_{1}(h)$ as:

[TABLE]

where $D_{n}=(n(2n-1))^{n}$. Since the function

[TABLE]

satisfies $F_{5}(0,0)=\partial_{w}(F_{5})(0,0)=0$, by Lemma 5.4 there exists an increasing function $C_{6}(h)$ such that

[TABLE]

As far as concerned the function $Z_{2}(h)$, by introducing

[TABLE]

it is immediate to see that

[TABLE]

By applying Lemma 5.5 to $K(h,W_{h})$, and by noting that $K(h,W_{h})=0$, we obtain

[TABLE]

Since we have that $-\partial_{t}(K)(h,W_{h})+\partial_{ww}(K)(h,W_{h})/2+c\Psi_{0,h}=0$, and that $K(h,W_{h})=\partial_{w}(K)(h,W_{h})=\partial_{ww}(K)(h,W_{h})=0$, by Jensen’s inequality, Lemma 5.4 and by applying the same techniques used for obtaining (40) we find that

[TABLE]

or, equivalently,

[TABLE]

with the obvious definition of the function $C_{9}(h)$.

Finally we have

[TABLE]

where $G_{4}(T)$ is given by (41).

6 Numerical examples

We show some numerical experiments which confirm the theoretical estimate proved in Section 5 and permit to study other properties of the new discretization methods introduced in Section 4.

We simulate the linear SDE (10) with coefficients $a=-2$, $b=10$, $c=10$ e $d=10$. The coefficients are such that $a+\frac{c^{2}}{2}>0$ with $a<0$. This means that the considered linear equation admits an equilibrium probability density with finite first moment and infinite second moment. The coefficient $d$ has been chosen big enough to put in evidence the noise effect.

We make a comparison between the Euler and Milstein methods applied directly to equation (10) and the new exact methods (22) and (23) with the constants $k=0$ and $k=\frac{-d}{c}=-1$. In particular we observe that when $k=-1$, the schemes (22) and (23) coincide. We calculate the following two errors:

• •

the weak error $E^{w}=|\mathbb{E}[X_{t}-X^{N}_{t}]|$,

• •

the strong error $E^{s}=\mathbb{E}[|X_{t}-X^{N}_{t}|]$.

The weak error is estimated trought the explicit expression

[TABLE]

for the first moment of the linear SDE solution, and by using Monte-Carlo method with $1000000$ paths for calculating $\mathbb{E}[X^{N}_{t}]$. The strong error is estimated by exploiting Monte-Carlo simulation of $X_{t}$ and $X^{N}_{t}$ with $1000000$ paths. In order to simulate $X_{t}$ we apply the Milstein method with a steps-size of $h=0.0001$, for which we have verified that it gives a good approximation of both $\mathbb{E}[X_{t}]$ and the equilibrium density for $t\rightarrow+\infty$. Since we use Monte-Carlo methods for estimating $E^{w}$ and $E^{s}$, the two errors include both the systematic errors of the considered schemes and the statistical errors of the Monte-Carlo estimate procedure.

In Figure 1 we report the weak and strong errors with respect to the maximum time of integration $t$ which varies from $0.1$ to $1$ and stepsize $h=0.025$. As predicted by Theorem 5.2, the error of the exact method for $k=-1$ remains bounded. It is important to note that for the exact method in the case $k=0$ (where Theorem 5.1 and Theorem 5.2 do not apply) the errors remains bounded too, while for Euler and Milstein methods the errors grow exponentially with $t$.

In Figure 2 we report the weak and strong errors with respect to the maximum time of integration $t$, which varies from $0.1$ to $1$, and stepsize $h=0.01$. In this situation also the errors of the Mistein method remain bounded. In other words $h=0.01$ belongs to the stability region of the Milstein method but not to the stability region of the Euler method.

In Figure 3 we plot the weak and strong errors with fixed final time $t=0.5$ and steps number $N=10,...,80$, where the stepsize $h=\frac{t}{N}$. Here we note that the weak and strong errors for the exact methods do not change with the stepsize. This means that with a stepsize of only $h=0.05$ the exact methods have weak and strong systematic errors less than the statistical errors. Instead for the Milstein scheme the errors grow and only with a stepsize equal to $h=0.0125$ the systematic errors are comparable with the statistical ones. Equivalently we can say that the stability region is $[0,0.0125]$. In the Euler case the systematic error is not comparable with the statistical one.

In Figure 4 we report the total variation distance between the empirical probabilities of $X_{t}$ and of $X^{N}_{t}$ obtained simulating $1000000$ paths. We note that there is a big difference between the exact method for $k=0$ and for $k=-1$. The discrepancy is due to the fact that the exact method with $k=0$ tends to overestimate the points with probability less then $\frac{-d}{c}$ more than the Euler scheme does.

Now we simulate the two dimensional linear SDE analized in Section 4 by

choosing $\alpha=-20\mbox{, }\beta=-0.5\mbox{, }\sigma=\sigma^{\prime}=5\mbox{, }\textbf{c}=\textbf{e}=\begin{pmatrix}0.1&0.1\end{pmatrix}^{T}$ and $\textbf{d}=\begin{pmatrix}1&1\end{pmatrix}^{T}$. Our choice of the parameters guarantees the existence of an equilibrium probability density.

We compare approximated solutions obtained with the Euler method and with our exact method using $h=0.01$. To this end we calculated both the strong and weak componentwise error

[TABLE]

where $X^{i}_{t}$ is the $i-$th component of the solution. This time our true solution is calculated using the Euler method with timestep $h=0.0001$. As in the previous example the error are estimated using a Montecarlo simulation, this time with $10000$ paths, both for the approximated and the true solution. Again we expect $E^{w}_{i}$ and $E^{s}_{i}$ to include both systematic and statistical errors.

In Figure 5 and Figure 6 we compare the strong and weak error of both components of the simulated solutions with respect to the maximum time of integration varying from $0.1$ to $1$. As can be seen the error from our new method is bounded at all times while the Euler method errors show an exponential growth with respect to the maximum time.

In Figure 7 and Figure 8 we compare the errors of both approximations for solutions with $T=1$ and timestep size varying between $0.1$ to $0.01$. As in the previous one-dimensional case we can see how the new exact method gives a good approximation of the true solution even with large timesteps, while the Euler method fails to achieve the same magnitude of error even using a significative smaller timesteps.

7 Appendix

In the proof of Theorem 5.1, by using Lemma 5.3 and the independence of Brownian increments, we can estimate the errors in a very explicitely way. In particular without exploiting Lemma 5.4. We show main steps and final expressions.

From (33) we obtain that

[TABLE]

with

[TABLE]

Since $M_{1}(0)=\partial_{h}M_{1}(0)=0$, then $|M_{1}(h)|\leq M_{2}(h)h^{2}$ with $M_{2}(h):=max_{k\in[0,h]}|\partial_{h}^{2}M_{1}(k)|,$ and, finally,

[TABLE]

where $G_{1}(T)$ is given by (35), according with (34).

From (37) we obtain

[TABLE]

where

[TABLE]

Since $M_{3}(0)=\partial_{h}M_{3}(0)=0$, we have that $|M_{3}(h)|\leq M_{4}(h)h^{2}$ with $M_{4}(h):=max_{k\in[0,h]}|\partial_{h}^{2}M_{3}(k)|,$ and

[TABLE]

according with (38).

The second term on the right-hand side of (36) becomes

[TABLE]

where we have used independence and we have set

[TABLE]

We can obtain

[TABLE]

and, since $M_{5}(0)=\partial_{h}M_{5}(0)=0$, that $|M_{5}(h)|\leq M_{6}(h)h^{2}$, where $M_{6}(h):=max_{k\in[0,h]}|\partial_{h}^{2}M_{5}(k)|.$ Being:

[TABLE]

by putting $M_{9}(h)=M_{7}(h)M_{8}(h)$, one can easily verify that

[TABLE]

(because $M_{7}(0)=\partial_{h}M_{7}(0)=M_{8}(0)=\partial_{h}M_{8}(0)=0$) and, therefore, $|M_{9}(h)|\leq M_{10}(h)h^{4}$, where $M_{10}(h):=max_{k\in[0,h]}|\partial_{h}^{4}M_{9}(k)|.$ Finally

[TABLE]

that is

[TABLE]

with

[TABLE]

from which we get:

[TABLE]

to be compared with (40).