The Euler scheme for stochastic differential equations with discontinuous drift coefficient: A numerical study of the convergence rate

The Euler scheme is one of the standard schemes to obtain numerical approximations of stochastic differential equations (SDEs). Its convergence properties are well-known in the case of globally Lipschitz continuous coefficients. However, in many situations, relevant systems do not show a smooth behavior, which results in SDE models with discontinuous drift coefficient. In this work, we will analyze the long time properties of the Euler scheme applied to SDEs with a piecewise constant drift and a constant diffusion coefficient and carry out intensive numerical tests for its convergence properties. We will emphasize on numerical convergence rates and analyze how they depend on properties of the drift coefficient and the initial value. We will also give theoretical interpretations of some of the arising phenomena. For application purposes, we will study a rank-based stock market model describing the evolution of the capital distribution within the market and provide theoretical as well as numerical results on the long time ranking behavior.

The existence and uniqueness of solutions of SDEs in the standard case, i.e. the case of sufficiently smooth coefficients, is well understood [18]. However, the standard theory on SDEs does not apply anymore in case of a discontinuous drift coefficient, e.g. a piecewise constant drift coefficient, and a special theory is needed to address the question of existence and uniqueness of solutions of such SDEs [19,40,41]. The same is true for the numerical analysis: The convergence behavior of approximation schemes needs to be reconsidered and "research on numerical methods for SDEs with irregular coefficients is highly active." ([25, p. 2]). In the case of a sufficiently smooth drift and a constant diffusion coefficient, the exact strong rate of convergence is 1 for the Euler scheme, see [5,20]. However, so far -to the best of our knowledge -there is no such result in the case of a discontinuous, e.g. piecewise constant, drift coefficient. Although results on the theoretical convergence rates have been established (see e.g. [25,32]), the sharpness of these error bounds is not clear yet, see for example the simulation study in [26].
In this work, we will focus on numerical approximations of SDEs in the presence of a piecewise constant drift and a constant diffusion coefficient. We will provide theoretical considerations on the long time behavior of approximated SDE solutions based on results from the theory of ergodic Markov chains, and we will provide further insight into the numerical behavior of these approximation schemes, in particular the Euler scheme, by analyzing the numerical convergence rates based on a reference solution. The speed of convergence heavily depends on the initial value and properties of the drift coefficient, e.g. drift direction or jump height. Our tests reveal that for a special class of drift coefficients the convergence rates are higher and independent of initial conditions due to the ergodicity of the Euler scheme and the underlying SDE. We also use the Euler scheme to verify the long time behavior of a rank-based stock market model [4], a prominent model in finance to describe the evolution of the capital distribution within the market.
The remainder of this manuscript is as follows: In Section 2, we will introduce some theoretical and numerical basics and establish the ergodicity of the Euler approximations in the case of an appropriate, piecewise constant drift coefficient. In Section 3, we will discuss numerical convergence properties and further findings of several numerical tests. We will conclude this work in Section 4 with the application from mathematical finance mentioned above, where SDEs of discontinuous type naturally arise.

Problem Description
In this section, we will introduce our basic setting, i.e. the type of SDE, we are interested in, and some basic terms for the numerical tests. Besides the Euler scheme and its long time properties in our setting, we will also briefly discuss the applicability and performance of some other numerical schemes.

