{'fraction_non_alphanumeric': 0.05646065316493536, 'fraction_numerical': 0.01939316941932591, 'mean_word_length': 4.166602316602317, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Positivity preserving logarithmic Euler-Maruyama type scheme for stochastic differential equations

In this paper, we propose a class of explicit positivity preserving numerical methods for general stochastic differential equations which have positive solutions. Namely, all the numerical solutions are positive. Under some reasonable conditions, we obtain the convergence and the convergence rate results for these methods. The main difficulty is to obtain the strong convergence and the convergence rate for stochastic differential equations whose coefficients are of exponential growth. Some numerical experiments are provided to illustrate the theoretical results for our schemes.


Introduction
In this paper we consider the numerical solution of the following scalar stochastic differential equation (SDE) dx(t) = b(x(t))dt + σ(x(t))dW(t), (1.1) where b and σ are locally Lipschitz continuous functions and the initial data x(0) = x 0 > 0. In many practical problems the quantities (such as stock price in finance, the population affected in epidemic and so forth) that we model take positive values. When we use such a stochastic differential equation to model positive quantity, first we need to make sure that the solution to the equation is positive. This property of the equation can be guaranteed by the famous Feller's (non-explosion) criteria (see next section for more discussion). Now the problem is that if the stochastic differential equation (1.1) has a positive solution, how to construct its numerical approximations so that they are all positive? There have been some works along this line. On one hand, we know many existing numerical algorithms cannot maintain the property of positiveness for the approximate solutions. In fact, in [1], authors proved that the approximate solutions from the Euler-Maruyama method cannot remain positive for any (one dimensional) SDE. On the other hand, there are several existing works on preservation of positivity. For instance, in [2], by using the Lamperti transformation, authors transformed the original SDE whose diffusion coefficient is strictly positive in a certain domain into the following SDE dx(t) = F(x(t))dt + c 0 dW(t).
The transformed SDE was discretized by a drift implicit Euler method and was transformed back so that an approximate solution is inside the domain of the original solution under appropriate assumptions. The convergence rate of their scheme was proved to be 1.0 and their scheme can be applied to practical models, such as the Cox-Ingersoll-Ross (CIR) model and the constant elasticity of variance (CEV) model, etc. However, sometimes an additional nonlinear equation needs to be solved in order to obtain the approximate solution, which is time consuming to some extent. A class of balanced implicit methods which can preserve positivity under an appropriate choice of step size and weight functions was shown in [3], but the convergence results were derived under the linear growth condition which is restrictive for some practical models. Since its importance in finance the following CIR model has received particular attention. In the case of unattainable boundary (2κλ ≥ θ 2 ), based on the Lamperti transformation and the semi-implicit Euler method, the following numerical scheme was proposed in [4] for this model: It is obvious that the approximate solutions are all positive. The convergence order of this scheme on some finite interval was obtained to be 1 2 in [5], and was improved to 1.0 in [6]. This idea was applied to SDEs (driven by Brownian motions, see [2]) and SDEs with poisson jumps (see [7]) which take positive values. Additionally, there are various modifications of explicit Euler-Maruyama method. For example, for the CIR process, the following explicit method was proposed in [8], and in [9] this scheme was proved to have strong convergent rate 0.5 in L p sense when 2κλ θ 2 > 3, 1 ≤ p < λ − 1. A symmetrized Euler method with L p convergence rate of 0.5 was obtained in [10] under a more restricted condition on the parameters 2κλ Some other works can be found in [11][12][13] and references therein. However, for general stochastic differential equations not much study to construct positive numerical approximations has been done. In this work we shall construct positive numerical approximations for general equations. Namely, once the equation (1.1) satisfies the Feller's condition so that its solution is positive, we shall construct its numerical approximations which are all positive themselves. More precisely, we shall propose an explicit logarithmic Euler-Maruyama type scheme such that it can keep positive unconditionally. This paper is organized as follows. In Section 2, we give a general idea of our method. In Section 3, we consider the almost sure convergence of a class of positivity preserving numerical methods obtained by combining the logarithmic transformation and some classical stochastic methods, such as the Euler-Maruyama method and the stochastic theta method. In Section 4, we prove the exponentially integrability of both the exact solution and numerical solution from the logarithmic truncated Euler-Maruyama method (logTE). Moreover, we obtain the strong convergence results of the logTE method under Assumption 4.2 with m < 2. Since there are models (for example, the CIR process) which do not satisfy this assumption, we weaken Assumption 4.2 to Assumption 4.16, and prove that the logTE numerical approximation is also strongly convergent. Next, in Section 5, we analyze the convergence rate of the logTE method under a slightly stronger Assumption 5.1. Some numerical experiments are carried out in Section 6 to illustrate our theoretical conclusions for some popular models. Finally, we make a brief conclusion in Section 7.

General idea
Since we consider only positive solution to the equation (1.1), we assume that the coefficients b and σ are real valued functions defined on the positive real line R + := (0, ∞) and we also assume they are locally Lipschitz continuous. From these coefficients, we introduce the scale function Assumption 2.1. We assume that the following Feller's condition holds true: According to Feller's criteria (see [14]) the condition (2.1) is equivalent to the fact that the SDE (1.1) has a unique strong solution x(t), t ≥ 0, and its sample path is contained in R + , i.e.
Now assume the condition (2.1) holds true, the solution is then positive for all time t ≥ 0. How we construct numerical approximations so that they are also positive? In the introduction we briefly mentioned some earlier works about positivity preserving numerical schemes. In this work, we want to find general numerical schemes to construct positive approximate solutions under the Feller's condition criteria (2.1). Our general idea is first to take the following logarithmic transformation The application of the Itô formula with y(t) = ln x(t) yields the following transformed stochastic differential equation for y(t): where f (y) = e −y b(e y ) − 1 2 e −2y σ 2 (e y ) , g(y) = e −y σ(e y ) .

(2.4)
Then we find a numerical solutionỹ(t) for (2.3). Andx(t) = e˜y (t) will be a positive approximation of the original solution x(t).
Then the nonexplosiveness of the solution y(t) can still be explained by the Feller's condition, i.e. Since y = ln x, integration by substitution yields Hence, we can derive that the solution of SDE (2.3) will stay in (−∞, ∞) with probability 1. Namely, with probability 1, the solution y(t) to (2.3) has no finite time explosion. Since the SDE (2.3) takes value in the whole space, we do not need to worry whether the numerical approximation is contained in R + (or a sub-domain of R) or not. However, after the logarithmic transformation, the coefficients of the transformed equation (2.3) are usually of super-linear growth or even of exponential growth. This brings new difficulty to the construction of numerical scheme. Take the CIR model (1.2) as an illustrating example, which was commonly used in financial and biological mathematics. Using the Feller condition, we can derive that when κλ ≥ θ 2 2 , x 0 > 0, the CIR process has a (unique) positive solution, see [5,15]. The use of the logarithmic transformation y(t) = ln x(t) will give The coefficients of the above equation increase exponentially (as y goes to −∞), so the classical strong convergence theories based on the global Lipschitz condition or polynomial growth condition are no longer applicable to the above equation.

Almost sure convergence
First, let us introduce some known conclusion about the almost sure convergence of some classical stochastic methods. To simplify the presentation, we fixed a finite time interval [0, T ] for any given fixed T > 0. Consider the partition π : 0 = t 0 < t 1 < · · · < t N = T , where N is a positive integer and denote ∆ = T/N as the uniform step size of the interval [0, T ].
There are several papers considered the almost sure convergence of stochastic methods for SDE with non-global Lipschitzian coefficients, see [16,17] and references therein. In [18], author proved the almost sure convergence of a family of numerical methods, which cover a number of classical stochastic schemes. The methods has the following form where δU i = U(t i+1 ) − U(t i ), U = Wor Z, and other notations are explained in the following assumption.
(4) For any positive integer p, it holds that Then, we can get the following theorem.
Theorem 3.3. If Assumption 2.1 holds and the coefficients of approximation process (3.1) satisfy Assumption 3.1, then x π ∆ (t i ) = e y π ∆ (t i ) converges almost surely to the solution of (1.1). More precisely, for any γ ∈ (0, 1 2 ), there exists a finite random variableη γ,T such that sup Proof. The conclusion of almost sure convergence for this theorem can be obtained by combining the mean value theorem and Lemma 3.2.
Example 3.4. Consider the Euler-Maruyama scheme According to analysis of [18], we know that Φ ∆ (x) is locally Lipschitz if f is. As a result, for any R > 0, Φ ∆ can be inverted on a ball B R = {x :| x |≤ R} for sufficiently small ∆, and Φ −1 ∆ is Lipschitz continuous on Φ ∆ (B R ) and , a.s. with order γ for any γ ∈ (0, 1/2). Besides, we also have y IE , as ∆ goes to 0. Therefore, the stochastic theta approximation y IE ∆ (t i ) will converge almost surely to y(t i ), and x ∆ (t i ) = e y IE ∆ (t i ) will converge almost surely to exact solution x(t i ) of the original SDE (1.1).
Remark 3.6. The random constantη γ,T in inequality (3.2) may not be integrable. The path-wise convergence cannot imply the strong convergence or weak convergence. Therefore, in the next section, we shall study another kind of convergence: the convergence in L p sense.

Strong convergence analysis
In the precedent section, we analyze the almost sure convergence of some numerical schemes. Path-wise convergence plays an important role in analysis of random dynamical systems [19], while sometimes we need more. For example, in mathematical finance, the pricing of complex derivative products is one of the main tasks, which means we have to compute Eφ(x), whereφ : C([0, T ]; R) → R represents the discounted payoff of the derivative. Weak convergence is usually enough for pricing purposes, but to minimize the computational complexity strong convergence results are required when using multilevel Monte Carlo schemes [20]. But the integrability of random constants in inequality (3.2) is still unknown, and it seems hard to analyze these constants. Ifφ is continuous and bounded, then Eφ(x ∆ (t)) → Eφ(x(t)) as ∆ → 0. However, since the p-th power function needed for the L p convergence is not bounded, the path-wise convergence results cannot imply strong or weak convergence. For example, in [21], in terms of the Euler scheme (3.3), authors prove the following lemma.
for all | y |≥ C. Then, there exists a constant c ∈ (1, ∞) and a sequence of nonempty events We already know that the Euler-Maruyama approximation converges almost surely to exact solution if coefficients are Local Lipschitzian. The above lemma states that the Euler-Maruyama approximation may diverge in both strong and weak L p sense.
Since we are concerned with the L p convergence, we need to deal with the exponential integrability of the numerical solutions from our transformation. Therefore, the existing strong convergence conclusions (see, [22,23] and references there in) are not applicable to (2.5) because of its exponentially growing coefficients and because of the additional requirement of the exponential integrability. We shall focus on the Euler-Maruyama method, which is favored in practice because of its simplicity. Even for this simplest scheme the exponential growth causes significant difficulty.
Hence, we need to introduce some additional structure so that the Euler-Maruyama method will converge to solutions of SDEs with exponentially growing coefficients. Our idea is inspired by [24,25]. First, let us make some assumptions on the coefficients of SDE (2.3).
Assumption 4.2. We assume that the coefficients of SDE (2.3) are locally Lipschitzian and satisfy the following condition: for any m ≥ 1, there exists a constant K = K m > 0 (independent of p and y) such that Let l > 0 be an arbitrarily large positive number. Consider the following stochastic differential equation: Since coefficients of SDE (2.3) satisfy condition (4.1), and coefficients of SDE (4.2) satisfy the global Lipschitz condition, so there exists a unique solution y l (t), which converges to the exact solution y(t) of SDE (2.3) in probability as l goes to infinity. If we apply the Euler-Maruyama method to (4.2), we obtain a numerical approximation y ∆,l (t) to (4.2) which will strongly converge to y l (t) in L p (∀p > 0) sense. Therefore, for each ∆, there is a l ∆ such that y ∆,l ∆ (t) will converge to y(t) in probability. We will prove that under Assumption 4.1, both the exact solution and the numerical solution are exponentially integrable and like in [24,26] for each ∆, we can find a suitable l ∆ that makes the y ∆,l ∆ (t) strongly converge to y(t). Now let us describe the construction of the truncated Euler-Maruyama method. First, we assume that there exists a function ψ : [1, ∞) → (0, ∞), which is strictly increasing and lim R→+∞ ψ(R) = +∞, such that We also assume that the inverse ψ −1 : [ψ(1), ∞) → (0, ∞) of ψ is also increasing. At the same time, we define another strictly decreasing function h : The above equations together with (4.3) imply that Then, applying the Euler-Maruyama method to (4.4) yields Here, we set y ∆ (0) = y 0 . Define the continuous form of (4.6) as Next, we will prove the exponential integrability of both the exact solution y(t) of SDE (2.3) and the numerical solution y ∆ (t). From now on, we always assume that C stands generic constants, whose values might change every time as it appears. Proof. An application of the Itô formula to y p (t) yields Here we only consider the case where p is a positive integer.
Using the Gronwall inequality, we obtain Besides, using the definition of τ R , we can infer As a result, it can be derived that According to the Jensen inequality, we have It is a known fact that for any q > 0, we have Finally, if m < 2, by means of the Stirling formula, we can derive This completes the proof of the theorem.
Remark 4.4. There are some papers concerned exponential integrability of exact solution. For example, different kinds of exponential integrability are discussed in [27,28] for SDE whose diffusion coefficient is a constant or grows linearly at most. The exponential integrability discussed in this paper is applicable to SDEs with exponentially growing drift and diffusion coefficients.
Next, we prove that the truncated Euler-Maruyama solutions are also exponentially integrable. But first, we need to state some lemmas.
Proof. This can be easily proved by combining (4.5) with the above lemma. (4.12) In the above inequalities, we have used the basic inequality (a + b) p ≤ 2 p−1 (a p + b p ), ∀a, b ≥ 0, the inequality (4.5), and the Burkholder-Davis-Gundy inequality of Itô stochastic integral.
Proof. By the Itô formula and property of stochastic Itô integral, we have By means of the Young inequality ab ≤ a p p + b q q , where a, b > 0, p > 1, 1 By using the Young inequality again, we have Using the Gronwall inequality we have E y ∆ (s) 2p ≤ (y 2p 0 + C p tp mp )e C pt .
Hence, we can get that for arbitrary t > 0, E y ∆ (t) 2p ≤ C p t,y 0 (p) mp . The Jensen inequality yields E | y ∆ (t) | p ≤ C p t,y 0 (p) mp 2 , m < 2. Now, we can get the exponential integrability of truncated Euler-Maruyama approximation for any ∆ ∈ (0, 1] by the Stirling formula. Defineτ R = inf{t > 0 :| y ∆ (t) |≥ R} . (4.14) Since we already have the conclusion that both the exact solution and the numerical approximation are exponentially integrable, we can get Then the truncated Euler-Maruyama approximation y ∆ (t) will strongly converge to the exact solution y(t) of SDE (2.3), i.e. lim More importantly, the numerical solution x ∆ (t) maintains the positiveness of the original solution.
Proof. By means of the mean value theorem, we obtain The convergence result can be derived by using exponential integrability results and the above lemma.
Moreover, for the CIR process (1.2), after using logarithmic transformation, we obtain SDE (2.5). The coefficients of (2.5) does not satisfy the condition (4.1). This can be easily seen from the following computations. Denote When y > 0, it follows that A 0 ≤ K p. But when y < 0, denote The derivative has the following form If κλ < p 2 θ 2 , then we have . It is easy to see that H (y 1 ) < 0. We assume H (y 0 ) = 0. Then there exists a y 2 ∈ (y 0 , y 1 ), such that where c 2 = θ 2 2κλ−θ 2 and c 1 = 2κλ 2κλ−θ 2 . To cover the case of CIR and other more general SDE, we introduce the following assumption.
where the positive constant K = K p depends only on p, but independent of y. And there exists a positive constant K 1 = K 1,p which also depends only on p, such that e 2y [ f (y) + p 2 g 2 (y)] ≤ K 1 (1 + e 2y ), ∀y ∈ R , p > 0. (4.23) From [24], we know that under condition (4.22), the truncated Euler-Maruyama numerical solution y ∆ (t) strongly converges to the solution of SDE (2.3) and P(τ ∞ > T ) = 0. Since it is easy to see that the condition (4.23) becomes which implies E x(t) p ) < ∞, ∀p > 0, and consequently, Ee py(t) < ∞. In order to get strong convergence result, we need to prove Ee py ∆ (t) < ∞, and this is what the next theorem says. Ee py ∆ (t) < ∞.
Proof. By applying Itô formula, we derive For the last item on the right side of the above inequality, we have for t ∈ [t k , t k+1 ). Thus, according to the Gronwall inequality, we derive Ee py ∆ (t∧τ R ) ≤ (e py 0 + 2K 1 te C p 2 )e C pt .
Since P(τ ∞ > t) = 0 for any t > 0, an application of the Fatou lemma yields Proof. This can be easily proved by using similar techniques to the proof of Theorem 4.12.
Remark 4.19. It is not difficult to see that the coefficients of SDE (2.5) satisfy Assumption 4.16, therefore, Theorem 4.18 can be applied to CIR process (1.2).

