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 , 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  which links the SDE in (1
) to a class of partial differential equations (PDEs) with distributional drift studied in. 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  and , 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 , where a (strong) -rate of is obtained for any , and  where analogous results are obtained in a multi-dimensional setting with respect to a (strong) -norm (here ‘strong’ refers to the fact that  and  consider the -norm on the space of continuous paths, under expectation). While we were completing our work another paper has appeared () 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 , ,  and  (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  where a -rate of convergence of is obtained for (possibly degenerate) multi-dimensional SDEs. Notice that in  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  to find a (strong) -rate of convergence of ; differently from , the convergence in  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  (multi-dimensional setting, rate ) and  (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  proves strong convergence in but with no rate and  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 , 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 ). However, we would like to mention that weak convergence with rate up to is obtained in , and  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  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 .
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 fieldssuch 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  and make the following assumption.
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. Hereis the mild solution in of the following parabolic Kolmogorov-type PDE
In Section 4 we will review some of the above mentioned results from , 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.
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  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). ∎
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.
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).
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.
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.
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 .
The constant is found explicitly, see (71).
Combining the theorems above we obtain the next corollary.
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 . 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).
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., ), 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..
Furthermore, [3, Lemma 21] also guarantees that
For future reference we define
for large enough so that the denominator is positive.
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
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).
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 .
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.
For any and any real-valued, continuous semi-martingale we have
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
The next lemma controls the approximation error between and . This result holds for any .
Now we provide a bound on the difference . This result holds for any .
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