The Equation.
In this manuscript, we will consider time-homogeneous SDEs with piecewise constant drift coefficient and additive noise: Here, we have s ∈ N, α j , σ, ξ ∈ R and disjoint (possibly infinitely many) intervals B j ⊂ R for all 1 ≤ j ≤ s and (W t ) t∈[0,T ] is a one-dimensional Brownian motion.
The existence and uniqueness of solutions to this type of SDEs is guaranteed by results of [40,41] and [19]. In [40], conditions on the drift and diffusion coefficient are derived under which the corresponding SDE has a unique strong solution. As emphasized therein, those conditions are in particular fulfilled for a bounded drift coefficient and a constant diffusion coefficient. Thus, the existence and uniqueness of a strong solution for SDEs of type (1) is ensured.
For the numerical analysis of SDEs with discontinuous drift and/or diffusion coefficient, the situation is more involved. In this manuscript, we will focus on the strong convergence rate of the Euler scheme, which, for a general SDE where f and g are such that a unique strong solution exists, is given by The underlying time discretization of the time interval [0, T ] is 0 = t 0 < t 1 < · · · < t n = T with corresponding step size ∆ := T n , where n + 1 is the number of grid points. While its behavior is well-known for SDEs with Lipschitz continuous coefficients f and g, much less is known in more general cases, even for SDEs with additive noise and a piecewise constant drift coefficient. The first contribution in this area is -up to the best of our knowledge -the work [12], where almost sure convergence of the Euler scheme has been established in the case of a one-sided Lipschitz drift coefficient, a locally Lipschitz diffusion coefficient and the existence of a Lyapunov function for the SDE. The results of [13] give strong convergence of the Euler scheme for SDEs with additive noise in the case of a discontinuous, but monotone drift coefficient, while [38] establishes the almost sure and strong convergence of the Euler scheme for SDEs with additive noise and drift of the form f (x) = −sign(x). The most recent contributions with respect to strong approximations of SDEs with discontinuous drift coefficient are a series of articles by Ngo and Taguchi [32,33,31] and Leobacher and Szölgyeny [25,26,27], respectively. The weak approximation of SDEs with discontinuous coefficients has been studied in [21], where an Euler-type scheme based on an SDE with mollified drift coefficient is analyzed.
In the case of SDE (1), the latest results on the strong convergence rate of the Euler scheme for the approximation of X T , i.e. the solution at time T , are: • an L 1 -convergence order of 1/2 − ε with ε > 0 arbitrarily small in [32], and • an L 2 -convergence order 1/5 in [25].
Note that none of these results has been tailor made for SDE (1) and the Euler scheme (2). In [25], the L 2 -convergence order 1/5 has been established for multidimensional SDEs with a globally Lipschitz diffusion coefficient and a piecewise globally Lipschitz drift coefficient, where the exceptional set satisfies a geometric smoothness condition. Dai and Taguchi analyze in [32] the Euler scheme for scalar SDEs with discontinuous drift coefficient, which can be appropriately approximated by smooth functions, and bounded Hölder continuous elliptic diffusion coefficients.
For a better comparison, note that in the standard setting of an SDE with additive noise, where the drift coefficient is sufficiently smooth, the Euler scheme has an exact strong convergence order of 1, see e.g. [5] and [20, p. 350f].
So to summarize: The Euler scheme for our non-standard setting of SDE (1) has at least L 2convergence order 1/4 − ε. However, already a small simulation study as in [26] reveals that this order is not sharp.

2.2.
Simulation studies and empirical convergence rates. As already mentioned, we are interested in the strong convergence rate of the Euler scheme. More precisely, we will focus on the root mean squared error (RMSE) at time T for the Euler scheme (1) with step size ∆ = T /n, i.e.
Since an explicit form of X T is unknown in general, we need to replace X T in our simulation studies by a numerical reference solution X num T , which is computed by the Euler scheme for an extremely small step size ∆ = T /N with a very large number of N + 1 grid points such that this approximation can be considered close enough to the true solution. Moreover, also the expectation E|X num T − x expE n | 2 is not known explicitly, so we will approximate this expectation by the empirical RMSE Note that N has to be chosen sufficiently large to generate the numerical reference solution and to avoid oscillations in e emp (n), which might occur if N and n are close. The number of repetitions M should also be large enough to have a good approximation of the expectation, i.e. a small Monte Carlo error.