Convergence rate of numerical scheme
Although we have obtained the strong convergence theorem for the logTE scheme, we still do not know how fast the numerical solution converges. In this section we impose stronger conditions on the coefficients than those in the previous section to get the convergence rate.
Assumption 5.1. Let the coefficients of SDE (2.3) satisfy the following conditions: there exist positive constants K = K m , K 1 = K 1,q and K 2 , where K 1 only depends on q, K 2 is independent of y 1 , y 2 , and K only depends on m, such that, for any y, y 1 , y 2 ∈ R, we have 2) The main theorem of this section is the following convergence rate of the numerical scheme constructed in the previous section.

Proof.
Denote An application of the Itô formula to e ∆ (t) 2q 1 yields where q 1 is positive integer satisfying q 1 < q. It follows from Assumption 5.1 that Besides, we can derive that Using the Young inequality twice, we have Denote π(y) = | y | ∧ψ −1 (h(∆)) y |y| . Obviously, we have | π(y) |≤| y |. We have On the other hand, using Lemma 4.8 and then the Young inequality, we get

Applications of the Gronwall inequality and the Jensen inequality yields
).
Taking p sufficiently large such that p( 1 −1)( −1) Finally, using the mean value theorem and the fact that both the exact solution and the numerical solution are exponentially integrable, we derive Although in the beginning of the proof, we assume q 1 > 2, other cases (0 < q 1 ≤ 2 and non integer case) can be easily dealt with by the Jensen inequality. The theorem is then proved.
In this case, since the diffusion coefficient is a constant, the fifth row of inequality (4.12) has a simple form: The above conclusions follow from similar argument in the proof of Theorem 5.2.
Example 5.4. Lwt return to Example 4.13. Since the coefficients of (4.19) satisfy and α ∈ ( 3 4 , 1), there exists a positive constant K, such that Then for ∀y 1 , y 2 ∈ R, we have Thus Theorem 5.2 can be applied to Example 4.13.
Example 5.5. Similarly, we can prove easily that if r > 4ρ − 3 the coefficients of (4.21) in Example 4.14 satisfy for all q > 2, and for all y 1 , y 2 ∈ R. This means that Theorem 5.2 can be applied to Example 4.14.
Example 5.6. Consider the scalar stochastic Ginzburg-Landau SDE The above SDE can be considered as a special case of the Ait-Sahalia process with (a −1 , a 0 , a 1 , a 2 , r, ρ) = (0, 0, λ + σ 2 2 , 3, 1), and satisfies r > 2ρ − 1. As mentioned above, this SDE has positive solution, and the boundaries are unattainable. We already know that the analytic solution of this SDE can be expressed as from which we can also see that the solution is actually positive, if x 0 > 0. After using the logarithmic transformation y(t) = ln x(t), we get dy(t) = [−e 2y(t) + λ]dt + σdW(t).
Apparently, the coefficients of the above SDE satisfy Assumption 5.1, and then Theorem 5.2 can be applied to this process.
The examples mentioned above take values in (0, +∞), so the new processes obtained by the logarithmic transformation take value in the whole line. But there are processes whose solutions take value in a certain domain (0, D), and in this case we can change the logarithmic transformation to an appropriate transformation accordingly.
Example 5.7. Consider the following stochastic susceptible-infected-susceptible (SIS) model (appeared in [30]) it is not difficult to get that I(t) + S (t) = M. Then we only need to study the following SDE with initial value I(0) = I 0 ∈ (0, M). The coefficients of SDE (5.8) satisfy the local Lipschitz condition, then we know this SDE will have unique maximal local solution on t ∈ [0, τ ε ), here the τ ε is the time to zero or M for the first time.
Denote ς = βM − µ − γ, the scale function has the following form where we used the known results where M 0 = 5000, x i (T ) and x i ∆ (T ) denote the number of sample paths, the ith exact solution and the ith numerical solution respectively. We show mean square errors of logTE method and TE method in (b) of Figure 1. Since the diffusion coefficient of the transformed SDE (5.6) is a constant, then according to Remark 5.3 we know that the convergence rate of logTE method is 1 in this case. Since the convergence order of TE method is close to 0.5, that is to say, for the GLE model, our method is superior to the TE method in terms of convergence rate. The result showed in (b) of Figure 1 is consistent with our theoretical results  Figure 2, from which we can see that approximation from our method can maintain the domain of the exact solution while the Euler-Maruyama approximation cannot. Since the explicit solution of this SDE can not be derived, the numerical solution with step size ∆ = 2 −14 is used as the reference solution. We compute five different numerical solutions of the logTE with step size ∆ = 2 −6 , , 2 −7 , · · · , 2 −10 , and show the estimated mean square convergence rate in (b) of Figure 2. From Figure 2, we can also find that the numerical results are consistent with our theoretical results given above. Besides, we also consider the CIR model (1.2). Set κ = 2, λ = 1, θ = 3 2 , x 0 = 2. Since the exact solution has no explicit expression, we use the drift-implicit square-root Euler method (SRE) proposed in [5] which can preserve the positiveness as a reference solution. We draw single trajectories of SRE with ∆ = 2 −14 and logTE with ∆ = 2 −6 in (a) of Figure 3. In addition, the coefficient of the SDE (2.5) does not satisfy Assumption 5.1, so we do not get the convergence rate of the logTE method for this model. We show the convergence order of this method in (b) of Figure 3 by six different numerical solutions with step size 2 −7 , 2 −8 , · · · , 2 −12 . For the last example we consider the CEV process (4.18). Set κ = 4, λ = 0.5, θ = 3, x 0 = 2, α = 7 8 . Similarly, since the exact solution of this model is difficult to get, we will use the approximation from semi-implicit method (denote as TSE) appeared in [2] as a reference solution. We compare single trajectories of the logTE method and the TSE method with same step size ∆ = 2 −6 in (a) of Figure 4. The numerical solution from TSE was proved to be able to maintain the positiveness of the original solution. From (a) of Figure 4, we can see that the trajectory of approximation obtained by logTE method is very close to that of TSE. The estimated mean square errors of logTE numerical solutions with six different step sizes ∆ = 2 −7 , 2 −8 , · · · , 2 −12 are shown in (b) of Figure 4, from which one clearly observes that mean square convergence rate is 1 2 for logTE.

Conclusion
In this paper, we propose a class of logarithmic Euler-Maruyama type methods for SDEs with positive solutions. Our idea is to first apply the logarithmic transformation y(t) = ln x(t), and then obtain a numerical approximation y ∆ (t) to the transformed process y(t), t ∈ [0, T ] by using Euler-Maruyama type method. And x ∆ (t) = e y ∆ (t) will be a positive approximation to the original solution x(t). We demonstrate some convergence and convergence rate conclusions for these methods in path-wise and L p sense. The using of logarithmic transformation makes the L p convergence analysis much more difficult, since the coefficients of the transformed process y(t) may grow exponentially and we need to deal with the exponential integrability of both y(t) and its approximated solution imposed by our transformation. We also note that for SDE whose solution takes value in a certain domain (0, D), approximate solution from our method is inside the domain of the original solution if the logarithmic transformation is changed accordingly, see Example 5.7.