Positivity-preserving truncated Euler and Milstein methods for financial SDEs with super-linear coefficients

The goal of this paper is to derive a positivity-preserving numerical method for scalar SDEs which takes values in $\mathbb{R}_{+}$ and have super-linear coefficients. Typical examples of such SDEs in mathematical finance and bio-mathematical applications are the Heston-3/2 volatility process [1]

with $c_{1},c_{2},\sigma>0$, and the Aït-Sahalia (AIT, for short) model [2]

where all constant parameters are nonnegative, $\kappa,\theta>1$, and $B(t)$ is a Wiener process. The SDEs appearing in such models are highly nonlinear and may contain the singularity in the neighbourhood of zero in the drift or diffusion coefficient. Such SDEs in almost all cases cannot be solved explicitly, and it has been and still is a very active topic of research to approximate SDEs with super-linear coefficients; see, e.g., [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], and the references therein.

Two significant challenges in constructing a numerical method for non-linear SDEs with positive solutions are to preserve positivity and to derive convergence with as high a rate as possible. In this section, we mainly discuss the positivity-preserving numerical methods for AIT (1.2). In fact, the techniques to handle the superlinearity and the singularity of the AIT coefficients can be applied to a more general SDE with positive solutions and super-linear coefficients. Table 1 gives a summary of the known results of some of the key methods which possess positivity preservation for AIT (1.2). We now discuss two main methods for the strong approximation of (1.2), the first is based on Lamperti/logarithmic transformation, and the second is direct approximation of AIT (1.2).

A classical transformation approach is to apply the Lamperti transformation $Y=X^{1-\theta}$ in order to obtain an auxiliary SDE in $Y$ with a globally Lipschitz diffusion, but a drift function which is unbounded when solutions are in a neighbourhood of zero. For Lamperti transformation methods that preserve positivity for (1.2) see, for example, Lamperti backward EM (BEM) [4], Lamperti projected EM [5], and Lamperti Semi-Discrete EM [21]. However, this type of polynomial transformation would shift all the nonsmoothness into the drift function. In recent years, much effort has focused on deriving the combination of logarithmic transformation and Euler/Milstein-type schemes for general SDEs which have positive solutions. For such SDEs, an effective explicit scheme that preserves positivity is the logarithmic truncated EM proposed by Yi et al. [11]. The method is to apply a logarithmic transformation to original SDE and to numerically approximate the transformed SDE by a truncated EM (an explicit method established in [25, 26]). Under the Khasminskii-type and the global monotonicity conditions on coefficients of transformed SDE, the method was shown to have strong $L^{p}$-order of convergence $1/2$. In a subsequent paper [16], a logarithmic truncated Milstein scheme was proved to have strong $L^{p}$-convergence of order $1$, see Table 1. However, when the logarithmic method in [11] is applied to AIT (1.2), the method requires a more restrictive condition: $\kappa>4\theta-3$ in order to satisfy the Khasminskii-type requirement for the transformed SDE, in contrast to the other positivity-preserving methods in Table 1. The same approach was also used by Tang and Mao in [22] to derive strong order $1/2$ for the SDEs with positive solutions, but weaker assumptions on the coefficients of original SDE was required.

The second main approach that preserve positivity of the numerical solution is to directly discretise (1.2) by using the implicit or truncated tricks. A backward EM (BEM) was shown to preserve positivity of the solution for AIT (1.2) in Szpruch et al. [27], where the authors proved the strong convergence of the method, but without revealing a rate. This gap was filled by Wang et al. in [23], where the authors successfully recovered the mean-square convergence rate of order $1/2$ for the theta EM applied to (1.2). However, the computational costs of solving implicit algebraic equations produced by classical BEM rise as the parameter $\kappa$ increases.

A truncated EM that seeks to directly approximate AIT was proposed by Emmanuel and Mao in [28], where the authors proved the strong $L^{p}$-convergence of the method, but they did not reveal a convergence rate. Recently, the $L^{p}$-convergence rate of order $\frac{1}{2p}$ of this method was obtained by Deng et al. [24]. At the same time, a novel semi-implicit plus tamed EM was introduced by Liu et al. in [19], where the authors established the strong $L^{p}$-convergence of rate $1/2$ based on an implicit treatment of the singular term $a_{-1}x^{-1}$ to ensure the property of positiveness and a suitable taming of the super-linear terms $-a_{2}x^{\kappa}$ and $bx^{\theta}$. It should be pointed out that this method is explicit, as the proposed scheme can be explicitly expressed by finding a positive root of a quadratic equation. Therefore, compared with the implicit schemes in [3, 4, 23] the computational costs in [19] can be significantly reduced. Similar approach was also used by Jiang et al. in [18] to derive strong order $1$ of semi-implicit projected Milstein for AIT process. See Table 1 for a comparison.