Other schemes.
A natural idea is of course to consider other schemes than the explicit Euler scheme and to compare them in our simulation studies.
2.3.1. The implicit Euler scheme. Implicit schemes have good stability properties, thus, they are a natural choice to consider. For an SDE with additive noise, where the drift coefficient is sufficiently smooth, the implicit Euler scheme has strong convergence order 1 (see e.g. [1] and [30]).
However, for SDEs of type (1), already the implicit Euler scheme is not well defined. To see this, consider the SDE . . , n−1, requires to solve, for fixed but arbitrary z ∈ R, the equation with respect to y ∈ R. This equation does not possess a solution if z ∈ (−α 1 ∆, −α 2 ∆), and hence an implicit Euler scheme is not well defined in this setting.
2.3.2. The Heun scheme. The Heun scheme is another scheme with strong order one for SDEs with additive noise under appropriate smoothness conditions on the drift coefficient. Adapted from [20, p. 373] for SDEs of type (1), it is defined by For a closer look at the behaviour of this scheme at a discontinuity assume that the drift coefficient is given by a(x) = ±sign(x). An increment of the Heun scheme with x Heun k = x is then given by So if no drift change occurs in the Euler step x + a(x)∆ + σ(W (k+1)∆ − W k∆ ), a Heun step and an Euler step coincide. However, if a drift change occurs in the Euler step, the Heun step reads as , it approximates the drift coefficient by zero and its dynamics are purely diffusion-based in this case.

2.3.3.
A Wagner-Platen type scheme. A strong order 1.5-scheme for SDEs with smooth drift coefficient and additive noise is given by a Wagner-Platen type scheme (see e.g. [20, p. 383]), which reads in our setting as Now, we look again at the case of a drift coefficient given by a(x) = ±sign(x). For a Wagner-Platen step with x Pla k = x, it depends now on whether have the same sign or not. If this condition is fulfilled, i.e., if x is sufficiently far away from the discontinuity, then, a Wagner-Platen step and an Euler step coincide. If the latter condition is not satisfied, then we have the dynamics So also here, the diffusive dynamic dominates the scheme when taking values close to the discontinuity.

Ergodicity and Stability of the Euler Scheme.
We will now address the long time properties of the Euler scheme based on results from the theory of ergodic Markov chains. For simplicity, we consider here a special case of SDE (1), namely and assume that i.e. a drift coefficient, which is pointing towards zero. Clearly, we have i.e., on average the solution is moving inwards. Moreover, following e.g. chapter 6 in [10], this SDE admits a unique invariant distribution with Lebesgue density and the law of large numbers It will turn out that the explicit Euler scheme with will recover these properties. (Here we also indicate the dependence on the initial value ξ in our notation.) The Euler scheme (8) corresponds to a time homogenous Markov chain with transition kernel and satisfies the discrete counterpart to (5), i.e.
Now, we will proof the of the existence of a unique stationary distribution for the Euler scheme. In particular, due to the discontinuity at zero, the following Proposition 1 is not covered by the standard references as e.g. [28] and [35] for Euler-type discretizations of ergodic SDEs. Note that the long time properties of (8) have also been heuristically studied in [37].
However, we can easily verify that V (x) = e τ |x| , x ∈ R, is an appropriate Lyapunov function for the above Markov chain, if τ > 0 is sufficiently small. This is a direct consequence of the well known form of the moment generating function for the folded normal distribution, i.e.
Note that the limit distribution is independent of the initial value, as for the underlying SDE.
Finally, an ergodic Theorem as e.g. Corollary 2.5 in [6] yields also the discrete counterpart to the law of large numbers (7): We have

Simulation Studies
This section is concerned with the numerical investigation of SDEs of type (1). For the remainder, we will choose T = 1, M = 10 5 , N = 2 14 and n = 2ñ withñ ∈ {4, . . . , 10} (unless otherwise mentioned). We then calculate the corresponding Euler approximation and the empirical RMSE e emp (n). For simplicity, we omit the upper index of the numerical approximation indicating that the approximation is based on the Euler scheme. The empirical convergence rate is given by the negative slope of the regression line, which we obtain when plottingñ versus log 2 e emp (n) . Here, we will focus on two types of drift coefficients: inward and outward pointing drift coefficients.
and outward pointing, if there exists x * ∈ R such that Our numerical investigations are based on several additional key characteristics: We consider the average number of drift changes. As the Euler scheme for SDE (1) is exact up to the first drift change, another quantity of interest is the number of paths with at least one drift change. To get further insight whether some paths are really far away from the true solution, we measure the largest error that occurs within the considered time interval (not necessarily in the end). Besides the error sizes themselves, it is interesting to see what proportion of errors at final time T is large, medium, or small and how this distribution of error sizes depends on the step size. Furthermore, we analyze the evolution of the error over time for a fixed step size. To underline the influence of the drift direction towards or away from the discontinuity, we generate plots of several solution sample paths. We will see that the rates of convergence heavily depend on whether the drift coefficient is inward or outward pointing. Whereas for the latter one, there is a dependency on the initial value of the SDE, rates in case of an inward pointing drift coefficient seem to be independent of the initial value, corresponding to Proposition 1. In addition, we analyze how the jump height (difference in drift values) influences the empirical convergence rate.
As representatives of the class of SDEs (1), we consider here the SDEs given in Table 1.

