On the weak convergence rate of an exponential Euler scheme for SDEs governed by coefficients with superlinear growth

TL;DR

This paper analyzes the weak convergence rate of an exponential Euler scheme for one-dimensional SDEs with superlinear coefficients, establishing a convergence rate of order one under certain conditions and supporting findings with numerical experiments.

Contribution

It introduces a semi-explicit exponential Euler scheme for SDEs with superlinear growth and proves its weak convergence rate of order one, extending existing methods.

Findings

Weak convergence rate of order one established

Scheme effective for SDEs with superlinear coefficients

Numerical experiments support theoretical results

Abstract

We consider the problem of the approximation of the solution of a one-dimensional SDE with non-globally Lipschitz drift and diffusion coefficients behaving as $x^\alpha$, with $\alpha>1$. We propose an (semi-explicit) exponential-Euler scheme and study its convergence through its weak approximation error. To this aim, we analyze the $C^{1,4}$ regularity of the solution of the associated backward Kolmogorov PDE using its Feynman-Kac representation and the flow derivative of the involved processes. From this, under some suitable hypotheses on the parameters of the model ensuring the control of its positive moments, we recover a rate of weak convergence of order one for the proposed exponential Euler scheme. Finally, numerical experiments are shown in order to support and complement our theoretical result.

Tables

Table 1. Table 1: Observed numerical weak error | 𝔼 ​ [ f ​ ( X T ) ] − 1 n ​ ∑ i = 1 n f ​ ( X ¯ T ​ ( ω i , 2 − p ) ) | 𝔼 delimited-[] 𝑓 subscript 𝑋 𝑇 1 𝑛 superscript subscript 𝑖 1 𝑛 𝑓 subscript ¯ 𝑋 𝑇 subscript 𝜔 𝑖 superscript 2 𝑝 |\mathbb{E}[f(X_{T})]-\frac{1}{n}\sum_{i=1}^{n}f(\overline{X}_{T}(\omega_{i},2^{-p}))| .

Weak Error with $\Delta t=2^{-p}$, for $p=2,\dots ,9$
Case 1: $(B_2,\sigma ,\alpha )=(2,\frac{1}{10},\frac{3}{2})$ and H5 is valid.
Test function	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	3.397e-2	1.606e-2	7.756e-3	3.823e-3	1.923e-3	1.033e-3	4.965e-4	3.199e-4
$f\left(x\right)=x^2$	2.147e-2	1.043e-2	5.102e-3	2.529e-3	1.277e-3	6.864e-4	3.297e-4	2.131e-4
$f\left(x\right)=e^{-x^2}$	1.94e-2	9.378e-3	4.568e-3	2.258e-3	1.135e-3	6.06e-4	2.874e-4	1.829e-4
Case 2: $(B_2,\sigma ,\alpha )=(3,1,\frac{5}{4})$ and H5 is not valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	2.179e-2	1.069e-2	5.07e-3	2.529e-3	1.32e-3	8.021e-4	2.598e-4	3.043e-4
$f\left(x\right)=x^2$	6.412e-3	3.243e-3	1.582e-3	8.397e-4	3.965e-4	2.148e-4	9.101e-5	5.065e-5
$f\left(x\right)=e^{-x^2}$	6.113e-3	3.07e-3	1.5e-3	7.5e-4	3.65e-4	1.868e-4	5.625e-5	5.439e-5
Case 3: $(B_2,\sigma ,\alpha )=(1,1,\frac{3}{2})$ and H5 is not valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	2.3e-2	1.219e-2	5.864e-3	2.893e-3	1.255e-3	9.507e-4	3.13e-4	3.14e-4
$f\left(x\right)=x^2$	9.408e-3	5.302e-5	1.749e-4	2.956e-4	8.41e-3	5.95e-4	2.574e-3	2.52e-4
$f\left(x\right)=e^{-x^2}$	1.485e-2	8.108e-3	4.162e-3	2.248e-3	1.31e-3	1.164e-3	3.78e-4	2.762e-4
Case 4: $(B_0,B_1,B_2,\sigma ,\alpha )=(1,1,\frac{2}{5},\frac{1}{10},3)$ and H5 is not valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	3.741e-3	3.292e-3	2.103e-3	1.364e-3	6.915e-4	2.936e-4	1.131e-4	5.497e-5
$f\left(x\right)=x^2$	2.687e-2	1.332e-2	7.476e-3	4.584e-3	2.302e-3	1.034e-3	4.423e-4	2.312e-4
$f\left(x\right)=e^{-x^2}$	5.027e-3	2.846e-4	2.372e-4	2.666e-4	1.545e-4	5.72e-5	2.363e-5	1.529e-5

Table 2. Table 2: Observed numerical weak error | 𝔼 ​ [ f ​ ( X T ) ] − 1 n ​ ∑ i = 1 n f ​ ( X ¯ T ​ ( ω i , 2 − p ) ) | 𝔼 delimited-[] 𝑓 subscript 𝑋 𝑇 1 𝑛 superscript subscript 𝑖 1 𝑛 𝑓 subscript ¯ 𝑋 𝑇 subscript 𝜔 𝑖 superscript 2 𝑝 |\mathbb{E}[f(X_{T})]-\frac{1}{n}\sum_{i=1}^{n}f(\overline{X}_{T}(\omega_{i},2^{-p}))| for Case 1 to Case 7 .

Weak Error with $\Delta t=2^{-p}$, for $p=2,\dots ,9$
Case 1: $(B_0,B_1,B_2,\sigma ,\alpha )=(0,0,2,\frac{1}{10},\frac{3}{2})$ and H5 is valid.
Test function	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	3.397e-2	1.606e-2	7.756e-3	3.823e-3	1.923e-3	1.033e-3	4.965e-4	3.199e-4
$f\left(x\right)=x^2$	2.147e-2	1.043e-2	5.102e-3	2.529e-3	1.277e-3	6.864e-4	3.297e-4	2.131e-4
$f\left(x\right)=\frac{1}{x}$	3.451e-1	1.521e-1	7.311e-2	3.556e-2	1.678e-2	9.122e-3	4.793e-3	2.428e-3
$f\left(x\right)=e^{-x^2}$	1.94e-2	9.378e-3	4.568e-3	2.258e-3	1.135e-3	6.06e-4	2.874e-4	1.829e-4
Case 2: $(B_0,B1,B_2,\sigma ,\alpha )=(0,0,3,1,\frac{5}{4})$ and H5 is not valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	2.179e-2	1.069e-2	5.07e-3	2.529e-3	1.32e-3	8.021e-4	2.598e-4	3.043e-4
$f\left(x\right)=x^2$	6.412e-3	3.243e-3	1.582e-3	8.397e-4	3.965e-4	2.148e-4	9.101e-5	5.065e-5
$f\left(x\right)=\frac{1}{x}$	2.208	9.48e-1	4.503e-1	2.215e-1	1.1e-1	5.93e-2	3.68e-2	1.526e-2
$f\left(x\right)=e^{-x^2}$	6.113e-3	3.07e-3	1.5e-3	7.5e-4	3.65e-4	1.868e-4	5.625e-5	5.439e-5
Case 3: $(B_0,B_1,B_2,\sigma ,\alpha )=(0,0,1,1,\frac{3}{2})$ and H5 is not valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	2.3e-2	1.219e-2	5.864e-3	2.893e-3	1.255e-3	9.507e-4	3.13e-4	3.14e-4
$f\left(x\right)=x^2$	9.408e-3	5.302e-5	1.749e-4	2.956e-4	8.41e-3	5.95e-4	2.574e-3	2.52e-4
$f\left(x\right)=\frac{1}{x}$	1.084	1.058	8.919e-2	4.341e-2	2.483e-2	1.224e-2	9.783e-3	4.861e-3
$f\left(x\right)=e^{-x^2}$	1.485e-2	8.108e-3	4.162e-3	2.248e-3	1.31e-3	1.164e-3	3.78e-4	2.762e-4
Case 4: $(B_0,B_1,B_2,\sigma ,\alpha )=(1,1,\frac{2}{5},\frac{1}{10},3)$ and H5 is not valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	3.741e-3	3.292e-3	2.103e-3	1.364e-3	6.915e-4	2.936e-4	1.131e-4	5.497e-5
$f\left(x\right)=x^2$	2.687e-2	1.332e-2	7.476e-3	4.584e-3	2.302e-3	1.034e-3	4.423e-4	2.312e-4
$f\left(x\right)=\frac{1}{x}$	4.374e-3	2.804e-3	5.505e-4	4.508e-4	2.383e-4	8.651e-5	2.313e-5	9.249e-6
$f\left(x\right)=e^{-x^2}$	5.027e-3	2.846e-4	2.372e-4	2.666e-4	1.545e-4	5.72e-5	2.363e-5	1.529e-5
Case 5: $(B_0,B_1,B_2,\sigma ,\alpha )=(0,0,10,\frac{1}{2},\frac{9}{8})$ and H5 is valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	2.735e-3	1.413e-3	7.122e-4	3.69e-4	1.873e-4	9.082e-5	5.053e-5	1.979e-5
$f\left(x\right)=x^2$	2.955e-5	1.718e-5	9.153e-6	4.892e-6	2.511e-6	1.226e-6	6.833e-7	2.82e-7
$f\left(x\right)=\frac{1}{x}$	1.175e+2	4.523e+1	2.025e+1	9.639	4.765	2.35	1.202	5.874e-1
$f\left(x\right)=e^{-x^2}$	2.952e-5	1.715e-5	9.128e-6	4.867e-6	2.486e-6	1.201e-6	6.583e-7	2.571e-7
Case 6: $(B_0,B_1,B_2,\sigma ,\alpha )=(0,0,\frac{1}{100},\frac{1}{10},\frac{5}{4})$ and H5 is not valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	8.51e-5	1.756e-5	6.374e-5	6.616e-5	1.943e-4	6.087e-5	1.69e-5	2.31e-5
$f\left(x\right)=x^2$	1.538e-4	4.004e-5	1.187e-4	1.334e-4	4.363e-4	2.421e-4	3.981e-5	4.204e-5
$f\left(x\right)=\frac{1}{x}$	5.64e-5	2.389e-4	4.294e-5	2.704e-4	1.394e-4	2.57e-4	1.724e-4	7.968e-5
$f\left(x\right)=e^{-x^2}$	1.623e-5	1.033e-4	1.396e-4	1.372e-4	3.481e-5	1.251e-4	7.891e-5	7.133e-5
Case 7: $(B_0,B_1,B_2,\sigma ,\alpha )=(0,0,\frac{2}{5},\frac{1}{10},3)$ and H5 is valid.
	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$	$p=9$
$f\left(x\right)=x$	1.698e-2	7.814e-3	3.881e-3	1.844e-3	9.058e-4	4.943e-4	2.266e-4	8.054e-5
$f\left(x\right)=x^2$	2.608e-2	1.207e-2	6.018e-3	2.864e-3	1.405e-3	7.688e-4	3.515e-4	1.489e-4
$f\left(x\right)=\frac{1}{x}$	2.902e-2	1.323e-2	6.337e-3	3.156e-3	1.387e-3	9.119e-4	3.679e-4	2.059e-4
$f\left(x\right)=e^{-x^2}$	1.436e-2	6.571e-3	3.232e-3	1.505e-3	7.1e-4	3.609e-4	1.341e-4	1.122e-5

Table 3. Table 3: Comparison of the weak approximation error for test function f ​ ( x ) = x 𝑓 𝑥 𝑥 f(x)=x . The comparison consider the following numerical schemes: exponential Euler, Symmetrized Euler and Milstein schemes, Tamed and Stopped Tamed Euler schemes.

	Observed weak Error with $\Delta t=2^{-p}$