The motivation of this paper is the following observation: in [24] the theoretical $L^{p}$-convergence rate of truncated EM for AIT is only $\frac{1}{2p}$, rather than the optimal $\frac{1}{2}$, and it should be further increased. As a consequence, this paper is concerned with the optimal strong convergence rate of the positivity-preserving truncated methods for general SDE (2.2) with solution and super-linear coefficients. We note that the error analysis in [24] is based on the higher moment estimates of the numerical solution and the interpolation of the discretization scheme to continuous time. According to the theory of Mao [25, 26], it is inevitable to estimate the probability of truncated EM solution escaping from the truncated interval $(1/R,R)$, i.e., $\mathbb{P}(\rho_{\Delta,R}\leq T)$. However, [24, Lemma 3.2] implies that the upper bound for this probability is of order $1/2$ at most, i.e.,

which prevents the truncated EM in [24] from achieving the optimal strong order $1/2$. In this work, we seek to bypass this issue.

Motivated by the approaches of [5, 29], we adopt the error analysis based on the higher moment and inverse moment estimates of the true solution and the discretised scheme rather than the continuous time extension of numerical solution. We first prove a key Lemma 2.4 which shows the global monotonicity for the truncated functions and implies the stochastic C-Stability of the truncated EM in the sense of [29, Definition 3.2]. By the finiteness of the $p_{0}$th moments and $p_{1}$th inverse moments of the true solution $X$, we then obtain the local truncations errors of $f(X)$ and $g(X)$ in $L^{p}$, which implies the stochastic B-Consistency of the truncated EM in the sense of [29, Definition 3.3], see Lemma 3.1. Combining Lemma 2.4 with Lemma 3.1 allows us to establish the optimal $L^{p}$-convergence of order $1/2$ of positivity-preserving truncated EM for SDE (2.2) in the similar proofs as in [30, Theorem 5.36]. According to the global error estimate (3.14) in Theorem 3.2, we show that the truncated numerical solution $Y_{\Delta}$ has some finite moments and inverse moments, see Proposition 3.4. By $L^{p}$-convergence Theorem 3.2 and moment boundedness Proposition 3.4, we apply truncated EM to directly approximate the 3/2 and the AIT models, see Applications 4.1 and 4.2. Based on the Lamperti transformation, we also apply this method to indirectly approximate the CIR process, see Application 4.3. Combining the Milstein scheme with the truncated strategy, we also prove the strong $L^{2}$-convergence order $1$ for the truncated Milstein, in the similar method of analyzing the convergence rate of the truncated EM applied to SDE (2.2).

The main contributions of this paper are as follows:

• 1.

  In contrast with the results of [24], our truncated EM for AIT obtains the optimal $L^{p}$-convergence of order $\frac{1}{2}$, which is far better than the rate $\frac{1}{2p}$ in [24].

• 2.

  Compared with the projected EM in [5], where the method requires that the diffusion coefficients of SDEs are globally Lipschitz continuous, our truncated methods can apply to a wide class of SDEs with positive solutions and super-linearly growing diffusion coefficients.

• 3.

  Our error analysis method does not relies on the high-order moments and inverse moments estimates, the continuous time extension of the numerical solution. Thus our truncated methods can be applied to approximate SDEs with super-linear coefficients directly and also with sub-linear coefficients indirectly, see Applications 4.

The structure of this paper is as follows. In Section 2 we give the form of the SDE with positive solution and super-linear coefficients, specify some regularity and integrability constraints placed upon the coefficients for the main strong convergence theorems. In Section 3 we set up the framework and prove the optimal strong convergence order as well as establish the moment boundedness results for the positivity-preserving truncated EM and Milstein schemes. In Section 4 we apply the proposed methods to some financial SDEs, such as the 3/2 and the AIT as well as the CIR models. In Section 5 we numerically compare convergence and efficiency of several commonly used methods. Finally, we have included in a small appendix a couple of proofs to make this paper self contained.