Drift coefficient
Corresponding In the remainder of this chapter, we will present and discuss some key results of the simulation studies.
3.1. Key results. The empirical convergence rates obtained by the Euler scheme are given in Table 3.1 (outward pointing drift coefficients highlighted in light gray, the discontinuity in gray):  The Euler scheme converges for every analyzed SDE and for every tested initial value with a better order than the theoretically guaranteed 1/4. However, our results show that • in general, we loose convergence order one, which the Euler scheme has under standard assumptions for SDEs with additive noise (and we even have a convergence order less than 1/2 in some cases), • and that a crucial factor is whether the drift coefficient is inward or outward pointing.
Furthermore, our numerical tests show that neither using the Heun scheme nor using the Platen scheme yields a different picture. In particular, convergence rates do not improve significantly, and the schemes do not yield a better resolution of the discontinuity (see Tables 3 and 4).

Drift direction and initial value.
For an outward pointing drift coefficient, the convergence order even seems to depend on the initial value and the spectrum of orders obtained for different initial values is very broad with values between 0.25 and infinity (see Table 3.1). Such a behavior of the convergence rates dependent on the initial values has been unknown so far in the literature. On the other hand, for an inward pointing drift coefficient, the convergence order seems to be independent of the initial value and the spectrum of orders obtained for different initial values and inward pointing drift coefficients is tight with values between 0.80 and 0.91 (see Table 3.1). The stability of the estimates is due to the ergodicity of the SDE and the Euler scheme in this case, see Subsection 2.4. The geometric convergence speed in Proposition 1 explains why the numerical tests for inward pointing drift coefficients yield such stable estimates, independently of the initial value: X num T and x n are, for a sufficiently large number of grid points n + 1, close to their unique stationary distributions, which stabilizes the Monte-Carlo estimates.
For the above equations, the monotonicity of the drift coefficient is directly related to the number of drift changes. An inward pointing drift coefficient results in many drift changes, while in the case of an outward pointing drift coefficient, only few drift changes occur. We can further observe that: (i) when starting away from the discontinuity, rates for outward pointing drift coefficients seem to be better than for inward ones; (ii) when starting close to the discontinuity, results suggest that outward pointing drift coefficients imply worse convergence rates than inward ones. So, in the latter case we obtain a positive correlation between the number of drift changes and the convergence rate, which implies that frequent drift changes are not necessarily bad for the quality of the approximation -quite the contrary seems to apply, which is surprising at first glance.
Hence, the type of monotonicity of the drift coefficient is of great importance. Intuitively, an inward pointing drift coefficient should lead to many drift changes, which suggests that individual drift changes are not of great importance. An outward pointing drift coefficient on the other hand pushes the solution away from the discontinuity implying a low number of drift changes.

Jump height.
The intensity of the effects related to inward and outward pointing drift coefficients depends on the jump height, i.e. the distance between assigned drift values. In case of elementary minus34, this distance amounts to 7 whereas it is 1.6 in case of elementary minus0.6 1. The empirical convergence rates in Table 3.1 show: The higher the jump height, the more pronounced are the effects described in Subsection 3.2. Exemplary, there is a difference of 0.8 in the empirical convergence rates for elementary minus34 for initial values 0 and 1 whereas this difference is only 0.06 for elementary minus0.6 1. This phenomenon is related to a scaling property. By enlarging the drift value, the influence of the diffusive part of the SDE is weakened: Consider e.g. the SDE dX t = α sgn(X t )dt + dW t with α ≥ 1. Using the new variable Y t = 1 α X t we have the dynamics with a reduced diffusion coefficient.