Cases $(B_2,\sigma ,\alpha$)	$p=1$	$p=2$	$p=3$	$p=4$	$p=5$	$p=6$	$p=7$	$p=8$
Case 1 $(2,\frac{1}{10},\frac{3}{2})$
exp-ES	7.866e-2	3.402e-2	1.6e-2	7.829e-3	3.939e-3	1.918e-3	8.92e-4	4.774e-4
SES	1.557e-4	2.866e-4	2.573e-4	1.328e-4	1.747e-4	4.402e-4	4.66e-4	1.602e-4
SMS	9.462e-2	2.272e-2	3.028e-2	1.065e+17	1.480e-2	5.901e-3	3.103e-3	9.497e-4
STES	6.939e-2	5.43e-2	2.969e-2	1.518e-2	7.906e-3	3.383e-3	2.111e-3	6.74e-4
TES	3.419e-2	2.12e-2	1.14e-2	6.091e-3	3.109e-3	1.759e-3	7.997e-4	3.626e-4
Case 2 $(3,1,\frac{5}{4})$
exp-ES	4.655e-2	2.237e-2	1.05e-2	5.263e-3	2.872e-3	1.489e-3	2.852e-4	2.458e-4
SES	3.815e-1	5.788e-2	4.519e-2	2.219e-2	1.048e-2	4.736e-3	2.434e-3	1.128e-3
SMS	2.813e-1	7.577e-2	3.312e-2	1.59e-2	7.82e-3	3.798e-3	1.86e-3	9.431e-4
STES	-	-	-	-	6.596e-3	3.411e-3	5.594e-4	1.895e-3
TES	-	-	-	-	4.292e-3	1.177e-3	8.831e-4	2.887e-4
Case 3 $(1,1,\frac{3}{2})$
exp-ES	1.195e-2	6.353e-3	2.94e-3	1.523e-3	8.059e-4	5.286e-4	3.573e-6	2.319e-5
SES	6.306e-3	6.306e-3	6.384e-3	6.571e-3	6.308e-3	6.416e-3	6.008e-3	5.995e-3
SMS	5.396e-2	1.567e-1	1.841e+6	1.411e-1	1.243e-1	1.178e-1	1.134e-1	1.116e-1
STES	-	-	-	-	-	5.979e-3	8.053e-3	1.758e-3
TES	-	-	-	-	-	-	-	2.715e-3
Case 5 $(10,\frac{1}{2},\frac{9}{8})$
exp-ES	4.980e-3	2.739e-3	1.416e-3	7.226e-4	3.707e-4	1.837e-4	8.583e-5	4.565e-5
SES	24.287	17.878	4.936e-3	3.262e-3	1.746e-3	9.092e-4	4.598e-4	2.655e-4
SMS	25.272	20.77	4.968e-3	3.556e-3	2.156e-3	1.45e-3	1.336e-3	1.654e-3
STES	5.49e-1	-	-	2.849e-3	1.652e-3	8.828e-4	4.905e-4	2.324e-4
TES	-	-	-	2.252e-3	1.21e-3	6.29e-4	3.17e-4	1.938e-4
Case 6 $(\frac{1}{100},\frac{1}{10},\frac{5}{4})$
exp-ES	4.728e-3	5.047e-3	4.656e-3	4.998e-3	5.406e-3	4.959e-3	4.587e-3	4.917e-3
SES	9.382e-1	7.76e-2	7.229e-2	5.053e-2	3.962e-2	3.5e-2	3.231e-2	3.109e-2
SMS	2.816e-1	7.682e-2	3.409e-2	1.674e-2	8.629e-3	4.738e-3	2.724e-3	1.783e-3
STES	8.355e-4	8.966e-4	6.832e-4	2.909e-4	1.913e-3	8.137e-4	2.142e-4	3.546e-4
TES	1.582e-4	1.222e-3	1.698e-3	3.282e-4	7.02e-5	8.031e-4	9.383e-4	9.324e-4

Equations

$$d{X}_t=b(X_t)dt+\sigma X_t^{\alpha }d{W}_t,X_0=x>0,$$

$$b(x)\le B_{1}x-B_2{x}^{2\alpha -1}+b(0),\forall x\in {\mathbb{R}}^{+},$$

$$d{X}_t=-B{X}_t^{2\alpha -1}dt+\sigma X_t^{\alpha }d{W}_t,X_0=x>0,$$

$$dX_{t}=X_{t}\big{(}\tfrac{b(X_{t})}{X_{t}}dt+\sigma{X_{t}^{\alpha-1}}dW_{t}\big{)},\quad X_{0}=x>0,$$

$$\overline{X}_{t_{n+1}}=\overline{X}_{t_{n}}\exp\{\sigma\overline{X}_{t_{n}}^{\alpha-1}\;(W_{t_{n+1}}-W_{t_{n}})+\big{(}\tfrac{b(\overline{X}_{t_{n}})-b^{+}(0)}{\overline{X}_{t_{n}}}-\frac{\sigma^{2}}{2}\overline{X}_{t_{n}}^{2(\alpha-1)}\big{)}\Delta t\}+b^{+}(0)\Delta t,$$

$$|f(x)-f(y)|\leq C\left(1+|x|^{\overline{\gamma}}+|y|^{\overline{\gamma}}+|x|^{-\underline{\gamma}}+|y|^{-\underline{\gamma}}\right)|x-y|,~{}\forall x,y\in\textsl{dom$(f)$-\{$0$\}}.$$

$$|f(x)-f(y)|\le C\left(1+|x{|}^{\gamma }+|y{|}^{\gamma }\right)|x-y|,\forall x,y\in \text{textsl}dom(f).$$

$$b(x)\le B_{1}x-B_2{x}^{2\alpha -1}+b(0).$$

$$\displaystyle\sup_{t\in[0,T]}\mathbb{E}\big{[}X_{t}^{2p}\big{]}\leq~{}C_{p}(1+x^{2p}),$$

$$\underset{t\in [0,T]}{\sup }\mathbb{E}[X_t^{-q}]\le C_q(1+x^{-q}),$$

$$\sup_{t\in[0,T]}\mathbb{E}\big{[}\exp\{\mu\int_{0}^{t}X_{s}^{2\alpha-2}ds\}\big{]}\leq\;C\;\left(1+x^{\frac{1}{2}+\frac{B_{2}}{\sigma^{2}}}\right).$$

$$\sup_{t\in[0,T]}\mathbb{E}\big{[}\exp\{\mu\int_{0}^{t}X_{s}^{2\alpha-2}ds\}\big{]}\leq\;C(1+x^{\frac{2\mu}{\sigma^{2}}})~{}\left(1+\exp\{C\mu x^{-1}\}\right).$$

$$dX_{t}=X_{t}\Big{(}\tfrac{b(X_{t})}{X_{t}}dt+\sigma X^{\alpha-1}dW_{t},\Big{)},~{}~{}X_{0}=x>0,$$

$$\widehat{X}_{t_{n+1}}=\widehat{X}_{t_{n}}\exp\Big{\{}\Big{(}\tfrac{b(\widehat{X}_{t_{n}})}{\widehat{X}_{t_{n}}}-\frac{\sigma^{2}}{2}\widehat{X}_{t_{n}}^{2(\alpha-1)}\Big{)}\Delta t+\sigma\overline{X}_{t_{n}}^{\alpha-1}\;(W_{t_{n+1}}-W_{t_{n}})\Big{\}},~{}~{}\widehat{X}_{0}=x,$$

$$d\widehat{X}_{t}=\widehat{X}_{t}\Big{(}\tfrac{b(\widehat{X}_{\eta(t)})}{\widehat{X}_{\eta(t)}}dt+\sigma\widehat{X}^{\alpha-1}_{\eta(t)}dW_{t}\Big{)},~{}~{}X_{0}=x>0,$$

$$\frac{b(x)}{x}\le B_1-B_2{x}^{2(\alpha -1)},\forall x\ge 0,$$

$$d\widehat{X}_{t}=b(0)dt+\widehat{X}_{t}\Big{(}\tfrac{b(\widehat{X}_{\eta(t)})-b(0)}{\widehat{X}_{\eta(t)}}dt+\sigma\widehat{X}^{\alpha-1}_{\eta(t)}dW_{t}\Big{)},$$

$$\displaystyle\widehat{X}_{t}=\widehat{X}_{\eta(t)}\exp\Big{\{}\sigma\widehat{X}_{\eta(t)}^{\alpha-1}\;(W_{t}-W_{\eta(t)})+\Big{(}\tfrac{b(\widehat{X}_{\eta(t)})-b(0)}{\overline{X}_{\eta(t)}}-\tfrac{\sigma^{2}}{2}\widehat{X}_{\eta(t)}^{2(\alpha-1)}\Big{)}{\delta(t)}+\int_{\eta(t)}^{t}\tfrac{b(0)}{\widehat{X}_{s}}ds\Big{\}},$$

$${\int }_{\eta (t)}^t\frac{b(0)}{{\hat{X}}_s}ds\approx \frac{b(0)}{{\hat{X}}_{\eta (t)}}\delta (t),$$

$${\int }_{\eta (t)}^t\frac{b(0)}{{\hat{X}}_s}ds\approx \frac{b(0)}{{\hat{X}}_t}\delta (t),$$

$$f\left(t,\check{X}_{t}\right)=\check{X}_{\eta(t)}\exp\big{\{}\sigma\check{X}_{\eta(t)}^{\alpha-1}\;(W_{t}-W_{\eta(t)})+\big{(}\tfrac{b\left(\check{X}_{\eta(t)}\right)-b(0)}{\check{X}_{\eta(t)}}-\tfrac{\sigma^{2}}{2}\check{X}_{\eta(t)}^{2(\alpha-1)}\big{)}{\delta(t)}\big{\}},$$

$$f(t,x)\approx x-b(0)\delta (t).$$

$$\overline{X}_{t_{n+1}}=b(0)\Delta t+\overline{X}_{t_{n}}\exp\Big{\{}\sigma\overline{X}_{t_{n}}^{\alpha-1}(W_{t_{n+1}}-W_{t_{n}})+\big{(}\tfrac{b\left(\overline{X}_{t_{n}}\right)-b(0)}{\overline{X}_{t_{n}}}-\tfrac{\sigma^{2}}{2}\overline{X}_{t_{n}}^{2(\alpha-1)}\big{)}\Delta t\Big{\}},$$

$$\overline{X}_{t}=b(0){\delta(t)}+\overline{X}_{\eta(t)}\exp\Big{\{}\sigma\overline{X}_{\eta(t)}^{\alpha-1}(W_{t}-W_{\eta(t)})+\big{(}\tfrac{b\left(\overline{X}_{\eta(t)}\right)-b(0)}{\overline{X}_{\eta(t)}}-\tfrac{\sigma^{2}}{2}\overline{X}_{\eta(t)}^{2(\alpha-1)}\big{)}{\delta(t)}\Big{\}},$$

$$d\overline{X}_{t}=\left(\overline{X}_{t}-b(0){\delta(t)}\right)\Big{(}\tfrac{b\left(\overline{X}_{\eta(t)}\right)-b(0)}{\overline{X}_{\eta(t)}}dt+\sigma\overline{X}_{\eta(t)}^{\alpha-1}dW_{t}\Big{)}+b(0)dt.$$

$$\sup_{t\in[0,T]}\mathbb{E}\big{[}{\overline{X}_{t}^{2p}}\big{]}\leq C_{p}(1+x^{2p}),\,x>0.$$

$$\displaystyle\mathbb{E}[\overline{X}_{t}^{2p}]=x^{2p}+2pb(0)\mathbb{E}\big{[}\int_{0}^{t\wedge\tau_{M}}\overline{X}_{s}^{2p-1}ds\big{]}$$

$$\displaystyle+p\mathbb{E}\Big{[}\int_{0}^{t}\overline{X}_{s}^{2p-2}(\overline{X}_{s}-b(0){\delta(s)})\Big{\{}2\overline{X}_{s}\tfrac{b(\overline{X}_{\eta(s)})-b(0)}{\overline{X}_{\eta(s)}}+(2p-1)\sigma^{2}\left(\overline{X}_{s}-b(0){\delta(s)}\right)\overline{X}_{\eta(s)}^{2(\alpha-1)}\Big{\}}ds\Big{]}$$

$$\displaystyle\leq x^{2p}+2pb(0)\int_{0}^{t}\mathbb{E}[\overline{X}_{s}^{2p-1}]ds+2pB_{1}\mathbb{E}\big{[}\int_{0}^{t}\overline{X}_{s}^{2p-1}(\overline{X}_{s}-b(0){\delta(s)})ds\big{]}$$

$$\le C(1+x^{2p})+C{\int }_0^t\mathbb{E}[{\overline{X}}_s^{2p}]ds.$$

$$\displaystyle|b^{(i)}(x)|\leq C\big{(}1+x^{\overline{\gamma}_{(i)}+1}+x^{-\underline{\gamma}_{(i)}}\big{)},~{}\mbox{ for }i=1,2,3,4,$$

$$b(x)\le B_{1}x-B_2{x}^{2\alpha -1}+b(0),\text{mbox}and{b}^{\prime }(x)\le B_1^{\prime }-B_2^{\prime }{x}^{2(\alpha -1)}.$$

Full text

Keywords.

*Stochastic differential equation, Numerical scheme, Polynomial coefficients, Weak convergence, Rate of convergence; *

MSC 2010 subject classifications. 65C20. 60-08. 91G60.

1 Introduction

Within the extensive literature on the numerical analysis of approximation schemes for Brownian-driven stochastic differential equations with non-Lipschitz coefficients, existing convergence results mainly deal separately with the singularity hypothesis on the drift coefficient or on the diffusion coefficient. More rarely, the Lipschitz property is dropped for both coefficients. In this paper, we propose a numerical scheme for one-dimensional stochastic differential equations (SDEs for short) having non-globally Lipschitz, polynomial drift and diffusion coefficients, and we analyze its convergence for the weak error. In this context, we present the first direct proof of the weak convergence with rate one, accompanied by an extendable methodology to analyze the $C^{1,4}$ regularity of the Feynman-Kac representation involving the exact process.

More precisely, we are interested in the numerical approximation of the solution to the following class of SDEs

[TABLE]