Let $(\Omega,{\mathcal{F}},\mathbb{P})$ be a complete probability space with a filtration $\{{\cal F}_{t}\}_{t\geq 0}$ satisfying the usual conditions (i.e., it is increasing and right continuous while $\cal{F}_{\textrm{0}}$ contains all $\mathbb{P}$-null sets). Let $\mathbb{R}=(-\infty,\infty)$ and $\mathbb{R}_{+}=(0,\infty)$. We write $|\cdot|$ and $\langle\cdot,\cdot\rangle$ respectively for the Euclidean norm and the inner produce on $\mathbb{R}$. Let $\{B(t)\}_{t\geq 0}$ be a one-dimensional Brownian motion with respect to the normal filtration $\{{\cal F}_{t}\}_{t\geq 0}$. Let $\mathbb{E}$ denote the expectation. Denote by $L^{p}(\Omega;\mathbb{R})$, for $p>0$, the set of random variables $Z$ such that $\|Z\|_{p}:=(\mathbb{E}[|Z|^{p}])^{1/p}<\infty$. We denote by $\mathbb{E}_{t}[X]:=\mathbb{E}[X|\mathcal{F}_{t}]$ the conditional expectation given the filtration $\mathcal{F}_{t}$. For two real numbers $a$ and $b$, $a\vee b:=\max(a,b)$ and $a\wedge b:=\min(a,b)$. W denote by $\mathcal{C}^{2}(\mathbb{R}_{+})$ the space of twice differential functions with continuous derivatives on $\mathbb{R}_{+}$. For any $x\in\mathbb{R}_{+}$, we denote $\lceil x\rceil$ the rounded up integer. For a set $G$, its indicator function is denoted by $\mathbf{1}_{G}$. In the following it will be convenient to introduce the abbreviation

Moreover, given $V(x)\in\mathcal{C}^{2}(\mathbb{R},\mathbb{R}_{+})$, we define a functional $\mathbb{L}V:\mathbb{R}\to\mathbb{R}$ by

In this work, we frequently apply the weighted Young inequality

Consider a one-dimensional SDE of the form

Throughout this paper, we assume that SDE (2.2) has a unique strong solution in $\mathbb{R}_{+}$. We now formulate the conditions on the drift and the diffusion coefficient functions.

In what follows, $C$ will be used to denote a positive constant depending on $K$, $T$, $\alpha$, $\beta$, $p_{0}$, $p_{1}$ and $X_{0}$, but whose value may be different from line to line. We denote it by $C_{p}$ if it depends on an extra parameter $p$.

Assumption 2.2 implies that the monotonicity condition is fulfilled on $\mathbb{R}_{+}$, i.e.,

and

see [22, Remark 2.1]. Assumption 2.3 imposes a condition on the moments of the true solution $X$ in terms of the locally Lipschitz exponents $\alpha$, $\beta$. Combining this assumption and Assumption 2.1 allows us to bound the local truncations errors for $f(X)$ and $g(X)$, see Lemma 3.1.

Let $\Delta\in(0,1]$, define a truncation mapping $\pi_{\Delta}:\mathbb{R}\to\left[\frac{1}{R(\Delta)},R(\Delta)\right]$ by

where $R(\Delta):=L_{1}\Delta^{-\gamma}$ is the size of truncated interval with $L_{1}\in\left[{\frac{1}{X_{0}}}\vee X_{0},\infty\right)$ and $\gamma\in\left(0,\frac{1}{2(\alpha\vee\beta)}\right]$ is a constant to be determined later. Define the truncated functions

The following lemma shows that the truncated mapping $\pi_{\Delta}$ is $1$-Lipschitz continuous on $\mathbb{R}$, and the truncated functions $f_{\Delta}$ and $g_{\Delta}$ preserve the monotonicity condition (2.6).

The proof of this lemma is postponed to Appendix 6.1. From Lemma 2.4, we have,

Moreover, by Assumption 2.1 and (2.8) as well as (2.10), we have

where $\varphi(\Delta)=C\Delta^{-\gamma(\alpha\vee\beta)}$ with

We now provide two lemmas, their proofs are postponed to Appendix 6.2 and 6.3. Lemma 2.5 reveals the Hölder continuity of the true solution to (2.2) with respect to the norm in $L^{p}(\Omega;\mathbb{R})$, while Lemma 2.6 shows the $L^{p}(\Omega;\mathbb{R})$-error due to truncating the true solution $X$ on $[1/R,R]$.