3.4.
Case study of an inward versus outward pointing drift coefficient. In this Subsection, we will analyze the pattern described in 3.2 in more detail, exemplary for the drift coefficients elementary4minus3 and elementary minus34.

3.4.1.
Drift changes. Figure 1 shows the average number of drift changes for both coefficients. The behavior goes along with the intuitive understanding described above. Hereñ is the exponent of the dyadic step size ∆ = 2 −ñ . Note that for step sizes 2 −4 to 2 −8 and elementary minus34 the number of drift changes stays below 2.   Figure 2 shows 100 sample paths of the numerical reference solution (∆ = 2 −14 ). The black line represents the discontinuity in the drift coefficient.
In the situation of Figure 2 (b), where the solution drifts away from the discontinuity, it is of tremendous importance whether a drift change is captured by the approximation or not: the solution does not stay close to the discontinuity and thus, there are not many chances for a drift correction to take place, see Figure 3. For the SDE with α 1 < 0 < α 2 and ξ > 0 the conditional probability p(ξ, θ, ∆) that the exact solution changes its drift over [0, ∆] given that the approximation x 1 at t = ∆ has value θ ≥ 0 (and thus has not changed its drift) satisfies see e.g. [11], page 169, so the (conditional) probability of missing drift changes is not negligible and even close to one for small ξ or θ.  Largest error. The latter observation is also reflected in the largest distance for 10 4 sample paths between the approximation based on step size 2 −10 and the numerical reference solution, see Table 5. The largest distances amount to 1.271 for elementary4minus3 and 4.508 for elementary minus34.  We can extract from Figure 4 at least two features: (i) The error stays constant or even decreases over time for elementary4minus3 -in contrast to a strong error accumulation over time for elementary minus34. (Note that the ordinate has a base-2 log scale.) (ii) In the inward pointing drift coefficient case, the error is by several magnitudes smaller than for an outward pointing drift coefficient. This illustrates again the stabilizing effect of an inward pointing drift coefficient and the importance of capturing the first drift changes correctly in case of an outward pointing drift coefficient.
3.4.5. Distribution of error sizes. Besides the empirical RMSE itself, the empirical distribution of the errors in t = T is of interest. The error at final time T is quantified by |x N − x n | for step size ∆ = 2 −ñ . The histograms in Figure 5 are based on M = 10 4 simulations for different step sizes and highlight again the different magnitudes of the empirical RMSE (abscissa with a base-2 logarithm scale). Another feature, which we can extract from the histograms, is a non-negligible part of simulated paths with an error of machine accuracy size for elementary minus34. We will discuss this feature in more detail in the next Subsection.

3.5.
Rare events and goodness of the regression fit. In case of an outward pointing drift coefficient and an initial value far away from the discontinuity, the empirical RMSE and the linear regression estimates become unreliable or at least questionable. Here, in the underlying SDE only very few drift changes occur (if at all), which means that, if the step size of the Euler scheme is significantly small, these changes are captured and the error drops drastically. Figure 6 illustrates this by comparing the regressions for an initial value ξ = 0 away from the discontinuity in 1.4 and an initial value ξ = 1, which is closer to the discontinuity. Note that the Euler scheme for (1) is always exact up to the time of the first drift change.
Moreover, for an outward pointing drift coefficient, the Euler scheme and the exact solution coincide with high probability, which explains e.g. the errors close to machine accuracy for the drift coefficient sign and the initial value ξ = 5. Note that in this setting, the number of paths with at least one drift change is even zero over all saved 10 4 solution paths.   To explain this phenomenon, consider again the SDE with α 1 < 0 < α 2 . An application of formula (5.13) in chapter 3.5.C in [18] gives Note that an initial value ξ = 0 is not a restriction as we analyze the case of an initial value far away from the discontinuity. So, for drift values −α 1 = α 2 = 1, and an initial value ξ = 5, the Euler scheme is exact with a probability of at least 1 − e −10 ≈ 0.99995460 . . .
To summarize: While Monte Carlo simulations for testing convergence rates are possibly unreliable here, the error of the Euler scheme can be easily controlled. For SDE (13) the above formula and give that , independently of the step size. So even a simple one-step Euler scheme can have a very small error for such SDEs.