where $(W_{t};0\leq t\leq T)$ is a standard Brownian motion on the probability space $(\Omega,\mathcal{F},\mathbb{P})$ equipped with its natural filtration $(\mathcal{F}_{t};\,0\leq t\leq T)$, and the power $\alpha$, characterizing the diffusion, is assumed strictly greater than one. The drift $b:[0,+\infty)\rightarrow\mathbb{R}$ is a locally Lipschitz continuous function allowing polynomial dependence in the Lipschitz constant (see Definition 1.1 below for a precise statement), and with a polynomial growth bound of the form:

[TABLE]

for some constants $B_{1},B_{2},b(0)\geq 0$. Since exponent $\alpha$ can be non-integer, some particular hypotheses under which the SDE (1.1) has a unique positive solution have to be made (see Proposition 2.1).

Convergence results in this particular setting of non-Lipschitz coefficients rarely deal directly on the weak error analysis. In this particular setting, Gyöngy [12] obtained pathwise almost surely convergence, with a convergence rate’s order of at most $\tfrac{1}{4}$, for the classical Euler-Maruyama scheme applied to SDEs with locally Lipschitz continuous drift and diffusion coefficients satisfying some Lyapunov condition. Such result immediately implies weak convergence for continuous bounded test functions but not $L^{p}$-strong convergence. Similarly, Higham, Mao and Stuart [14] established $L^{2}$-strong convergence of the classical Euler-Maruyama scheme for locally Lipschitz coefficients but assuming a priori the control of some $p$ th-moments ($p>2$) of the continuous solution of the SDE and of its approximation. A rate of strong convergence of order $\tfrac{1}{2}$ was also established for the time-implicit split-step backward Euler-Maruyama scheme, when the diffusion is globally Lipschitz, and the drift satisfies a one-sided Lipschitz condition and locally Lipschitz condition. Staying in the framework of the Euler-Maruyama scheme, Yan [27] obtained the weak convergence with diffusion and drift coefficients continuous only almost everywhere and having at most linear growth.

With superlinear growth coefficients, classical Euler-Maruyama scheme may present some degenerated behavior. Hutzenthaler, Jentzen and Kloeden [17] established the $L^{p}$-strong divergence, for $p\in[1,+\infty)$, related to the Euler-Maruyama scheme for SDEs with both drift and diffusion satisfying some superlinear growth condition. In particular, the authors obtained the divergence of the moments of the Euler approximation. Later in [18], the authors proposed a time-explicit tamed-Euler scheme to overcome this divergence problem of the Euler approximation, based on renormalized-increments to the scheme. Recently Hutzenthaler and Jentzen [16] proved the $\tfrac{1}{2}$ rate of the $L^{p}$-strong convergence for the tamed-Euler scheme for a family of SDE that includes some locally Lipschitz cases for both continuous drift and diffusion coefficients.

In the same vein of explicit in time alternative scheme to the $L^{p}$-strongly divergent Euler-Maruyama scheme, Sabanis [26] obtained the $L^{p}$-strong convergence for a scheme with renormalized coefficients under some superlinear growth condition, and recovered the $\tfrac{1}{2}$-$L^{p}$-strong convergence rate under Lispchitz diffusion and one-sided global Lispchitz drift.

Other time-explicit numerical schemes have been proposed over the years to solve the approximation problem of SDEs with locally Lipschitz continuous coefficients. For instance, Lamba, Mattingly and Stuart [21] proposed an adaptive Euler algorithm based on the control of the drift coefficient, and proved the $L^{2}$-strong convergence assuming the control of some moments of the solution and of its approximation. Chassagneux, Jacquier and Mihaylov [7] considered the case of globally Lipschitz diffusion and locally-Lipschitz drift function satisfying a one-sided Lipschitz condition and proposed a modified explicit Euler numerical scheme, for which, under suitable assumptions on the control of some moments of the solution, an $L^{2}$-strong convergence with explicit rate is proven.

Fewer results dealing with weak convergence are available. Milstein and Tretyakov [23] established a weak convergence result for a class of SDEs with non-globally Lipschitz coefficients, based on existing schemes with known rate of weak convergence for Lipschitz and smooth coefficients, and on the rejection of the approximated trajectories that go out a sphere of given radius. But the relation between the level of error, the radius of the rejection sphere and the time step threshold to be used in order to observe the convergence is not explicit, making the algorithm difficult to use in practice.

In this paper, we propose a new scheme for SDEs with smooth coefficients under some superlinear growth condition. The scheme is designed to ease the upper bound control of some moments of the approximated process and we prove the optimal convergence rate of order one for the weak error. The convergence analysis extends the methodology introduced in Bossy and Diop [4] to establish the regularity of the associated backward Kolmogorov PDE.

Our motivating problem

Our interest for the numerical approximation of (1.1) was initially motivated by the simulation/calibration problem for the instantaneous turbulent kinetic energy model issued from the Lagrangian description of a non inertial particle dynamics within a turbulent fluid flow (see [22, Chapter 1]). Such model can be described by a SDE having the prototype form:

[TABLE]

where $\alpha>1$. To your knowledge, no weak convergence rate for this model are available. Only strong convergence results are proposed. Equation (1.3) is a particular case of (1.1) and can be seen as a generalized Constant Elasticity of Variance (CEV) model (see e.g. Delbaen and Sirakawa [9]). In particular, the transformation $r_{t}=\tfrac{X_{t}^{2(1-\alpha)}}{4\sigma^{2}(\alpha-1)^{2}}$ applied to the solution of (1.3) produces the so-called CIR process (Cox, Ingersoll and Ross [8]) classically used for modeling short interest rate dynamics, and for which various schemes have been considered over the years. For the $L^{p}$-strong convergence of some proposed explicit schemes for CEV models, we refer to Bossy and Olivero [5] and the reference therein; for implicit proposed schemes, we refer to Dereich, Neuenkirch and Szpruch [10], Alfonsi [2] and the references therein.

Alternatively, the transformation $Y_{t}=\tfrac{X_{t}^{(1-\alpha)}}{\sigma(\alpha-1)}$ produces a Bessel process for which we can use an Explicit Projected Euler scheme proposed in [7], obtaining a strong rate of convergence of order $\frac{1}{6}$ provided that we control up to the $4(\alpha-1)$-th moments of the process $(X_{t};0\leq t\leq T)$ (or higher rate of convergence by controlling higher moments).

Exponential scheme

The keyword exponential scheme refers to generic semi-linear integration methods and is of main importance in numerical analysis. Methods for ODEs proposing integration schemes based on the semi-linear integration of equations are classics (see e.g. Pope [24], Hochbruck and Ostermann [15] and the references therein). Extend this methodology for SDEs is straightforward (particularly in dimension one where affine diffusions allow exact scheme), but establishing the weak rate of convergence results in the context of non-globally Lipschitz coefficients is much more demanding. With the same appellation, semi-linear integration methods are proposed for PDEs or SPDEs and concern schemes based on a mild formulation of the underlying equations (see e.g. Beccari et al [3] for SPDE problems with superlinear coefficients).

We would like to stress out that we were looking for a scheme that, potentially applied to prototype model (1.3), allows weak convergence rate of order one to set up an efficient calibration method for the model. Motivated by this problem for the model (1.3), the requirement of stability condition on the moments brings us to the variant scheme (1.4) below, as a remedy for the divergence problem of the Euler-Maruyama scheme, and an alternative to the tamed-Euler scheme (for which the weak convergence rate of order one is not established).

The proposed numerical approximation, which will be referred, from now on, to as the exponential-Euler scheme (exp-ES, for short), originates from rewriting the SDE (1.1) into

[TABLE]

and semi-linear integration produces, for an homogeneous $N$-partition of the time interval $[0,T]$ with time-step $\Delta t=t_{n+1}-t_{n}$, the approximation algorithm:

[TABLE]

(where $b^{+}(0)$ stands for $b(0)\vee 0$) that preserves the positiveness of the solution. We refer the reader to Section 3 for a detailed construction of (1.4).

The exponential Euler scheme (1.4) can be applied to large family of SDEs with non-globally Lipschitz coefficients, having strictly positive solution. The range of possible applications of our results includes some meaningful financial models such as the generalized CEV model, the non-linear mean reversion model (see e.g., Ait-Sahalia [1], Higham et al. [13]) and the Chan-Karolyi-Longstaff-Sanders model [6] among others.

As it will be established later on, a main advantage of the exp-ES scheme is that it preserves the control of the moments of the continuous model, assuming a superlinear growth condition on the drift coefficient (see Proposition 3.3 and Lemma 3.2).

Weak convergence and $C^{1,4}$-regularity of the Kolmogorov PDE associated to (1.1)

Our main result, stated in Proposition 3.3, exhibits an optimal theoretical rate of convergence of order one under hypotheses that are introduced in Section 3, and for bounded $C^{4}$ test functions.

Although some space of improvement are identified, the hypotheses in Section 3 are stated in order to balance the control moments of the exact and approximated processes with the moments and exponential moments required for the flow-derivative process used to establish the regularity of the Feynman-Kac formula. Indeed, the key point of the convergence rate analysis is to estimate the regularity of the solution to the backward Kolmogorov PDE associated with the representation $\mathbb{E}[f(X_{T}^{x})]$, where $(X_{t}^{x};0\leq t\leq T)$ denotes the flow of diffeomorphisms with initial condition $x>0$.

The technique presented in this paper for the analysis of the Kolmogorov PDE can be derived for a larger family of SDEs. Adapted from [4] which was dealing with the particular situation where $\tfrac{1}{2}<\alpha<1$, this methodology allows to control the successive derivatives of the Feynman-Kac representation up to the order four, by bypassing the difficulty of deriving the flow process more than one time, through a change of measure technique (see Sections 5.1 and 5.3 for dedicated results and details on this main point).

The paper is organized as follows. Conditions for the well-posedness of the generic SDE (1.1) as well as on the finiteness of the positive, negative and exponential moments of its solution are stated in Section 2. In Section 3 we construct the exponential-Euler scheme (1.4) and we present our main Proposition 3.3 on the weak rate of convergence when applied to model (1.1) and when applied to model (1.3) (Corollary 4.1). Section 4 presents some numerical experiments in order to show the effectiveness of the theoretical rate of convergence of the proposed exponential Euler scheme. We also compare exp-ES with the classical Symmetrized Euler scheme, the Symmetrized Milstein scheme, and the (stopped)-tamed scheme. Section 5 is devoted to the analysis of the regularity of the backward Kolmogorov PDE (Proposition 5.1) and Section 6 presents the proof of the weak error estimate.

For additional comments on the presented results and proofs, and complements on numerical experiments related to this work, we refer to [22, Chapters 1 and 2].

1.1 Notation

Throughout this paper, $T>0$ will refer to an arbitrary finite time horizon, $C$ will denote a positive constant, possibly depending on the parameters of the dynamic, which may change from line to line. Any process $(Z_{t},t\in[0,T])$ will be simply denoted $Z$. For any measurable function $f$, $x\mapsto f^{+}(x)$ denotes its positive part.

For any $a,b\in\mathbb{R}$, $a\vee b$ and $a\wedge b$ denote respectively the maximum and minimum between $a$ and $b$. Given the non-negative discrete time-step parameter $\Delta t$, we set $\eta(t)=\Delta t\lfloor\frac{t}{\Delta t}\rfloor$, and $\delta(t)=t-\eta(t)$. In order to shorten the writing of some expressions, we will use $\mathbb{E}^{\beta}[Z]$ for $\left(\mathbb{E}[Z]\right)^{\beta}$ and $f^{(k)}$ for the $k$-th derivative (whenever $k>1$) of a function $f:\mathbb{R}\rightarrow\mathbb{R}$.

We introduce the notion of locally Lipschitz continuity property used in this paper from the formalism previously used in [14] and [7]. The following definition specifies power-dependencies involved in the local Lipschitz factor.

Definition 1.1**.**

Let $f$ be a real-valued function and $\textsl{dom}(f)$ denotes its definition domain. We say that $f$ is $(\overline{\gamma},\underline{\gamma})$-locally Lipschitz if there exists a non-negative constant $C$, $\overline{\gamma}$ and $\underline{\gamma}$ such that

[TABLE]

When $\underline{\gamma}=0$, $f$ is said to be $\overline{\gamma}$-locally Lipschitz continuous and

[TABLE]

With this definition, a Lipschitz function is [math]-locally Lipschitz, and it is included in the set of all $\gamma$-locally Lipschitz functions, $\gamma\geq 0$.

The following lemma formalize the straightforward link between the locally Lipschitz property of a function and its derivative.

Lemma 1.2**.**

Let $f$ be a real-valued function, continuously differentiable, with $f^{\prime}$ being $(\overline{\beta},\underline{\beta})$-locally Lipschitz continuous in the sense of Definition 1.1. Then, $f$ is $(\overline{\alpha},\underline{\alpha})$-locally Lipschitz continuous with $\overline{\alpha}\leq\overline{\beta}+1$, and $\underline{\alpha}\leq\underline{\beta}$.

