A numerical scheme for stochastic differential equations with distributional drift

by   Tiziano De Angelis, et al.
University of Leeds

In this paper we present a scheme for the numerical solution of stochastic differential equations (SDEs) with distributional drift. The approximating process, obtained by the scheme, converges in law to the (virtual) solution of the SDE in a general multi-dimensional setting. When we restrict our attention to the case of a one-dimensional SDE we also obtain a rate of convergence in a suitable L^1-norm. Moreover, we implement our method in the one-dimensional case, when the drift is obtained as the distributional derivative of a sample path of a fractional Brownian motion. To the best of our knowledge this is the first paper to study (and implement) numerical solutions of SDEs whose drift cannot be expressed as a function of the state.



There are no comments yet.


page 1

page 2

page 3

page 4


Taming singular stochastic differential equations: A numerical method

We consider a generic and explicit tamed Euler–Maruyama scheme for multi...

A New Discretization Scheme for One Dimensional Stochastic Differential Equations Using Time Change Method

We propose a new numerical method for one dimensional stochastic differe...

Semi-implicit Taylor schemes for stiff rough differential equations

We study a class of semi-implicit Taylor-type numerical methods that are...

Multi-dimensional Avikainen's estimates

Avikainen proved the estimate E[|f(X)-f(X)|^q] ≤ C(p,q) E[|X-X|^p]^1/p+1...

A generalized Avikainen's estimate and its applications