The Euler scheme for the Atlas model
In this Section, we will use the Euler scheme to simulate the so-called Atlas model, which is a particular first-order market model [4]. In such models, the asset dynamics depend on the size (measured in terms of market capitalization) of the corresponding firm, which results in an SDE model with discontinuous coefficients.
Consider now stocks for which the market capitalizations are given by X 1 , . . . , X d , where the index i ∈ {1, 2, . . . , d} indicates the name of the firm, and that follow the dynamics Here, W 1 , ..., W d are independent Brownian motions and the growth rates γ i : [0, ∞) → R and volatilities σ i : [0, ∞) → (0, ∞) are given by The ranks r i (t) for the stock X i (t) at time t arise from the reverse order-statistics: Ties in the ranking are resolved by giving the firm with a lower index i the better ranking. So, in such a model the k-th largest firm is assigned a growth rate of γ + g k and a volatility of σ k over the whole time horizon.
According to [4], the simplest among the first-order models is the so-called Atlas model, which was introduced in [7, Ex. 5.3.3]. Within the setting of (17) and (18), choosing leads to the Atlas model. Here, only the smallest stock in the market -called the Atlas stockhas a nonzero but positive growth rate (for its log-dynamics).
i.e., all stocks in the market asymptotically spent at each rank approximately the same amount of time. Similar ergodic relations also hold for general first-order market models.

4.2.
Numerical results. For simulations of the Atlas and general first-order models one has to rely on discretization schemes such as the Euler method. In this Subsection, we test whether the Euler scheme is able to recover the long time behavior (22), i.e., whether the discrete occupation rates 1 T where r i is the discretized counterpart of (19) based on the Euler scheme and T /∆ ∈ N, converge to the analytical value.
Here, we consider a three-dimensional model with initial log-capitalizations Y (0) = [3.4, 4.1, 5.7] and Y (0) = [1.2, 3.5, 10.8], γ = 0.1 as market drift and σ = 0.09 as market volatility 1 . Table  6 presents the discrete occupation rates (averaged over M = 10 3 repetitions) for ∆ = 2 −14 and different values of T as well as the sum of the squared deviations from the analytical asymptotic occupation rate. As hoped, the discrete occupation rates converge to the analytical asymptotic occupation rate of 1/d = 1/3 with increasing time horizon.
Furthermore, results suggest that similar initial capitalizations imply that the numerical values are closer to the analytical result already for shorter time horizons, which coincides with the intuitive understanding. We also simulated the above scenarios with ∆ = 2 −10 instead of ∆ = 2 −14 : all occupation times where equal with an accuracy of four digits and one third of the 90 occupation rates differed in the fifth digit. This suggests that -as soon as the step size is small enough -a further refinement of the step size is no longer beneficial and the crucial simulation parameter is T , the endpoint of the considered time horizon.

Conclusion and Outlook
We have seen that the numerical approximation of solutions of SDEs with discontinuous drift coefficients is a challenging task, where several particularities arise. We were able to identify two main classes of discontinuous drift coefficients: outward and inward pointing drift coefficients. For the latter class, we analyzed stability properties. It turned out that the main difficulty in the numerical approximation is how to appropriately capture drift changes. We tested two higherorder numerical schemes, that are frequently used in a setting where coefficients are sufficiently smooth. However, both schemes did not lead to an improved convergence behavior. A remedy might be the usage of adaptive numerical schemes. This has already been done successfully for SDEs with additive noise and a smooth drift coefficient (see e.g. [9,14,15]). Subject of future work will be the definition and implementation of an adaptive step size control close to the discontinuity. We expect to thereby get a better resolution of the discontinuity and, thus, to improve the overall quality of the approximation.