2 Strong wellposedness for the solution to SDE (1.1)

Some sufficient conditions ensuring the strong well-posedness and control of moments of the solution to (1.1) are now exhibited. Hereafter we assume the following hypotheses on the SDE (1.1):

Hypothesis H1**.**

$\alpha>1$, $\sigma>0$, and the (deterministic) initial condition $x>0$.

Hypothesis H2**.**

The drift $b$ is $2(\alpha-1)$-locally Lipschitz continuous (in the sense of Definition 1.1) and $b(0)\geq 0$.

Hypothesis H3**.**

There exist some finite constants $B_{i}\geq 0$, with $i=1,2$, such that, for all $x\geq 0$,

[TABLE]

The proofs of the following proposition and lemmas are technical and by itself not directly in relation with the convergence analysis of our scheme. They are given in Appendix C.

Proposition 2.1**.**

Assume H1, H2 and H3. Then there exists a unique (strictly) positive strong solution $X$ to the SDE (1.1). In addition, for all exponent $p$ such that $0\leq 2p\leq 1+\tfrac{2B_{2}}{\sigma^{2}}$, we have

[TABLE]

where the non-negative constant $C_{p}$ may depend on $p$, but does not depend on $x$.

Remark 2.2**.**

In all the results of this section, Assumption H2 can be weaken to a $\gamma$-locally Lipschitz property on $b$, provided that $\gamma\geq 2(\alpha-1).$ The limit case $\gamma=2(\alpha-1)$ is considered in order to simplify the computation. The effect of the non-negative constants $B_{1}$ and $b(0)$ on the SDE (1.1) is to push the solution away from zero, while the constant $B_{2}$ enables to counter the growth of the solution due to the diffusion term. The assumption on $B_{2}\geq 0$ can so be weaken into $B_{2}\geq-\tfrac{\sigma^{2}}{2}$ in order to preserve the well-posedness of SDE (1.1). However, in that case, only the moments of order strictly less than one are controlled.

We complement the preceding proposition by the control of the negative and exponential moments.

Lemma 2.3**.**

Assume H1, H2 and H3. Then, for all $q>0$, the solution $X$ to (1.1) satisfies

[TABLE]

where the non-negative constant $C_{q}$ may depend on $q$, but does not depend on $x$.

Lemma 2.4**.**

Assume H1, H2, and H3. Then, when $b(0)=0$, there exists a constant $C$, independent on $x$, such that for all $\mu\leq\tfrac{(\sigma^{2}+2B_{2})^{2}}{8\sigma^{2}}$,

[TABLE]

Otherwise, when $b(0)>0$, if we assume in addition that $\alpha>\frac{3}{2}$, then for all $\mu<B_{2}\sigma^{2}$,

[TABLE]

In the above upper-bounds, the non-negative constant $C$ ca be bounded uniformly in $\mu$.

3 The exponential Euler scheme and its rate of convergence

Given the possible non-integer power value for $\alpha$ in the diffusion term, we seek for an appropriate numerical approximation preserving the positiveness of the process and exponential form is a good candidate for this purpose. By rewriting the SDE (1.1) as

[TABLE]

and given $\{t_{0}=0,t_{1},\ldots,t_{N-1},t_{N}=T\}$, a $N$-partition of the time interval $[0,T]$ with time-step $\Delta t=t_{n+1}-t_{n}$, we consider first the approximation $(\widehat{X}_{t_{n}},n\geq 1)$ given by

[TABLE]

and its continuous version given by the interpolation in time:

[TABLE]

where $\eta(t):=\sup\{t_{i}:t_{i}<t\}$. Ensuring the strict positivity of the approximation at all times, the scheme (3.1) enables also to counterbalance the rapid growth of the diffusion $\widehat{X}^{\alpha-1}_{\eta(t)}$ by the drift contribution $\tfrac{b(\widehat{X}_{\eta(t)})}{\widehat{X}_{\eta(t)}}$ subject to H3 and H2. The scheme (3.1) is also sensitive to the value of $b$ near zero: when $b(0)=0$, H3 yields to

[TABLE]

and, combined with H2, enables to prove the existence of some positive moments for $(\widehat{X}_{t};0\leq t\leq T)$ (replicating for instance the last proof steps of Proposition 2.1). But when $b(0)>0$, numerical instabilities can be observed when $\widehat{X}$ comes close to zero. More specifically, we haven’t been able to find a threshold $\xi$ such that $\mathbb{P}(\widehat{X}_{t}\leq\xi)$ decays in $\Delta t$, nor to control some positive moments in that case.

To overcome such instabilities, the continuous version of the scheme (3.1) can be modified by adding and subtracting $b(0)$ as follows:

[TABLE]

or equivalently, defining ${\delta(t)}:=t-{\eta(t)}$,

[TABLE]

for which we need to discretize the integral appearing in the right hand-side to turn it in a numerical algorithm. The approximation

[TABLE]

makes the corresponding scheme comes back to (3.1) for which we do not control –a priori– positive moments. In contrast, the approximation

[TABLE]

produces the following implicit numerical scheme:

[TABLE]

where $f(t,x)=x\exp\{-\frac{b(0){\delta(t)}}{x}\}$, for which control of positive moments for $(\check{X}_{t};0\leq t\leq T)$ are obtained under H3 (see [22, Chap. 2]).

To turn (3.2) in a numerical scheme, we combine it with an approximation method for $x\mapsto f^{-1}(x)$, by considering a Taylor expansion of first order:

[TABLE]

With this last approximation, we define the scheme $(\overline{X}_{t_{n}},n\geq 1)$, that we now refer to as exp-ES, for * Exponential-Euler Scheme*, by $\overline{X}_{0}=x$, and

[TABLE]

admitting the continuous version

[TABLE]

driven by the SDE:

[TABLE]

Remark 3.1**.**

By construction, due to the exponential form in (3.4), $\overline{X}_{t}-b(0){\delta(t)}>0$, provided that $x>0$. In particular, for all $0\leq t\leq T$, $\overline{X}_{t}>0$ .

For the exp-ES $\overline{X}$, we bound the same order of $2p$ th-moments than for $X$ in Proposition 2.1.

Lemma 3.2**.**

Assume H1, H2 and H3. For all exponent $0<2p\leq 1+\frac{2B_{2}}{\sigma^{2}}$, there exists a non-negative constant $C_{p}$, depending on $p$ but not on $x$, such that

[TABLE]

Proof.

Applying Itô’s formula to the process $\overline{X}$ (for simplification, we omit the localization argument previously used in the proof of Proposition 2.1), we get

[TABLE]

The proof ends by applying Gronwall’s inequality. ∎

3.1 Main results

Under the following hypotheses we state below the weak rate of convergence of order one for (1.1) associated with the exponential-Euler scheme (3.3):

Hypothesis H2’****.

(For the regularity of the Kolmogorov PDE (5.3).)* The function $b$ is $2(\alpha-1)$-locally Lipschitz, and $b(0)\geq 0$. In addition, $b$ is of class $\mathcal{C}^{4}(\mathbb{R}^{+})$, with derivatives $b^{(i)}$ being $({\overline{\gamma}}_{(i)},{\underline{\gamma}}_{(i)})$-locally Lipschitz continuous, for $i=1,\ldots,4$.*

H2’ implies in particular that

[TABLE]

where, according to Lemma 1.2, $2\alpha-3\leq\overline{\gamma}_{(1)}$, $\overline{\gamma}_{(i)}\leq\overline{\gamma}_{(i+1)}+1$, and $\underline{\gamma}_{(i)}\leq\underline{\gamma}_{(i+1)}$ for $i=1,2,3$.

Hypothesis H3’****.

(For the exponential moments of $X$.)* There exist a set of constants $B_{i},B_{i}^{\prime}\geq 0$, with $i=1,2$, such that*

[TABLE]

Hypothesis H4**.**

(For the weak convergence rate derivation.)* The powers $\underline{\gamma}_{(i)}$ in H2’ satisfy:*

[TABLE]

Hypothesis H5**.**

(For the regularity of the Kolmogorov PDE in Proposition 5.1).* The constants $B_{2}$, $B_{2}^{\prime}$, $\alpha$ and the $(\overline{\gamma}_{(i)},i=1,\ldots,4)$ in H2’ and H3’ satisfy*

[TABLE]

where $\overline{\overline{\beta}}=3(\overline{\gamma}_{(2)}+1)\vee(\overline{\gamma}_{(2)}+\overline{\gamma}_{(3)}+2)\vee(\overline{\gamma}_{(4)}+1)$.

In addition whenever $b(0)>0$, we assume that $\alpha>\tfrac{3}{2}$ and we modify the constraint on $B_{2}$ as

[TABLE]

Proposition 3.3**.**

Let $f$ be a bounded $\mathcal{C}^{4}(\mathbb{R}^{+})$ function with bounded derivatives up to order 4. Consider the process $X$ solution to the SDE (1.1) with deterministic initial condition $x>0$, together with its approximation $\overline{X}$ in (3.4). Assume H1, H2’,H3’, H4, and H5. Then, there exists a constant $C>0$ depending on the parameters $B_{i},B_{i}^{\prime},\alpha,\sigma$ and possibly on $T$ and $x$, but independent on $\Delta t$, such that

[TABLE]

The hypotheses in Proposition 3.3 are all sufficient conditions, considered in order to simplify the analysis of the regularity associated with the solution of the backward Kolmogorov PDE.

Precisely, Hypotheses H3’ and H2’ are considered in order to obtain polynomial bounds for the derivatives of the solution to the backward Kolmogorov PDE (Proposition 5.1). Later, in the computation of the weak error, we use H5 specifically to control the resulting positive moments of the exp-ES process (see the proof of Proposition 3.3 in Section 6), and by considering H4 we seek to avoid the need to control the negative moments of the approximation scheme arising also from the estimation of these derivatives (see for instance the inequality (6.6) below).

We also emphasis that the analysis exposed in this paper can be easily adapted to the case of locally bounded $\mathcal{C}^{4}(\mathbb{R}^{+})$-function $f$ with locally bounded derivatives.

4 Numerical Experiments

This section illustrates with some experiments the theoretical rate of convergence in Proposition 3.3. More experiments are shown in Appendix D and explore the hypotheses on the coefficient parameters by testing the convergence through different cases, fulfilling or not some of the hypotheses.

In particular we illustrate the fact that hypothesis H5 do not correspond to a necessary condition on the parameters $B_{2},B_{2}^{\prime},\alpha,\sigma$ involved in the model. First we restrict the set of parameters by considering the explicit model

[TABLE]

for which $B^{\prime}_{1}=1$ and $B^{\prime}_{2}=(2\alpha-1)B_{2}$. Proposition 3.3 can be shapely adapted in this particular situation as follows:

Corollary 4.1**.**

We consider the solution $X$ to (4.1). When $B_{0}=0$, assume $\alpha>1$ and

[TABLE]

When $B_{0}>0$, assume in addition that $\alpha>\tfrac{3}{2}$ and

[TABLE]

Then, for $f\in\mathcal{C}_{b}^{4}(\mathbb{R}^{+})$, there exist a constant $C>0$ depending on the parameters $B_{i},\alpha,\sigma$ and possibly on $T$ and $x$, but independent on $\Delta t$, such that

[TABLE]

Numerical parameters.

For all the presented numerical experiments, we consider a unit terminal time $T=1$, the initial condition $x=1$ and the time step $\Delta t=1/2^{p}$, for $p=1,\ldots,9$. In addition, the empirical mean of the scheme approximating $\mathbb{E}f(\overline{X}_{T})$ is estimated by a Monte Carlo approximation, involving $n=10^{5}$ independent trajectories.

Test functions.

Along this section, we consider four different test functions, not all bounded,

[TABLE]

Model cases.

Denoting $\kappa$ as the left-handside of (4.2) or (4.3), we consider the following cases, determined by the data $(B_{0},B_{1},B_{2},\sigma,\alpha)$

[TABLE]

with two of them, Cases 2 and 3, that are not satisfying H5.

Moreover, Case 3 satisfies assumptions of Theorem 2.1 in [17] that states that the approximated moments by the Euler-Maruyama scheme and the strong $L^{p}$-error associated to moment-approximations diverges. Also in [17], the authors prove the divergence in the weak sense for the $p$-th moments of the Euler-Maruyama scheme in that case.

From Proposition 2.1 and Lemma 2.3, the triplets $(B_{2},\sigma,\alpha)$ in Case 1 to Case 4, guarantee the finiteness of the expectation $\mathbb{E}f(X_{T})$ for each test function in the set $\{x,x^{2},\exp(-x^{2})\}$.

Computation of the reference values.

For both test functions $f(x)=x$ and $f(x)=x^{2}$, reference values of $\mathbb{E}[X_{T}]$ and $\mathbb{E}[X_{T}^{2}]$ are computed analytically for Case 1 to Case 3, (see details in [22, Chapter 1] and Appendix D). For the others cases, the reference values are computed based on a $n_{0}$-Monte Carlo method combined with the scheme (3.1), $n_{0}=10^{7}$ and $\Delta t_{\text{ref}}=2^{-14}$:

[TABLE]

Numerical results are shown in Table 1, where we can observe the rate of convergence of order one, except for Case 3 with test function $f(x)=x^{2}$, in all the rows corresponding to the selection of bounded/unbounded test functions: the error is divided by 2 when going from left to right, even if some saturation can be observed for the smallest error values $(p=8,9)$ when Monte Carlo error starts to be dominant.

This behavior is also illustrated in Figure 1, plotting the obtained results in a log-log scale. This confirms that our proofs can certainly be extended for a larger class of test functions, and model parameters. In particular, we highlight Case 3 that converges weakly with order one for $f(x)=x,\exp\{-x^{2}\}$, even if H5 is not fulfilled, and even moreover we know that the classical Euler-Maruyama scheme is strongly diverging (as stated in [17]) in this case.

Additional numerical results are proposed in Appendix D where, in particular, some comparisons on the exp-ES with other dedicated numerical schemes are presented.

5 Analysis of the backward Kolmogorov PDE related to (1.1)

This section is devoted to the regularity analysis on the solution of the backward Kolmogorov PDE related to (1.1). Stochastic analysis is used here to establish key estimates on the solution of the PDE. We consider the flow process $(X_{t}^{x};0\leq t\leq T)$ starting from $x>0$:

[TABLE]

According to the Feynman-Kac representation theorem (see e.g. [20, Chap. V]), provided that $x\mapsto X_{T-t}^{x}$ and a given $f$ are smooth enough, the function

[TABLE]

is a natural candidate to be the classical solution to the backward Kolmogorov PDE:

[TABLE]

Proposition 5.1**.**

Let $f\in\mathcal{C}^{4}_{b}(\mathbb{R}^{+})$. Then, assuming H1,H2’,H3’ and H5, the function $u(t,x)$ in (5.2) is a solution to the PDE (5.3) of class $\mathcal{C}^{1,4}([0,T]\times\mathbb{R}^{+})$, and there exists a finite constant $C$ such that

[TABLE]

where $\overline{\gamma}_{(i)}$ and $\underline{\gamma}_{(i)}$ are as in H2’, $\overline{\overline{\beta}}$ as in H5, and $\overline{\beta}$, $\underline{\beta}$ are $\underline{\underline{\beta}}$ are given by

[TABLE]

5.1 Main lines for the proof of Proposition

5.1

The proof of Proposition 5.1 follows the methodology used in [4] that combines some adequate successive changes of measure in the Feynman-Kac formula for $u$ and in its derivatives in order to kill some unsuitable term in the obtained expression for $\frac{\partial^{i}u}{\partial x^{i}}$ before to derive it again. More precisely, the main hypothesis $\alpha>1$ allows to derive at most one time the diffusion coefficient $x\mapsto x^{\alpha}$, and by extension $x\mapsto(X^{x}_{t};\,0\leq t\leq T)$, before to potentially produce some negative power term for higher order derivative. In contrast, Lemma 2.4 embed a local Novikov condition which allows to control the exponential martingale for the first derivative of the diffusion only (the power $2(\alpha-1)$ corresponding to the quadratic variation of resulting from this derivative). We present here, briefly and formally, how we can combine derivatives and change of measure to overcome higher order derivative of the diffusion before detailing the proof in the rest of this section and in the Appendix B.

Following [4], we introduce the family of processes $X^{x}(\lambda)$, with parameter $\lambda>0$, as the solution of the SDE:

[TABLE]

For each $\lambda>0$, Equation (5.5) can be seen as a modification of (5.1) with a drift component $b^{\lambda}$ given by

[TABLE]

where $B_{2}^{\lambda}:=B_{2}-\lambda\sigma^{2}$. Due to the locally Lipschitz property of the coefficients (and their derivatives) in (5.5), the process $X^{x}(\lambda)$ is continuously differentiable w.r.t. $x$ (see Protter [25, Thm V.39]) for $\mathbb{P}$-almost all $\omega\in\Omega$. Therefore, we can define the derivative of the flow with respect to the initial condition $x$:

[TABLE]

as the solution to the SDE, for $t\in(0,T]$,

[TABLE]

Whenever the process $(\int_{0}^{t}\left({X_{s}^{x}}(\lambda)\right)^{\alpha-1}dW_{s};~{}0\leq t\leq T)$ is a square integrable martingale, $J^{x}(\lambda)$ admits the following exponential form (see e.g. [25, Theorem V.52])

[TABLE]

Now, we may identify the first order derivative of $u(t,x)=\mathbb{E}[f(X^{x}_{T-t}(0))]$ as:

[TABLE]

Before computing the second derivative, we change the measure in the expectation above in order to eliminate $J_{T-t}^{x}(0)$ and avoid the problem of the a priori control of $\tfrac{d}{dx}J_{T-t}^{x}(0)$: Consider the Radon-Nikodym density

[TABLE]

for $(B_{t}^{\alpha};\,0\leq t\leq T)$ the standard $\mathbb{Q}^{\alpha}$-Brownian motion given by ${B_{t}^{\alpha}=W_{t}-\alpha\sigma\int_{0}^{t}(X_{s}^{x}(0))^{\alpha-1}ds}$. Then

[TABLE]

From the explicit form of the process $J^{x}(0)$, we recognize

[TABLE]

and hence, J_{T-t}^{x}(0)\,\mathcal{Z}_{T-t}^{(0,\alpha)}=\exp\big{\{}\int_{0}^{T-t}b^{\prime}(X_{s}^{x}(0))ds\big{\}}. Moreover, from the identification

$\textit{Law}^{\mathbb{Q}^{\alpha}}(X^{x}(0))=\textit{Law}^{\mathbb{P}}(X^{x}(\alpha))$, we can rewrite $\tfrac{\partial u}{\partial x}$ as

[TABLE]

which is continuously differentiable in $x$, with a derivative that depends on the derivative of $f^{\prime}$ and $b^{\prime}$ only. The same procedure can be iterated for other higher order derivative of $x\mapsto u(t,x)$. The price to pay here resides in the assumptions needed to get, at each iteration step, the estimates stated in Proposition 2.1 and Lemmas 2.3 and 2.4, which are more constraining to satisfy according to the increasing value $\lambda$ introduced after each successive changes of measure, and corresponding to successive derivatives. Typically, while we add some unbounded term $\lambda\sigma^{2}x^{2\alpha-1}$ to $b$, the constant $B_{2}$ in H3, ensuring the wellposedness for (1.1), must be strengthened to also ensure the wellposedness and the finiteness of the moments of the solution to (5.5), as well as some moments of $J^{x}(\alpha)$. This strengthening on $B_{2}$ is summarized in the following corollary, combining the results of Proposition 2.1 and Lemmas 2.3 and 2.4:

Corollary 5.2**.**

Assume H1, H2 and H3. Then, for any $\lambda>0$ such that $\lambda<\tfrac{B_{2}}{\sigma^{2}}$, (and so $B_{2}^{\lambda}=B_{2}-\lambda\sigma^{2}>0$), there exists a unique positive strong solution to (5.5). For all $q>0$ and $p>0$ such that $p\leq\frac{1}{2}+\frac{B^{\lambda}_{2}}{\sigma^{2}}$, this solution further satisfies:

[TABLE]

The non-negative constants $C_{p}$, $C_{q}$ do not depend on $x$, but $C_{p}$, respectively $C_{q}$, may depend on $p$, respectively $q$.

If $b(0)=0$, then for all $\mu\leq\frac{(\sigma^{2}+2B_{2}^{\lambda})^{2}}{8\sigma^{2}}$,

[TABLE]

The estimate (5.11) still holds when $b(0)>0$ under the restriction that $\alpha>\frac{3}{2}$, and $\mu<B_{2}^{\lambda}\sigma^{2}$.

In the rest of Section 5, we address rigorously the formal steps described above, validating first the representation (5.8) in the next section, detailing the change of measure leading to (5.9), before finally completing the proof of Proposition 5.1

5.2 Interchanging derivative and expectation

The (sufficient) conditions for the possible interchange between expectation $\mathbb{E}$ and derivative with respect to the initial condition $\tfrac{\partial}{\partial x}$ are stated in following:

Proposition 5.3**.**

Let $\Phi,g,h\in\mathcal{C}([0,+\infty))$ some continuously differentiable functions such that $\Phi$ is bounded and has bounded first derivative in $[0,+\infty)$, $g$ is bounded from above, and the functions $g^{\prime}$ and $h$ satisfy the following growth conditions

[TABLE]

[TABLE]

for all $x>0$ and some non-negative constants $\rho_{i}$, with $i=0,1,\ldots,5$. Assume H1, H2’, H3’ and

[TABLE]

Then for any $\lambda\geq 0$, the function defined by

[TABLE]

is continuously differentiable in $x$, with

[TABLE]

The proof of Proposition 5.3 can be summarized as follows: the goal is to show (omitting the dependence to $\lambda$ in the notation during a few lines) that $\mathbb{E}[\tfrac{\Phi(X_{t}^{x+\epsilon})-\Phi(X_{t}^{x})}{\epsilon}]$ tends to $\mathbb{E}[\Phi^{\prime}(X_{t}^{x})J_{t}^{x}]$ when $\epsilon$ tends to [math]. Introducing the process

[TABLE]

we start with the following decomposition

[TABLE]

The convergence $\mathbb{E}[\left|J_{t}^{x,\epsilon}-J_{t}^{x}\right|]\rightarrow 0$, when $\epsilon$ tends to 0, is thus a crucial step of the proof and we first establish this result. The detailed proof of the proposition is postponed, after having obtained some dedicated estimates on processes $J^{x,\epsilon}(\lambda)$ and $J^{x}(\lambda)$ in the following subsection.

5.2.1 Preliminary estimations

The process $J^{x,\epsilon}(\lambda)$ given by $J_{t}^{x,\epsilon}(\lambda)=\tfrac{X_{t}^{x+\epsilon}(\lambda)-X_{t}^{x}(\lambda)}{\epsilon}$ satisfies the linear SDE

[TABLE]

where we have defined

[TABLE]

As $J_{t}^{x,\epsilon}(\lambda)>0$ a.s., these auxiliary processes may also write

[TABLE]

The two following lemmas, whose proofs are postponed to Appendix A, assert the finiteness of the moments of the processes $J^{x,\epsilon}(\lambda)$ and $J^{x}(\lambda)$ as well as the $L^{p}$ continuity of $x\mapsto X^{x}_{t}(\lambda)$ (Lemma 5.4) and the $L^{2}$ continuity of $x\mapsto J^{x}_{t}(\lambda)$ in Lemma 5.5.

Lemma 5.4**.**

Assume H1, H2, H3’ and

[TABLE]

Then, for all $q>0$ such that $(2\alpha-1)\lambda+\tfrac{\alpha^{2}}{2}(q-1)\leq\tfrac{B_{2}^{\prime}}{\sigma^{2}}$, the processes $J^{x,\epsilon}(\lambda)$ and $J^{x}(\lambda)$, respective solutions to (5.13) and (5.6), satisfy

[TABLE]

and

[TABLE]

Lemma 5.5**.**

Assume H1, H2’, H3’ and

[TABLE]

for the constants $B_{2},B_{2}^{\prime},\sigma,\alpha,\overline{\gamma}_{(1)}$ as in H2’ and H3’. Then, for all $t\in[0,T]$,

[TABLE]

5.2.2 Proof of Proposition 5.3

To simplify notation we omit the dependence on $\lambda$ in the processes and prove the result for $h\equiv 0$. In order to prove the interchange between expectation and operator $\frac{\partial}{\partial x}$ we must show the equality

[TABLE]

introducing the process $(Y_{t}^{x}:=\int_{0}^{t}g(X_{s}^{x})ds;\;0\leq t\leq T)$, with derivative $\frac{dY_{t}^{x}}{dx}=\int_{0}^{t}g^{\prime}(X_{s}^{x})J_{s}^{x}ds$. Following the decomposition (5.12), we rewrite the difference

[TABLE]

and we analyze separately the limit, when $\epsilon$ tends to 0, of each term in the right-hand side of (5.2.2). Notice that $g$ being bounded from above, all the exponential terms above are bounded.

Since $|\Phi^{\prime}|$ and $\mathbb{E}[|J_{t}^{x}\exp\{Y_{t}^{x}\}|]$ are bounded, the second term $\mathbb{E}[|B^{\epsilon}_{t}|]$ is uniformly integrable. Moreover, $\Phi^{\prime}$ is continuous, and, under the hypotheses ($\Lambda$.1), according to Lemma 5.4, $X_{t}^{x+\epsilon}$ converges in $L^{4}$ and then in probability to $X_{t}^{x}$ when $\epsilon$ tends to 0. Therefore, we can apply the Lebesgue dominated convergence Theorem, obtaining (up to a subsequence still denoted by $\epsilon$) that $\lim_{\epsilon\rightarrow 0}\mathbb{E}[|B^{\epsilon}_{t}|]=0$.