Avikainen provided a sharp upper bound of the difference E[|g(X)-g(X)|^q...

Randomized derivative-free Milstein algorithm for efficient approximation of solutions of SDEs under noisy information

We deal with pointwise approximation of solutions of scalar stochastic d...

Fréchet derivatives of expected functionals of solutions to stochastic differential equations

In the analysis of stochastic dynamical systems described by stochastic ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

The aim of our paper is to obtain a numerical algorithm capable of approximating the solution of a multi-dimensional SDE of the form


where is a distribution (in a subspace of Schwartz distributions which will be specified later) and is a -dimensional Brownian motion. Existence and uniqueness of a solution to (1) was first derived in [3], where the authors give a mathematical meaning to the term by introducing a new concept of solution, so-called virtual solution. The latter is indeed needed since is a distribution and cannot be evaluated in . Further work on equations of a similar kind can be found, for example, in [2, 7, 8].

The literature on numerical methods for SDEs with low regularity in the drift coefficient is vast and we will review some of the most recent results later in this introduction. However, to the best of our knowledge our work is the first one to address the question for a class of SDEs whose drift is so irregular that it cannot be expressed as a function of the state . This improvement on the existing results hinges on the concept of virtual solution given by [3] which links the SDE in (1

) to a class of partial differential equations (PDEs) with distributional drift studied in

[6]. It is worth emphasising that our algorithm does not require a numerical solution of the PDE and instead it deals directly with the SDE in (1). Hence, a simpler version of the methods that we use here can be adopted to complement/extend the existing studies on numerical schemes for SDEs whose drift is a function with low regularity.

In this paper we devise a 2-step numerical algorithm for the solution of (1). In the first step we approximate the distributional drift with a sequence of functions , , that converges in a suitable sense to . In the second step we adopt a Euler-Maruyama (EM) approximation of the SDEs associated to the ‘more regular’ drift (we denote by the solutions of such SDEs). As a direct consequence of this approach, relying on results contained in [3] and [4], we obtain convergence in law of the approximating process to the solution of (1) (see Proposition 3.1).

In order to obtain a rate of convergence we restrict our attention to the one-dimensional case. In particular, one of the main results of the paper (Theorem 3.3) is to show that the -rate of convergence of to is controlled by the rate of convergence of to with respect to a suitable Sobolev norm. There are two main ingredients in our proof: (i) the concept of virtual solution and (ii) a bound on the local time of continuous semi-martingales (see Lemma 5.1 and Proposition 5.4). The latter is a technical result of independent interest based on an application of Itô-Tanaka formula.

In the one dimensional setting we also obtain a (weak) -rate of convergence of the EM scheme (applied to ) for SDEs with bounded and piece-wise Lipschitz drift (Theorem 3.5). In particular, given a mesh with time points we find a convergence rate of

for some and arbitrary , where is the EM approximation of . This result is in line with recent work [16], where a (strong) -rate of is obtained for any , and [18] where analogous results are obtained in a multi-dimensional setting with respect to a (strong) -norm (here ‘strong’ refers to the fact that [16] and [18] consider the -norm on the space of continuous paths, under expectation). While we were completing our work another paper has appeared ([17]) in which a (weak) -rate of up to can be found for EM schemes on 1-dimensional diffusions with possibly discontinuous drift (with Sobolev-Slobodeckij type regularity).

Our method is different from the ones already employed in the literature because we rely on the link between virtual solutions and PDEs (hence complementing existing results). We do not investigate in this paper (strong) -rates of convergence for the EM scheme because the convergence in the first step of our algorithm (that is, convergence of to ) is in and we prefer to maintain this symmetry for clarity of exposition. However, a deeper study of the ideas that we develop in Section 6 seems in order and we leave it for future work.

One important aspect of our paper is that, due to the distributional nature of our drift, the actual numerical implementation of the scheme is non-trivial and, in particular, the choice of the approximating functions , in the first step of the algorithm, needs to be addressed carefully. It turns out that a convenient choice is to use Haar wavelets to construct the sequence , . The main reasons for this choice are: (i) these wavelets form a basis for the Sobolev spaces of negative order which are needed to accommodate the original drift ; (ii) they enjoy the so-called multi-resolution property, which improves the computational efficiency of the algorithm (see further details in Section 7); (iii) since Haar wavelets are piecewise constant functions, their values can be stored exactly in a computer, hence adding no error to the computations. To illustrate the fine intricacies of our numerical scheme we discuss in detail, in Section 7, the actual implementation of the algorithm when the drift is obtained as the distributional derivative of a sample path of a fractional Brownian motion. To the best of our knowledge, this is the first time a numerical scheme is implemented for SDEs with distributional drift.

The literature on EM approximation of SDEs with irregular drift is very vast and here we only provide a short overview. In the few paragraphs above we have discussed the contributions given in the papers [16], [17], [18] and [4] (the latter in particular derives convergence in law of the scheme and uses it as a tool for a proof of strong existence of solutions to a class of SDEs). Related results can also be found in [13] where a -rate of convergence of is obtained for (possibly degenerate) multi-dimensional SDEs. Notice that in [13] the EM scheme is applied to a process obtained as a suitable transformation of the solution of the SDE. Similar ideas were also used in [14] to find a (strong) -rate of convergence of ; differently from [13], the convergence in [14] is for the approximation of the original SDE. In the case of non-degenerate SDEs with irregular coefficients (strong) rates of convergence can be found in [19] (multi-dimensional setting, rate ) and [20] (one-dimensional setting, rate ).

Going back a few years we find other contributions to the study of EM schemes for multi-dimensional SDEs with discontinuous drift. For example [5] proves strong convergence in but with no rate and [29] proves convergence in law and obtains a (weak) -rate of convergence under the assumption of Hölder coefficients in a one-dimensional setting. A scheme in two steps is analysed in [10], where authors first regularise the drift of their SDE and then apply EM scheme to the more regular process (our approach follows the same spirit). Of course there are also numerous results in the case of SDEs with smooth coefficients. A detailed review is difficult and outside the scopes of our paper (for a general introduction one may refer to [9]). However, we would like to mention that weak convergence with rate up to is obtained in [1], and [15] contains further results in that direction.

Our paper is organised as follows. In Section 2 we introduce the necessary notation and give a rigorous meaning to the SDE (1). In Section 3 we present the main results of the paper, whose proofs are then provided in Sections 5 and 6. Background material on SDEs with distributional drift, which is needed to understand our arguments of proof, is presented in Section 4. Section 7 is devoted to implementing our numerical scheme in the case when the drift is obtained as the distributional derivative of one path of a fractional Brownian motion. The paper is completed by two technical appendices that account for important properties of Haar wavelets and their use in the simulation of a fractional Brownian motion and of its sample-path’s derivative (in distributional sense).

2. Notations and setting

In this section we introduce the theoretical framework in which equation (1) is well defined and we recall useful results from [3] on its solution. Throughout the paper we will use and for the spatial gradient and Laplacian of a function, respectively, and for its partial derivative with respect to time.

For all and , we denote by the fractional Sobolev spaces (or Bessel-potential spaces) defined as the images of through fractional powers of , that is . These spaces are Banach spaces when equipped with the norm

For more details see [25].

We observe that if then does actually contain distributions, while when it only contains functions. For we have the special case . These spaces have the inclusion property , for , and we will use the notation . This implies that if then is an element of all spaces for . We will also consider

-dimensional vector fields

such that each component . These elements belong to the space which we denote simply by for ease of notation. Analogously we denote their norm by .

For any Banach space we denote by the space of -valued continuous functions of time. This is again a Banach space when endowed with the norm . For future reference we also introduce on the equivalent norm


We denote by the closure of Schwartz functions with respect to the norm . For simplicity of notation we just write . These spaces will be used to track the regularity in the space variable . Further, we will also use the space of differentiable functions with -Hölder first derivatives for , that is the space

where the norm is defined as

where is the Euclidean norm in (when no confusion shall arise we will simply write ). Again for simplicity of notation we write in place of . Finally, we denote

We will work in the framework of [3] and make the following assumption.

Assumption 1.

Let for , , and .

It was shown in [3, Theorem 28], under assumptions more general than our Assumption 1, that for every there exists a unique in law virtual solution to (1). The definition of virtual solution was originally introduced in [3, Definition 25] and it is given in terms of a stochastic basis and a -adapted, continuous stochastic process (shortened as ) such that the integral equation


holds for all

, with probability one. Here

is the mild solution in of the following parabolic Kolmogorov-type PDE


with . In [3, Theorem 14] the authors show existence and uniqueness of a mild solution of equation for and in a suitable set of parameters, denoted (see (8) for details).

In Section 4 we will review some of the above mentioned results from [3], which are needed in our work. Here we just notice that the stochastic integral that appears in (3) is well-defined thanks to fractional Morrey’s inequality ([24, Thm. 2.8.1, Remark 2]) which guarantees the embedding with (in our case we will have ).

It is worth noticing that the concept of virtual solution follows a Zvonkin-type transformation based on heuristic application of Itô’s formula to

. This allows to replace the drift term in (1) with the terms in (3) depending on and . The reader might have noticed that the PDE (4) and the virtual SDE (3) depend on an extra parameter , while the original SDE (1) does not. This is due to a technical step in the proof, that leads to good properties of . However, it is possible to show that the virtual solution is independent of , as shown in [3, Section 3.3] (see also Section 4 below).

3. Main theoretical results

Our numerical scheme for (1) is based on two subsequent approximations. First the distributional coefficient is replaced by a function and then we apply a generalised Euler-Maruyama scheme.

To fix notation, let us consider a function such that the SDE


admits a strong solution in the classical sense (see more detailed assumptions in Proposition 3.1). Let and let us consider an equally-spaced partition of with intervals, that is for . Further, let us define

Then the Euler-Maruyama approximation of the solution is given by


Notice that from the point of view of simulation one uses


are i.i.d. standard Gaussian random variables.

Next we give a preliminary result on the convergence (in law) of the 2-step scheme presented above.

Proposition 3.1.

Let Assumption 1 hold and let be the virtual solution of (1). Consider a sequence in such that for some and is bounded on compacts. Then the approximated stochastic process defined in (6) converges in law to as and . In particular, for every continuous bounded function we have


where the limits are taken in order.


The proof is in two steps.

(Convergence of to ). The fact that converges in law to the virtual solution follows from [3, Proposition 29]. Notice that assumption (i) in [3, Proposition 29] is only used in Step 1 therein to ensure existence of a strong solution of (5) and to guarantee that it is also a virtual solution. These two facts remain true under our assumptions. Indeed [11] ensures existence of a unique strong solution to (5) and [3, Proposition 26, part (ii)] shows that is also a virtual solution.

(Convergence of to ). Now we fix and let . Using [4, Theorem 2.4] we have

for any . The latter implies convergence in law for the process.

Combining the two steps above we obtain (7). ∎

Remark 3.2.

The above result continues to hold, by the same arguments of proof, if we assume (continuous with bounded first derivative). However, from the point of view of numerical simulations the choice of the specific sequence will in general affect the computational efficiency.

In Section 7 below we will construct a sequence that satisfies the hypotheses of Proposition 3.1 by using Haar wavelets, and we will implement it numerically in the special case .

The main theoretical result of the paper is about the rate of convergence of the numerical scheme in the one dimensional case () with respect to a -norm. We first find the rate of convergence of to in terms of the rate of convergence of to (Theorem 3.3). Then, for fixed , we obtain the rate of convergence of the Euler-Maruyama scheme (Theorem 3.5).

A solution to (4) can be found for (see details in Section 4), where


The set is drawn in Figure 1 for the reader’s convenience and it is not empty thanks to Assumption 1. For our proofs below we will need to assume and , where we recall that

relates to the Hölder space which the solution of (4) belongs to. Then, given , we define


Is is now easy to verify that if and only if .

Below we state the result about the convergence rate of to . The proof builds on a number of lemmas and we give it in Section 5.

Theorem 3.3.

Let Assumption 1 hold, let and and let , where is given in (9). Then, for any such that , there exists such that


for all and some .

Remark 3.4.

The constant is found explicitly, see (43).

Next we establish a rate of for the convergence of to as in the Euler-Maruyama scheme. To this aim, we further require that is piecewise Lipschitz and bounded, in line with existing literature. As explained in the Introduction, here we propose a method which is of independent interest and holds for the Euler-Maruyama approximation of any one-dimensional SDE with unitary diffusion coefficient and whose drift is piecewise Lipschitz and globally bounded.

The proof of the next theorem is given in Section 6.

Theorem 3.5.

Let and let be a finite partition of the real line. Fix and assume that is piecewise Lipschitz on the partition intervals, that is for all and . Assume moreover that is uniformly bounded on by . Then, for any and sufficiently large we have


for some constant depending on but independent of .

Remark 3.6.

The constant is found explicitly, see (71).

Combining the theorems above we obtain the next corollary.

Corollary 3.7.

Let Assumption 1 hold, let and and let , where is given in (9). Assume that in as and that for each there exists a finite partition of such that satisfies the assumptions of Theorem 3.5. Fix such that . Then, there exists such that for all and sufficiently large we have


for all .

4. Background material on virtual solutions

As anticipated, the proofs of our main results (Theorem 3.3 and Theorem 3.5) rely upon a few technical lemmas. To set out clearly our arguments and keep the exposition self-contained it is convenient to review some results from [3]. A reader familiar with that material can move on directly to Sections 5 and 6.

As explained in Section 2, a virtual solution for (1) is understood in terms of (3). The latter requires existence and uniqueness of the solution (in a suitable sense) for the Cauchy problem (4). For the numerical scheme illustrated in Section 3 we also need to consider the approximating PDE


where, for each , is an actual function that belongs to a suitable subset of (see Proposition 3.1).

We will now review the arguments that guarantee existence, uniqueness and regularity of the solutions to (4) and (13). First of all we need to restrict (recall (8)).

Figure 1. The set . Figure taken from [3].

Then, [3, Theorem 14] guarantees that, for each there exists a unique solution to (4). Since the time derivative and the second spatial derivative of are not well defined, is a so-called mild solution (for details see, e.g., [6]), and it is obtained as a fixed point in the space equipped with the norm , with sufficiently large.

Using fractional Morrey’s inequality [24, Thm. 2.8.1, Remark 2] it is possible to embed the fractional Sobolev space in smoother spaces. In particular we have that the solution with . Analogously, (13) admits a unique solution (regularity of could be upgraded by virtue of higher regularity of but this will not be needed for our purposes).

Next, [3, Lemma 20] gives useful bounds for the gradient of and . We give a statement which is adapted to our notation111We note that there is a typo in the statement of [3, Lemma 20]. Indeed it can be easily checked from the proof that the condition is not needed therein..

Lemma 4.1.

Let . There exists such that, given any , letting and be the mild solutions in to the corresponding problems (4) and (13), respectively, we have


Furthermore, [3, Lemma 21] also guarantees that


The next result is a refined statement of [3, Lemma 23]. In particular our equation (16) is contained in the final part of the original proof in [3].

Lemma 4.2.

Let and let and be the mild solutions in to (4) and (13), respectively. Then, if in , there is a constant such that


for any that is sufficiently large to guarantee that the denominator above is positive.

For future reference we define


for large enough so that the denominator is positive.

Remark 4.3.

Notice that in Lemma 4.1 we can choose , sufficiently large, so that depends only on and , because in . Then, in Lemma 4.2 we can choose so that the denominator in (16) is positive and (as needed for the fixed point in [3, Thm. 14]).

From now on we will simplify our notation and set , for some sufficiently large so that Lemma 4.1 holds. In order to solve equation (3) and find a virtual solution of (1), one has to transform the SDE (3) into a more standard one. This is achieved by setting where


Notice that thanks to (15). Moreover by Lemma 4.1 is invertible for each fixed (see also [3]), with its inverse denoted by


By Lemma 4.1 is 2-Lipschitz, uniformly in . Then, solving (3) is equivalent to solving the standard SDE for below


where . Existence of a weak solution for (20) is guaranteed by [23, Theorem 10.2.2] since its coefficients and are bounded continuous with uniformly non-degenerate (see [3, Proposition 27] for details).

Likewise, letting , and , the analogue of (20) for the approximated SDE (5) is given by a SDE for . That is


Moreover, is 2-Lipschitz, uniformly in , by Lemma 4.1.

Finally, we recall that in [3, Section 3.3] the authors prove that the virtual solution is independent of .

Remark 4.4.

Notice that in one dimension () both equations (20) and (21) admit a unique strong solution if because the diffusion coefficient is -Hölder continuous, whilst the drift is Lipschitz. This result is used in the proof of Theorem 3.3

to justify the use of the same Brownian motion when estimating


5. Convergence rate of

In this section we show our first main result, which we anticipated in Theorem 3.3. It turns out that in order to show the convergence rate of to stated in Theorem 3.3 we must provide an upper bound for the local time at zero of . Recall that for any real-valued continuous semi-martingale , the local time is defined as


for all . Now we derive a bound on (22) that will be needed later on.

Lemma 5.1.

For any and any real-valued, continuous semi-martingale we have


Figure 2. The function .

For and we define (see Figure 2)

Straightforward calculations allow to show that and it is semi-concave, in the sense that is concave. Moreover, we have


Now, an application of Itô-Tanaka formula gives


where denotes the left/right limit of the derivative at zero. Rearranging terms, taking expectations and using (24)–(27) gives (23). ∎

The next lemma controls the approximation error between and . This result holds for any .

Lemma 5.2.

Let and let be the mild solutions to (4) and (13), respectively. Then, for and as in Remark 4.3, and all we have




given in (17) and .


Let . Since , by the fractional Morrey’s inequality we have and we can find such that for all it holds


Then, recalling (2) and (16) we easily obtain


Combining (31) and (32) gives (29).

Now we provide a bound on the difference . This result holds for any .

Lemma 5.3.

Take and as in Remark 4.3 and as in Lemma 5.2. Then we have


Recall that was defined in (18). For any , denoting the scalar product in , we have

where the first inequality uses Cauchy-Schwartz and the final one Lemma 4.1.

Now, taking and in the above inequality gives

The latter implies

where the final equality uses . By definition of and and (29) we also obtain