Similarly, since $|\Phi|$ and $\mathbb{E}[|\tfrac{dY_{t}^{x}}{dx}|]\leq\mathbb{E}[\int_{0}^{t}|g^{\prime}(X_{s}^{x})J_{s}^{x}|ds]$ are bounded according to $\mathbb{E}[|(J_{s}^{x})^{2}|]$ and $\mathbb{E}[|(X_{s}^{x})^{2\rho_{2}}|]$, the sequence $D^{\epsilon}_{t}$ is uniformly integrable. Then, by the convergence in probability of $Y_{t}^{x+\epsilon}$ towards $Y_{t}^{x}$ obtained from the following bound

[TABLE]

we obtain (again up to a subsequence) that $\lim_{\epsilon\rightarrow 0}\mathbb{E}[|D^{\epsilon}_{t}|]=0$.

For $\mathbb{E}[|A^{\epsilon}_{t}|]$, we use Cauchy-Schwartz inequality, from which we get

[TABLE]

Therefore, by Lemma 5.5 and Lebesgue’s Theorem, $\mathbb{E}[|A^{\epsilon}_{t}|]$ converges to 0 when $\epsilon$ tends to 0.

Finally, for $\mathbb{E}[|C^{\epsilon}_{t}|]$,

[TABLE]

and we conclude that $\lim_{\epsilon\rightarrow 0}\mathbb{E}[|C^{\epsilon}_{t}|]=0$ with similar arguments to those used for $A^{\epsilon}_{t}$.

Coming back to (5.2.2), we obtain the convergence to zero, up to a subsequence $\{\epsilon_{k}\}_{k\geq 1}$, of its right-hand side. Therefore, by uniqueness of the limit, we deduce the convergence in (5.16).

The condition ($\Lambda$.1) makes the intersection of Conditions ($\Lambda$.2) and ($\Lambda$.3) with the highest $p$-moment order needed to obtain the convergence in (5.16). When $h$ is not reduce to zero, it is sufficient for the proof to control

[TABLE]

by the moments $\sup_{0\leq t\leq T}\mathbb{E}[|X_{t}^{x}|^{2(\rho_{2}+\rho_{4})}+|X_{t}^{x}|^{2\rho_{0}}]$ which are finite when $\max\{\rho_{0},\rho_{2}+\rho_{4}\}\leq\frac{1}{2}+\frac{B^{\lambda}_{2}}{\sigma^{2}}$.

5.3 Change of measure

Considering a generic expression coming from the application of Proposition 5.3, typically of the form

[TABLE]

we introduce the change of probability measure that allows to remove the term $J_{t}^{x}(\lambda)$ in the expression above.

Let us consider the process $(B_{t}^{\lambda+\alpha};\;0\leq t\leq T)$ defined as $B_{t}^{\lambda+\alpha}=W_{t}-\alpha\sigma\int_{0}^{t}(X_{s}^{x}(\lambda))^{\alpha-1}ds$. Then, using Lemma 5.6 below and Girsanov’s Theorem, we can construct the probability measure $\mathbb{Q}^{\lambda+\alpha}$ under which $(B_{t}^{\lambda+\alpha};\;0\leq t\leq T)$ is a standard Brownian motion, by introducing the Radon-Nikodyn density

[TABLE]

Lemma 5.6 below gives a sufficient condition for the process $\mathcal{Z}^{(\lambda,\lambda+\alpha)}$ to be martingale for a given $\lambda\geq 0$. From the explicit form of the process $J_{t}^{x}(\lambda)$ in (5.7), we recognize

[TABLE]

Hence, by H3’ and ($\Lambda$.1),

[TABLE]

Moreover, we can easily check the identity $\textit{Law}^{\mathbb{Q}^{\lambda+\alpha}}\left(X^{x}\left(\lambda\right)\right)=\textit{Law}^{\mathbb{P}}\left(X^{x}\left(\lambda+\alpha\right)\right)$, so that $\tfrac{\partial v}{\partial x}$ can be rewrite

[TABLE]

The following lemma is a direct consequence of the Novikov’s criterion whose fulfillment is ensured by applying Corollary 5.2.

Lemma 5.6**.**

Assume H1, H2 and H3’. Assume in addition that $B_{2},\alpha$ and $\sigma$ in H3’ satisfy

[TABLE]

Then the process $(M_{t}^{x}(\lambda);\;0\leq t\leq T)$ defined by

[TABLE]

is a $\mathbb{P}$- martingale.

5.4 Proof of Proposition 5.1

Let us first note that the hypotheses considered in Proposition 5.1, in particular H5, allow to apply Proposition 5.3 up to $\lambda=3\alpha$.

Estimates on $u$ and $\tfrac{\partial u}{\partial t}$.

The uniform boundedness of $u(t,x)=\mathbb{E}[f(X^{x}_{T-t}(0))]$ is an immediate consequence of the bounded of $f$. Applying Itô’s formula and since $X^{x}(0)$ has finite $2\alpha$-th moment,

[TABLE]

From this expression, we can deduce that $\tfrac{\partial u}{\partial t}$ is continuous in $[0,T]\times\mathbb{R}^{+}$ with

[TABLE]

where $C$ is a positive constant depending on $\alpha,\sigma,b,\|f^{(i)}\|_{\infty}$ for $i=0,1,2$. The $2(\alpha-1)$-locally Lipschitz continuity of the drift $b$ gives us

[TABLE]

Applying Corollary 5.2 for the control of $\sup_{t\in[0,T]}\mathbb{E}|X_{t}^{x}|^{2\alpha}]$ granted by the condition $\alpha\leq\frac{1}{2}+\frac{B_{2}}{\sigma^{2}}$ in H5, we obtain

[TABLE]

Estimates on $\tfrac{\partial u}{\partial x}$ and $\tfrac{\partial^{2}u}{\partial x^{2}}$.

The differentiability up to order $4$ of $x\mapsto u(t,x)$ will rely on the iterative use of Proposition 5.3 and on rewriting (whenever its is necessary) the function $\frac{\partial^{j}u}{\partial x^{j}}(t,x)$, $j=0,1,2,3$, as

[TABLE]

for some continuous differentiable functions $g_{j}$ and $h_{j}$ with locally bounded spatial derivatives in $[0,+\infty)$, with $g_{j}$ bounded from above, and $f_{j}$ some bounded continuously differentiable functions with bounded derivative.

In order to prove the identity (5.9) for $\tfrac{\partial u}{\partial x}$, we apply Proposition 5.3 for $f_{0}=f$ and $g_{0}=h_{0}\equiv 0$ and a first change of measure. So we need Condition ($\Lambda$.1) to be satisfied for $\rho_{i}=0$ and $\lambda=0$ and we need also the hypotheses of Lemma 5.6 satisfied for $\lambda=0$. From (5.9), we immediately get that

[TABLE]

Next from (5.9), we identify the form (5.4) with $f_{1}=f^{\prime}$, $g_{1}=b^{\prime}(x)$ (bounded from above, and with $|g^{\prime}_{1}(x)|\leq C(1+|x|^{\overline{\gamma}_{(2)}+1}+|x|^{-\underline{\gamma}_{(2)}})$ and $h(x)\equiv 0$. Applying Proposition 5.3 again (with ($\Lambda$.1) to be satisfied for $\rho_{2}={\overline{\gamma}_{(2)}+1}$ and $\lambda=\alpha$) we obtain that $\frac{\partial u}{\partial x}$ is continuously differentiable in $x$ with derivative given by

[TABLE]

Notice that by means of the Markov property and time homogeneity of the process $(X_{s}^{x}(\alpha);0\leq s\leq T-t)$ we have

[TABLE]

where the last equality is obtained from (5.9). Then, we get

[TABLE]

Substituting the last equality in (5.22),

[TABLE]

Introducing the change of measure \tfrac{d\mathbb{Q}^{2\alpha}}{d\mathbb{P}}\big{|}_{\mathcal{F}_{t}}:=\frac{1}{\mathcal{Z}_{t}^{(\alpha,2\alpha)}} where $\mathcal{Z}^{(\alpha,2\alpha)}$ is defined in (5.18) (under the conditions (5.19) and ($\Lambda$.1) applied for $\lambda=\alpha$), we can observe that (assuming $2B^{\prime}_{2}\geq\alpha(2\alpha-1)\sigma^{2}$ in ($\Lambda$.1))

[TABLE]

So changing the measure in (5.24), with the observation that $\textit{Law}^{\mathbb{Q}^{2\alpha}}(X^{x}(\alpha))=\textit{Law}^{\mathbb{P}}(X^{x}(2\alpha))$, and by the boundedness of the functions $\frac{\partial u}{\partial x}$ and $f^{(2)}$ we obtain

[TABLE]

We then apply Corollary 5.2 to conclude on the estimate for the second derivative:

[TABLE]

under the condition that $\overline{\gamma}_{(2)}+1\leq 1+\frac{2B_{2}^{2\alpha}}{\sigma^{2}}$.

We end this proof by iterating the derivative estimations. This (technical) step is postponed to Appendix B.

6 Proof of Proposition 3.3

Introducing the notation

[TABLE]

for which we have

[TABLE]

we can rewrite the dynamics (3.5) as

[TABLE]

We associate to it, the differential operator

[TABLE]

Then, applying Itô’s formula to the $C^{1,4}$ function $u$ along $\overline{X}$ in the time interval $[0,T]$, we obtain

[TABLE]

Lemma 3.2 under H5 allows to control the $2p$-th moments of the exp-ES process $\overline{X}_{t}$ up to the order $2p:=6\alpha+2\overline{\overline{\beta}}\vee(\overline{\overline{\beta}}+2\alpha)$. By Proposition 5.1 we have for ach $k=1,\ldots,N$

[TABLE]

Moreover, since $u$ is solution to the Cauchy problem (5.3), we decompose the error in two contributions:

[TABLE]

Notice that the functions $\overline{b}$ and $\overline{\sigma}$ are continuously differentiable with respect to $x$, and piecewise continuously differentiable with respect to $t$ on each subintervals $[t_{k},t_{k+1})$. These linear functions in $x$ and $t$ produce constant values as derivatives, only parametrized by $y$: for all $t\in(0,T)$, for all $(\theta,x,y)\in(\eta(t),t)\times(0,+\infty)\times(0,+\infty)$,

[TABLE]

Then, observing that $I_{i}(\eta(s),\eta(s))=0$, we apply Itô’s formula a second time in the interval $[\eta(s),s]$ and we obtain the two following decompositions for each $s\in[t_{k-1},t_{k}]$ with $k=1,\ldots,N$:

[TABLE]

[TABLE]

We use again the backward Kolmogorov PDE (5.3) to compute the time derivatives

[TABLE]

From H2’ we have for all $x\geq 0$ and $i=1,2,3,4$

[TABLE]

The key of the proof is to upper bound each $I^{k}_{j}$ by combining the estimates of the derivatives of $u$, the polynomial growth of the drift and diffusion coefficients and its derivatives, with upper-bounds of moments of the exp-ES process obtained from Lemma 3.2. By considering the following Young inequality for arbitrary $m,n\geq 0$,

[TABLE]

we get

[TABLE]

for all $j=1,2$, $k=1,2,4$ and some $\beta_{k,j}\in[0,2\beta]$. Then, substituting (6.5) in (6) we recover the rate of order one for the weak approximation error:

[TABLE]

We detail the analysis of the first term $|I_{1}^{1}|$:

[TABLE]

where, from Proposition 5.1 and (3.6), (6.3), we have

[TABLE]

and

[TABLE]

Remaining with the biggest $\pm$ exponent (using from Lemma 1.2 that $\underline{\gamma}_{(i)}\leq\underline{\gamma}_{(i+1)}$, $i=1,2,3$),

[TABLE]

where $\overline{\beta}=2(\overline{\gamma}_{(2)}+1)\vee(1+\overline{\gamma}_{(3)})$ and $\underline{\beta}=2\underline{\gamma}_{(2)}\vee\underline{\gamma}_{(3)}\vee(\underline{\gamma}_{(2)}+3-2\alpha)$. Therefore, we get

[TABLE]

Since we do not have a priori control of negative moments of $\overline{X}$, we must impose $\underline{\gamma}_{(2)}\leq 1$ and $\underline{\gamma}_{(3)}\leq 2\alpha+1$. Thus, from (6.4), we obtain as desired

[TABLE]

The remaining terms can be bounded similarly. Explicitly, we get the following bounds

[TABLE]

that ensure that the stochastic integrals are martingales, and

[TABLE]

In the previous inequalities, we observe that H4 eliminates all possible negative moments in the $I^{i}_{j}$: for $|I_{1}^{2}|$, H4 imposes $\underline{\gamma}_{(1)}=0$. Similarly, for $|I_{2}^{1}|$ and $|I_{2}^{4}|$ and the definition of $\underline{\beta}$, $\underline{\underline{\beta}}$ in Proposition 5.1, H4 imposes $\underline{\gamma}_{(2)}\leq 1,~{}\underline{\gamma}_{(3)}\leq 2$ and $\underline{\gamma}_{(4)}\leq 4$, respectively. Further, the terms $|I_{2}^{1}|$ and $|I_{2}^{3}|$ contain the highest moments to be controlled, $4\alpha+\overline{\overline{\beta}}$ and $6\alpha+2\overline{\beta}$, both are less than the moment order $6\alpha+2\overline{\overline{\beta}}\vee(\overline{\overline{\beta}}+2\alpha)$ imposed by H5.

Appendix A Proof of Lemmas 5.4 and 5.5

To simplify notation, we omit the dependence on $\lambda$ in the processes $X^{x},X^{x,\epsilon},J^{x},J^{x,\epsilon}$.

Proof of lemma 5.4.

We consider first the process $J_{t}^{x,\epsilon}=\tfrac{X_{t}^{x+\epsilon}-X_{t}^{x}}{\epsilon}$ satisfying the linear SDE (5.13). From Corollary 5.2 with Condition $\Lambda$.2,

[TABLE]

Similarly, using also the $2(\alpha-1)$-locally Lipschitz property of $b$ in H2

[TABLE]

Therefore, $(\int_{0}^{t}\psi^{\epsilon}_{s}d\theta dW_{s};\;0\leq t\leq T)$ is a square integrable martingale, Equation (5.13) admits a unique strong solution given by the following exponential form (see e.g. [25, Thm V.52])

[TABLE]

In turn, $J_{t}^{x,\epsilon}\geq 0$ yields the increasing property of the flow, $X_{t}^{x+\epsilon}\geq X_{t}^{x}$, and the increasing property of the map $\theta\mapsto\left(X_{t}^{x}+\epsilon\theta J_{t}^{x,\epsilon}\right)$ in $[0,1]$, from which we obtain the following relation

[TABLE]

From the exponential form (A.1), for all $q>0,$

[TABLE]

where the last inequality is obtained from H3’: $\xi_{t}^{\epsilon}\leq\;B_{1}^{\prime}-B_{2}^{\prime}\phi_{t}^{\epsilon}$. Since $(\psi_{t}^{\epsilon})^{2}\leq\phi_{t}^{\epsilon}$, by choosing $q>0$ such that $(2\alpha-1)\lambda+\tfrac{\alpha^{2}}{2}(q-1)\leq\tfrac{B_{2}^{\prime}}{\sigma^{2}}$, we get

[TABLE]

Since the process in the right hand-side is a supermartingale, we have for all $t\in[0,T]$,

[TABLE]

and from this, we deduce that \mathbb{E}\big{[}|X_{t}^{x+\epsilon}-X_{t}^{x}|^{q}\big{]}\leq\epsilon^{q}\exp\{qB_{1}^{\prime}T\} tends to [math] with $\epsilon$.

The exponential form (5.7) allows to replicate the computation above for the process $J^{x}$, and conclude similarly that for all $t\in[0,T]$, \mathbb{E}\big{[}(J_{t}^{x})^{q}\big{]}\leq\exp\{qB_{1}^{\prime}T\}.

Proof of Lemma 5.5.

Consider the difference $\mathcal{E}_{t}^{x,\epsilon}:=J_{t}^{x}-J_{t}^{x,\epsilon},$ which can be rewritten, using the SDEs (5.6) and (5.13), as

[TABLE]

Introducing the stopping time $\tau_{M}:=\{0\leq t\leq T:J_{t}^{x}-J_{t}^{x,\epsilon}\geq M\}$, with $M>0$, Itô’s formula yields

[TABLE]

Using (A.2), there exists a non-negative constant $C$ independent on $\epsilon$ and $M$ such that

[TABLE]

and

[TABLE]

Condition ($\Lambda$.3) on $B^{\prime}_{2}$ allows to bound up to $\mathbb{E}[\left(J_{t}^{x,\epsilon}\right)^{6}]$, and on $B^{\lambda}_{2}$ allows to bound up to $\mathbb{E}[|X_{t}^{x}|^{4(\alpha-1)}]$. Similarly, from the $(\overline{\gamma}_{(1)},\underline{\gamma}_{(1)})$-locally Lipschitz property of $b^{\prime}$, there exists a constant $C\geq 0$ independent on $\epsilon$ and $M$ such that

[TABLE]

and $\mathbb{E}[|X_{t}^{x}|^{2\overline{\gamma}_{(1)}}]$ is bounded under ($\Lambda$.3).

In (A), summing separately the three terms multiplying $(\mathcal{E}_{s}^{x,\epsilon})^{2}$, using H3’ and next the Condition $\Lambda$.3 on $B^{\prime}_{2}$ we get

[TABLE]

Coming back to (A) using inequalities (A.4) and (A.5), and $\mathcal{E}^{x,\epsilon}\leq J_{t}^{x,\epsilon}+J_{t}^{x}$,

[TABLE]

and applying Gronwall’s Lemma. From which we obtain $\mathbb{E}(\mathcal{E}_{t\wedge\tau_{M}}^{x,\epsilon})^{2}\leq C~{}{\epsilon}$ for all $t\in[0,T]$. We end this proof by taking limits $M\rightarrow+\infty$ and $\epsilon\rightarrow 0$.

Appendix B

Final step for the proof of Proposition 5.1

We define for some non-negative integers $n,m$, the function $\widetilde{b}_{n,m}$ as

[TABLE]

with, for $j=0,1,2,3$, and using the fact that $2\alpha-1\leq\overline{\gamma}_{(j+1)}+j$ (see H2’),

[TABLE]

with the help of Lemma 1.2, From (5.24) and (5.25), we get

[TABLE]

We identify (B.2) with the form (5.4) with

[TABLE]

Using (B.1) on $\widetilde{b}$, and the control of $\tfrac{\partial^{2}u}{\partial x^{2}}$ in (5.26), we determine the powers involved in the upper bound of $|h^{\prime}_{2}|$ (that will coincide with the moments to bound for the control of $\tfrac{\partial^{3}u}{\partial x^{3}}$) by evaluating

[TABLE]

(using the Hardy symbol $\preceq$ as asymptotic notation) and hence $|h_{2}^{\prime}|(x)\leq(1+x^{\overline{\beta}}+x^{-\underline{\beta}})$, with

[TABLE]

Therefore, we apply Proposition 5.3 with $\rho_{0}=\overline{\beta}$, $\rho_{2}=\rho_{4}=\overline{\gamma}_{(2)}+1$, that must satisfy Condition $\Lambda$.1 for $\lambda=2\alpha$:

[TABLE]

where we write $X^{x}$ and $J^{x}$ for $X^{x}(2\alpha)$ and $J^{x}(2\alpha)$.

Estimates on $\tfrac{\partial^{3}u}{\partial x^{3}}$ and $\tfrac{\partial^{4}u}{\partial x^{4}}$.

We apply the same technique as for the second derivative, namely, we rewrite the second and fourth terms of the sum in (B), using the Markov property and time homogeneity of the process $(X_{s}^{x}(2\alpha);0\leq s\leq T-t)$ for the second term in $f^{(2)}$ in (B):

[TABLE]

We also use an integration by part in the second line of (B):

[TABLE]

where again we write $X^{x}$ and $J^{x}$ for $X^{x}(2\alpha)$ and $J^{x}(2\alpha)$. Then, substituting in (B) we get

[TABLE]

We consider the change of measure $\mathbb{Q}^{3\alpha}$ through the density $\mathcal{Z}_{t}^{(2\alpha,3\alpha)}$ (assuming (5.19) with $\lambda=2\alpha$, for which we observe that

[TABLE]

Therefore, using again that $\textit{Law}^{\mathbb{Q}^{3\alpha}}(X^{x}(2\alpha))=\textit{Law}^{\mathbb{P}}(X^{x}(3\alpha))$, we obtain

[TABLE]

Notice that $\widetilde{b}_{3,3}$ is bounded from above assuming $B^{\prime}_{2}\geq\alpha\sigma^{2}(2\alpha-1)$. By means of the boundedness of $\tfrac{\partial u}{\partial x}$ and $f^{(i)}$, we stay with

[TABLE]

Now using Corollary 5.2, with $\overline{\beta}\leq 1+\tfrac{2B_{2}^{3\alpha}}{\sigma^{2}}$, we get

[TABLE]

In order to apply Proposition 5.3 a last time, we identify in (B) the form (5.4) with

[TABLE]

Again, using (5.21), (5.26) and (B.5), we estimate the powers involved in the expression on $|h^{\prime}_{3}|$ (that will coincide with the moments to bound for the control of $\tfrac{\partial^{4}u}{\partial x^{4}}$) by evaluating

[TABLE]

with

[TABLE]

Then, assuming $\overline{\overline{\beta}}\leq\tfrac{1}{2}+\tfrac{B_{2}^{3\alpha}}{\sigma^{2}}$, we apply Proposition 5.3, obtaining $u\in\mathcal{C}^{1,4}([0,T]\times\mathbb{R}^{+})$. Using the Markov property and the time homogeneity of the process $(X_{s}^{x}(3\alpha);0\leq s\leq T-t)$, we deduce the following form (with $X_{s}^{x}$ understood as $X_{s}^{x}(3\alpha)$)

[TABLE]

Considering the change of measure $\mathbb{Q}^{4\alpha}$ with density $\mathcal{Z}_{t}^{(3\alpha,4\alpha)}$, we have

[TABLE]

with $\widetilde{b}_{4,6}(x)=4b^{\prime}(x)+6\alpha\sigma^{2}(2\alpha-1)x^{2(\alpha-1)}$ bounded from above according to

[TABLE]

and the assumption that $B_{2}^{\prime}>\tfrac{6}{4}\alpha\sigma^{2}(2\alpha-1)$. Therefore we start to bound $|\tfrac{\partial^{4}u}{\partial x^{4}}|$ with

[TABLE]

Combining this with the previous polynomial bounds for the derivatives and the control of moments for the process $X^{x}(4\alpha)$ in Corollary 5.2, under H5, we get

[TABLE]

We have obtained that $u\in\mathcal{C}^{1,4}([0,T]\times[0,+\infty))$ with partial derivatives satisfying (5.4). In view of the polynomial growth property of the maps $x\mapsto\frac{\partial u}{\partial x}(t,x)$, $\tfrac{\partial^{2}u}{\partial x^{2}}(t,x)$, $b(x)$, $x^{\alpha}$ and the appropriate control of the $\overline{\overline{\beta}}$-th moment of the flow, one can easily adapt the proof in Friedman [11, Ch. 5, Th 6.1] to show that $u(t,x)$ solves the Kolmogorov PDE (5.3).

We end this proof by reporting the conditions required on $B_{2},B^{\prime}_{2}$, $\sigma$, $\alpha$, $\overline{\gamma}_{(i)}$, $\underline{\gamma}_{(i)}$ in order to get all the controls to be applied in the previous steps, the combination of which forming H5:

$\bullet$ At most we used the upper-bound on the moment $\sup_{t\in[0,T]}(\mathbb{E}|(X_{t}(4\alpha))^{\overline{\overline{\beta}}}|+\mathbb{E}|(X_{t}(3\alpha))^{2\overline{\overline{\beta}}}|)$, by applying Corollary 5.2 with the double constrain that \tfrac{\sigma^{2}}{2}\big{(}\overline{\overline{\beta}}+8\alpha-1\big{)}\leq B_{2} and \tfrac{\sigma^{2}}{2}\big{(}2\overline{\overline{\beta}}+6\alpha-1\big{)}\leq B_{2}, knowing that $\overline{\overline{\beta}}:=3(\overline{\gamma}_{(2)}+1)\vee(\overline{\gamma}_{(2)}+\overline{\gamma}_{(3)}+2)\vee(\overline{\gamma}_{(4)}+1)$.

$\bullet$ We also have to justify the use of Girsanov transform, by applying Lemma 5.6 under the sufficient condition that

[TABLE]

$\bullet$ We have bounded the terms involving $J_{t}$ coming after the Girsanov transform, and at most the term

[TABLE]

by assuming $B_{2}^{\prime}>\sigma^{2}\alpha(3\alpha-\tfrac{3}{2})$.

$\bullet$ Finally, we considered the necessary condition on $B_{2}^{\prime}$ in order to apply Proposition 5.3 up to $\lambda=3\alpha$: ${B_{2}^{\prime}}\geq\sigma^{2}\alpha(\tfrac{17}{2}\alpha-3)$.

Appendix C Proof of Proposition 2.1 and related Lemmas

Proof of Proposition 2.1.

From the Lamperti-type transformation $\widetilde{X}=X_{t}^{-2(\alpha-1)}$, with $\alpha>1$, the well-posedness of (1.1) essentially relies on the existence and uniqueness of a positive solution to the one-dimensional SDE

[TABLE]

where the drift function $\widetilde{b}$, defined as

[TABLE]

is a $(\overline{\gamma},\underline{\gamma})$-locally Lipschitz, for some strictly positive $\overline{\gamma}$ and $\underline{\gamma}$. Then pathwise uniqueness holds for the solution of (C.1) (see e.g. Ikeda-Watanabe [19, Theorem 3.1]).

Feller test for non explosion.

From the regularity of the coefficients of (C.1), a weak solution up to an explosion time exists (in the sense of Definition 5.1 in Karatzas-Shreve [20]). Considering

[TABLE]

we use a Feller test to show that $\mathbb{P}(S=+\infty)=1$, or equivalently that $v(0^{+}):=\lim_{x\rightarrow 0^{+}}v(x)=v(+\infty):=\lim_{x\rightarrow+\infty}v(x)=+\infty$, where $x\mapsto v(x)$ is defined from the scale function $x\mapsto p(x)$ as

[TABLE]

with $c>0$ (e.g. [20, Theorem 5.29]). From the $2(\alpha-1)$-locally Lipschitz continuity of $b$, we have, for all $x>0$,

[TABLE]

Multiplying both sides of the inequality by $-x^{\tfrac{1}{2(\alpha-1)}}$, it follows that

[TABLE]

and

[TABLE]

In the same way, using H3, we get the lower bounds

[TABLE]

[TABLE]

From the lower and upper bounds on $\frac{\tilde{b}(x)}{x(\alpha-1)}$, we derive the following estimates for $p(0^{+})$ and $p(+\infty)$:

[TABLE]

where $a\mapsto\beta(a)=\tfrac{{2a}/{\sigma^{2}}+2\alpha-1}{2(\alpha-1)}>1$ for $a\geq 0$. Therefore we have that $p(0^{+})=-\infty$ implies $v(0^{+})=+\infty$ (see e.g [20], Problem 5.27).

To check that $v(+\infty)=+\infty$, we distinguish the two cases: $b(0)>0$ and $b(0)=0$ in the upper-bound (C.3). If $b(0)>0$, the dominant term in (C.3) corresponds to the integral of $z\mapsto\exp\{\frac{2b(0)}{\sigma^{2}(2\alpha-1)}z^{\frac{2\alpha-1}{2\alpha-2}}\},$ from which we obtain $p(+\infty)=+\infty$ implying $v(+\infty)=+\infty$. If $b(0)=0$, we must analyze the behavior of $v$ at infinity, for which we have

[TABLE]

where we are allowed to choose $\widetilde{\beta}>\beta(C)+1$ such that $\widetilde{\beta}-2>0$. Applying an integration by parts on the inner integral, we have

[TABLE]

and consequently

[TABLE]

The second integral being finite and the first one diverging to $+\infty$, we can conclude that $v(+\infty)=+\infty$ and thus there exists a unique strictly positive strong solution to (C.1) for all $t\in[0,T]$. Using reversely the Lamperti transformation, this immediately implies that $X_{t}={\widetilde{X}_{t}}^{\frac{1}{2(1-\alpha)}}$ satisfies the SDE (1.1) on $(0,T]$ and the pathwise uniqueness of strictly positive solution is also granted.

Moment controls for the solution to (1.1).

Applying Itô’s formula to the stopped process $({X}_{t\wedge\tau_{M}};0\leq t\leq T)$ with $\tau_{M}=\inf\{t\in[0,T]:{X}_{t}\geq M\},$ and, using H3, we get: for all $p\geq\tfrac{1}{2}$

[TABLE]

Then, for all $p$ such that $1\leq 2p\leq\tfrac{2B_{2}}{\sigma^{2}}+1$, the last term is non-positive and so, for some constant $C$ depending on the parameters $p,B_{i},b(0)$ and $T$, for all $0\leq t\leq T$, we get

[TABLE]

From Gronwall’s inequality, we conclude on the $2p$ th-moment control of $X$. We extend the result for an arbitrary exponent $0\leq 2p\leq 1$, by simply applying Hölder’s inequality.

Proof of Lemma 2.3.

Let $\widetilde{X}$ be the previous Lamperti transform, solution to (C.1). By means of the Itô’s lemma (for simplification, we omit the localization argument previously used in the proof of Proposition 2.1) and the upper bound (C.2), we have for all $t\in[0,T]$, and for all power $q>0$

[TABLE]

where $x\mapsto\varphi(x):=2C+\sigma^{2}+2Cx-2b(0)x^{\frac{2\alpha-1}{2\alpha-2}}$ achieves its maximum at the point $x^{*}=(\tfrac{2C(\alpha-1)}{b(0)(2\alpha-1)})^{2\alpha-2}$. It follows that

[TABLE]

By Gronwall’s lemma, we deduce the finiteness of all positive moments of $\widetilde{X}$, and retrieving $X$ through the transformation $X_{t}={\widetilde{X}_{t}}^{\frac{1}{2(1-\alpha)}}$, we conclude on the lemma.

Proof of Lemma 2.4.

When $b(0)=0$, we consider CIR process $Y$, unique strictly positive strong solution to the SDE

[TABLE]

(the strict positivity is granted as two times the drift at point zero is greater than the square of the diffusion parameter, see e.g. [4]). Applying the classical comparison principle for one-dimensional SDEs (see e.g [19, Chp.VI, Th 1.1]), we have

[TABLE]

The bound (2.2) directly follows from the application of Lemma A.2 in [4] (with $\nu=\frac{\sigma^{2}+2B_{2}}{2\sigma^{2}(\alpha-1)}$).

We consider now the case $b(0)>0$. By Itô’s formula

[TABLE]

Using the growth condition of $b$ in H3, we get

[TABLE]

Hence, applying Hölder’s inequality for $q,p$ such that $1=\frac{1}{q}+\frac{1}{p}+\frac{\mu}{B_{2}\sigma^{2}}$, we have

[TABLE]

(using the introduced notation $\mathbb{E}^{\beta}$). Expanding the last term using series and using Jensen’s inequality, we have for all $t\in[0,T]$ and all $\mu<B_{2}\sigma^{2}$

[TABLE]

where $\theta=\frac{2\mu b(0)p}{\sigma^{2}}$ and the inequality in (C.6) is obtained from the Fubini-Tonelli Theorem.

It remains to prove the finiteness of the series in the inequality above. Yet estimate on the negative moments of $X$ in Lemma 2.3 (or positive moments for $\widetilde{X}$) does not permit to directly conclude as we need to explicit the way the order $k$ appears in the estimation. To overcome this difficulty we come back to estimate (C.4) and balance the dependence on the power as follows: for all $t\in[0,T]$ and all power $q>1$, by Young’s inequality, we have $\widetilde{X}_{t}^{q-1}\leq\frac{1}{q}+\widetilde{X}_{t}^{q}$ and $q\widetilde{X}_{t}^{q-1}\leq q^{q-1}+\widetilde{X}_{t}^{q}$ from which we deduce that for any $a>0$ and $b>0$

[TABLE]

Applying these to the right hand side of (C.4) gives

[TABLE]

where $a_{1}=(\alpha-1)\varphi(x^{*})$, $a_{2}=\sigma^{2}(2\alpha-2)+\varphi(x^{*})$. Next, applying Gronwall’s inequality, we get

[TABLE]

and this bound can be extended for arbitrary $q>0$ using Jensen’s inequality. The estimation of the $(-k)$-th moment for $X$ is then obtained from the transformation $X_{t}^{-k}={\widetilde{X}_{t}}^{\frac{k}{2(\alpha-1)}}$:

[TABLE]

Coming back to (C.6), we use the above estimate as follows

[TABLE]

by setting $\beta(k)=\tfrac{1}{k!}\left(\tfrac{\theta~{}T\exp\{\frac{a_{2}}{2}T\}}{(2\alpha-2)^{1/(2\alpha-2)}}\right)^{k}~{}k^{\frac{k}{2(\alpha-1)}}$. We observe that

[TABLE]

which converges to zero when $\alpha>\frac{3}{2}$, and thus $\sum_{k=0}^{\infty}\beta(k)$ is finite. This ends the proof by substituting (C.8) in (C.6).

Appendix D Some complementary numerical experiments

We complement Section 4 by testing the obtained theoretical convergence rate through a larger set of model cases (4.1), where the parameters of the model satisfy all the assumptions, and some other cases where not all the hypotheses are fulfilled. In particular we illustrate the fact that hypothesis H5 do not correspond to a necessary condition.

We also report on the comparison between exp-ES and other proposed schemes of the literature (see Section D.1).

Numerical parameters.

For all the presented numerical experiments, we consider a unit terminal time $T=1$, the initial condition $x=1$ and the time step $\Delta t=1/2^{p}$, for $p=1,\ldots,10$. In addition, the simulation of the ensemble average of the scheme $\mathbb{E}f(\overline{X}_{T})$ is performed by a Monte Carlo approximation, according to $n=10^{5}$ independent trajectories.

Test functions.

Along this section, we consider four different test functions, not all bounded,

[TABLE]

Model cases.

Denoting $\kappa$ as the left-handside of (4.2) or (4.3), we consider the following cases, determined by the data $(B_{0},B_{1},B_{2},\sigma,\alpha)$

[TABLE]

with three of them, Cases 2, 3 and 6, that are not satisfying H5.

The model in Case 3 satisfies assumptions of Theorem 2.1 in [17] that states that the approximated moments by the Euler-Maruyama scheme and the strong $L^{p}$-error associated to moment-approximations diverges. Also in [17], the authors prove the divergence in the weak sense for the $p$-th moments of the Euler-Maruyama scheme in that case.

Computation of the reference values.

From Proposition 2.1 and Lemma 2.3, it is easy to check that for Case 1 to Case 7, the parameters guarantee the finiteness of the expectation $\mathbb{E}f(X_{T})$ for each test function in the set $\{x,x^{2},\frac{1}{x},\exp(-x^{2})\}$.

For both test functions $f(x)=x$ and $f(x)=x^{2}$, it is possible to compute the reference values of $\mathbb{E}[X_{T}]$ and $\mathbb{E}[X_{T}^{2}]$, for Case 1 to Case 6, through Lamperti transformation leading to CIR process for which a closed form for negative moments is available (see details in [22, Chapter 1]). Thus, the first and the second moment reduces to

[TABLE]

The choices of the power $\alpha$ in the Cases listed above are made in order to correspond to some explicit values of $\Gamma(\frac{1}{2\alpha-2})$ and $\Gamma(\frac{1}{\alpha-1})$.

For the other cases, $f(x)=\frac{1}{x}$ and $f(x)=\exp(-x^{2})$, we compute the reference values based on a Monte Carlo method combined with the scheme exp-ES (3.3)

[TABLE]

with $n_{0}=10^{6}$ or $10^{7}$, depending on the case, and $\Delta t_{\text{ref}}=2^{-14}$.

Numerical results are shown in Table 2. Except for Case 6 in red cells and Case 3 with test function $f(x)=x^{2}$, we can observe the rate of convergence of order one in all the rows corresponding to the selection of bounded/unbounded test functions: the error is divided by 2 when going from left to right, even if some saturation can be observed for the smallest error values $(p=8,9,10)$ when Monte Carlo error starts to be dominant. This behavior is also confirmed in Figure 1, illustrating the results obtained in Table 1 in a log-log scale, with an additional line (in black) representing the reference convergence order 1. This confirms that our proofs can certainly be extended for a larger class of test functions, and model parameters.

In particular, we highlight the parameter Case 3 $(0,0,1,1,3/2)$ that converges weakly with order one for $f(x)=x,\frac{1}{x},\exp\{-x^{2}\}$, even if H5 is not fulfilled, and even moreover we know that the classical Euler-Maruyama scheme is strongly diverging (as stated in [17]) in this case.

D.1 Some comparison with other schemes

One of the significant advantages of the proposed exp-ES is that it easily addresses the control of the moments of the numerical approximation (see Lemma 3.2). We seek now to (numerically) observe the stability of the exponential-Euler scheme. To this aim, we compare the exp-ES with the four others following schemes proposed in the literature:

• •

Symmetrized Euler scheme (SES) (see e.g. [4] and the reference therein), defined by

[TABLE]

which is the closest form of the classical Euler scheme to be applied to SDE (4.1).

• •

Symmetrized Milstein scheme (SMS) (see e.g. [5] and the reference therein), defined by

[TABLE]

• •

Tamed Euler scheme (TES, see [18]), defined by

[TABLE]

• •

Stopped tamed Euler scheme (STES, see [16] and the reference therein) defined by

[TABLE]

Results are shown in Figure 3 and Table 3. Table 3 reports on the stability of the exp-ES, in comparison with the other schemes. In particular, we experiment some instability with the tamed schemes when $\Delta t$ is not small enough (marked in Table 3 as - for the missing values). We also observe abnormally large level of errors for SES (Case 1) and SMS (Cases 1, 3, 5) when $\Delta t$ is not small enough as well.

In terms of convergence rate, the scheme exp-ES behaves very well, in the average of the other schemes, and even better in Cases 2, 3, 5. On the contrary, bf Case 3 (where the explicit Euler scheme is strongly diverging) is particularly unstable for the SMS and TES, STES. The same behavior with a smaller impact is observed in bf Case 5.