The Euler-Maruyama approximations for the CEV model

The CEV model is given by the stochastic differential equation $X_t=X_0+\int_0^t\mu X_sds+\int_0^t\sigma (X^+_s)^pdW_s$, $\frac{1}{2}\le p<1$. It features a non-Lipschitz diffusion coefficient and gets absorbed at zero with a positive probability. We show the weak convergence of Euler-Maruyama approximations $X_t^n$ to the process $X_t$, $0\le t\le T$, in the Skorokhod metric. We give a new approximation by continuous processes which allows to relax some technical conditions in the proof of weak convergence in \cite{HZa} done in terms of discrete time martingale problem. We calculate ruin probabilities as an example of such approximation. We establish that the ruin probability evaluated by simulations is not guaranteed to converge to the theoretical one, because the point zero is a discontinuity point of the limiting distribution. To establish such convergence we use the Levy metric, and also confirm the convergence numerically. Although the result is given for the specific model, our method works in a more general case of non-Lipschitz diffusion with absorbtion.

The aim of this paper is to prove validity of approximations of expectations of some functionals of these diffusions by the Euler-Maruyama scheme. For example, in finance this model represents a price process and it is important to evaluate expectations of payoffs depending on past prices EG(X [0,T ] ); as well as for evaluation of the ruin probability P(τ ≤ T ) by simulations. In population modeling the diffusion with p = 1 2 represents the size of a population and is known as Feller's branching diffusion. In this case τ represents the time to extinction of the population. When exact theoretical expressions of such functionals exist then a comparison of simulated and exact results show how good the approximation is. However, in many cases exact expressions are not available. Then one needs to justify such approximations. Weak convergence established below assures convergence of expectations of bounded and continuous functionals of simulated values to the exact ones.
To allow a larger family of functionals to be approximated we show the weak convergence of approximations in the Skorokhod metric rather than uniform. We also consider approximation of the ruin probability P(τ ≤ t) which is the expectation of the past dependent and discontinuous (in this metric) functional I {τ ≤T } .
The fact that the diffusion coefficient is singular and non-Lipschitz makes the analysis non-standard. Feller in [9] showed for the case p = 1/2 that the solution of the Fokker-Plank equation exists and is unique and gave its fundamental solution. This fact is used for evaluation of ruin probability by assuming that for any T , X T has a density on (0, ∞), and in particular its distribution function is continuous at any positive point, moreover it is right-continuous (as a distribution function), lim ε→0 P(X T ≤ ε) = P(X T = 0).
Previous various results on the Euler-Maruyama algorithm for different type of diffusion models can be found in Kloeden and Platten [20] and Milstein and Tretyakov [22]. They are concerned mostly with the approximation of the final value X T , rather than the whole trajectory X t , 0 ≤ t ≤ T . When both the drift and diffusion coefficients are Lipschitz and the diffusion coefficient is nonsingular, the standard theory applies. There is a body of literature on the topic of approximations, see e.g. Bally and Talay [2], Bossy and Diop [4], Gyöngy [10], Gyöngy and Krylov [11], Halidias and Kloeden [12], Higham, Mao and Stuart [15], Hutzenthaler and Jentzen [16], Jentzen, Kloeden and Neuenkirch [17]. In contrast to the above mentioned papers we always use continuous time setting, even for discrete time models as in [27] since the use of the Skorokhod metric allows to avoid some technical assumptions. Zähle [27] proved the weak convergence in terms of solutions to martingale problems in one dimension. He has shown that the solution of the martingale problem corresponding to (1) can be approximated by solutions of the corresponding time-discrete martingale problems under some conditions. We propose a new method of proof by adapting (with a help of a specific approximation) a standard technique for weak convergence of semimartingales (e.g. [21]) that applies to cases of Hölder continuous diffusion and also when the the limiting process can be absorbed at zero. Our approach avoids the technical condition (2.3) in [27], moreover, the suggested approximations can be used in high dimensions as well.
As a rule, the cases considered in the literature are such that the drift and diffusion coefficients exclude the absorbtion effect even when the diffusion coefficient is singular, such as in Bessel diffusion.
Our method applies to a more general situation of non-Lipschitz diffusion with absorbtion. We chose, however, to give the results for the particular model of CEV for the sake of transparency, and because this model is of significance in applications. The reader will note that the proofs are sufficiently involved already in this special case.
Recall next the Euler-Maruyama scheme for the diffusion process X t . Taking for simplicity equidistant partitions of [0, T ], 0 ≡ t n 0 < t n 1 < . . . < t n tn−1 < t n n ≡ T, t n k − t n k−1 ≡ T n , it is defined by the following recursion where (ξ k ) k≥1 is an i.i.d. sequence of (0, 1)-Gaussian random variables.
The process X n t n k can be considered as a discrete time semimartingale (see, e.g., Zähle [27]) and X n t as a continuous time semimartingale with discontinuous paths. The paths of the process X n = (X n t ) t∈[0,T ] are right continuous piece-wise constant functions with left limits belonging to the Skorokhod space D = D [0,T ] . Consider a metric space (D, d 0 ) endowed with the Skorokhod metric d 0 : if x x x, y y y ∈ D, then where Φ is a set of strictly increasing continuous functions ϕ = (ϕ(t)) 0≤t≤T with ϕ(0) = 0, ϕ(T ) = 1. Denote by Q n the probability measure on the space (D, D), where D is the Borel σ-algebra, of the distribution of X n = (X n t ) t∈[0,T ] , and by Q the distribution of X = (X t ) t∈[0,T ] . Since X is a continuous process, Q(C) = 1, where C denotes the subspace of continuous functions.
Evaluations of functionals by simulations is justified by the convergence of measures in the Skorokhod space for any bounded and continuous in the metric d 0 function f (x x x). The aim of the paper is precisely to show this convergence. Note that such a function f (x x x) need not be continuous in the uniform metric (x x x, y y y) = sup We remark on the use of the function x + next. Since the process X = (X t ) t∈[0,T ] is nonnegative, the notation X + t in (1) is used formally. However, it must be used for the approximating process X n+ t in (2), because X n t can (and does) become negative. Since X n t has piece-wise constant paths the stopping time much more likely corresponds to the negative value rather than the zero value of X n τ n . The plus notation X n+ t in (2) enables us to exclude negative values of X n t up to the first time it becomes negative. Once the process becomes negative it stays at that value for all further times t. Thus, in simulations of f (X) we must use f (X n+ ) and stop at the first time it becomes negative. Next, the function g(x x x) := f (x x x + ) inherits the properties of f (x x x) and is bounded and continuous in the metric d 0 , moreover, when applied to simulations it converges to the desired limit We introduce now the continuous approximation X n t used in the proof of weak convergence The above remark applies equally to the process X n t . Denote Q n the distribution of ( Theorem 1.1. For any T > 0, the Euler-Maruyama approximation for equation (2) converges weakly in the Skorokhod metric d 0 to the limit process The proof is done in three steps: first we show that the distance in the uniform metric between the Euler-Maruyama and the continuous semimartingale approximations converges to zero in probability, secondly we show that the continuous approximation converges in the Skorokhod metric to the solution (X t ) t∈[0,T ] , and finally we deduce convergence of the Euler-Maruyama approximation. Thus the three steps in the proof are The proof of "step 1" uses standard stochastic calculus for semimartingales. The proof of "step 2" follows from the results on diffusion approximation for semimartingales (see [21,Ch.8]). The proof of "step 3" is done as a standard application of We comment on the rate of weak convergence and why defer this important matter to another study. Our focus lies in computer simulation and Monte-Carlo technique. In computer realizations all numerical trajectories are independent objects. Therefore in practical Monte -Carlo simulations these trajectories can not be put on the same probability space. Stronger modes of convergence analyze the rate of decay to zero of the difference X n − X between the pre-limit process X n and the limit process X in various norms as a function of n → ∞ (see e.g. recent papers [17], [14], [26]). The rate of convergence in a weak mode, in particular in the Skorokhod metric, could not be analyzed in terms of the difference X n − X since random processes are defined on different probability spaces. What we have is in fact convergence of pre-limit measures Q n to the limit measure Q, defined on the Skorokhod space, being distributions of X n and X respectively. It is, of course, possible to compare Q and Q n in the total variation Q−Q n tv or Prokhorov's or Wasserstein's metrics. The choice of any of these metrics requires a deep analysis (see e.g. [18] for the total variation metric). There are known results on approximation of Ef (X n T ) at the final point T with a rate of convergence of the |Ef (X T ) − Ef (X n T )| termed as the weak error (see reviews of Kloeden and Platen [20] and Talay [25]). These rates depend on the smoothness of f and regularity of diffusion coefficients. Even for the case of weak convergence at the final point the analysis of the rate is involved, for example Bally and Talay [2] use the Malliavin calculus (stochastic calculus of variation) to establish convergence rate of the distribution function. A recent paper of Alfonsi [1] is devoted entirely to weak error analysis for various approximation schemes for f ∈ C ∞ K (compactly supported infinitely differentiable functions on D), satisfying various conditions which applies to CIR model. Therefore due to complexity of the subject of the rate of convergence, and to preserve the clarity of ideas in the present paper, we postpone the study of this question.
The paper is organized as follows. In Section 2 we give preliminary results used in the proof, with some being of independent interest. The proof of the Theorem 1.1 is given in Section 3. Section 4 gives simulations for the ruin probability P(τ ≤ T ).
To be self contained, the existence and uniqueness of solutions of the stochastic differential equation for the CEV model (1) The process M t is a square integrable martingale for a suitable filtration. So, the Doob maximal Next we give the result useful in the analysis of products of random variables.
Lemma 2.1. Let X n , Y n be sequences of random variables such that X n converges to zero in probability and expectations of Y n are uniformly bounded, sup n E|Y n | = r < ∞. Then X n Y n converges to zero in probability. Proof.
By the Chebyshev inequality P |Y n | > C ≤ r C . Hence Next we need the following result of convergence to zero the expectation of double-maximum the Brownian motion increments over partitions.  Proof. We calculate k-th moment of M n , k ≥ 1. Since M n ≥ 0 The inequality is obtained since (W * i ) 1≤i≤n are i.i.d. random variables, therefore for any x ≥ 0 where we have used the well-known law of the maximum of Brownian motion. Thus where C k depends only on k, (C k = 2 k+2 E|ξ| k with ξ ∼ N (0, 1)). Hence M n converges to zero in L k for any k > 2. But since convergence in L p for a p > 1 implies convergence in L k for any k ∈ [1, p], the statement is proved.
The next elementary inequality is new and is instrumental in the proof.
For x = y = 0 it is obvious.
Consider x > 0 and y < 0. Then, taking into account the proved inequality |x + − y + | ≤ |x − y| (the statement of this lemma for p = 1 2 ) and the fact that y + = 0, we obtain Clearly, the inequality remains true for x < 0 and y > 0.
If both x and y are positive and x > y, then and, in turn, |x p − y p | ≤ |x − y| p . It is easy to see the inequality remains true for y > x.
3. Proof of Theorem 1.1. In the following result we show that the maximum of discrete approximations have uniformly bounded seconds moments. Henceforth r denotes a generic positive constant independent of n with different values at different appearances.
Lemma 3.1. Let X n t n k be the Euler-Maruyama approximation defined in (2).
Proof. First we bound from above max We use repeatedly the following bound of (X n+ Next, for sufficiently large n, there exists r such that Hence for sufficiently large n, we obtain from (6) by using (8) the recurrent inequality Iterating it, for k ≤ n we obtain Thus max 1≤k≤n E|X n t n k | 2 ≤ r.
From the definition of the scheme (2) we obtain by iterations Next we use the Cauchy-Schwartz inequality ( The Cauchy-Schwartz inequality n j=1 |X n t n j−1 | 2 ≤ n n j=1 (X n t n j−1 ) 2 coupled with the bound on the largest second moment shown above (9) implies that the second term in (10) is bounded by a constant r independent of n. The last term is bounded by using the Doob's maximal inequality (4) and the bound in (7) E max Using (9) we can now see that the last term in (10) is bounded by r.
We proceed now to prove the steps in (3).

3.1.
Step 1. The distance between the two approximations in the sup norm converges to zero in probability and in L 2 .
Proof. From the definitions of X n t and X n t it follows that at the points of the partitions both approximations coincide and at intermediate points the following holds By the formula (11), Consequently, (X n , X n ) = sup The first term converges to zero in L 2 using Lemma 3.1 as E max 1≤k≤n |X n t n k−1 | ≤ r. For the second term use Hölder's inequality with parameters 2 p and 2 2−p :
Proof. The proof rests on a general result on the weak convergence of semimartingales to a diffusion [21], Theorem 1, Ch. 8, §3. This theorem states that for weak convergence to a diffusion it is enough to check convergence of the drifts and quadratic variations evaluated at the pre-limit processes. The processes X and X n are semimartingales with following decompositions where we have denoted above (B) B t (X) and B n t ( X n ) are drifts (M) M t (X) and M n t ( X n ) are continuous martingales with predictable quadratic variations The above mentioned Theorem 1 of [21], adapted to the present setting, states that the weak convergence takes place if the following three conditions hold.
(a) Equation (1) has a unique (at least weak) solution Proof. The result follows from the triangular inequality Since convergence in the uniform metric implies convergence in the Skorokhod metric, by Lemma 3.2 we have, in view of f is continuous in this metric, Taking now the lim sup in (13) and using the previous result of weak convergence of X n to X proves the step 3.
4. Evaluation of ruin probability by simulations. In this section we evaluate numerically a ruin probability P(τ ≤ T ), for a finite positive T , where τ = inf{t : X t = 0} by Euler-Maruyama approximations. The basis for analysis is an obvious formula P(τ ≤ T ) = P(X T = 0). It allows us to deal with the distribution function F (x) := P(X T ≤ x) of X T instead of a harder to compute distribution function of τ , P(τ ≤ T ). Notice that In view of this uncertainty, we give approximations for lower and upper bounds of F (0) by using the Lévy metric (see e.g. [13]): It is known that weak convergence of distributions implies convergence in the Lévy Although lim n→∞ L(F n , F ) = 0 does not catch the size of the atom it helps to localize its size. Namely we take x = 0 and a small suitable ε n , in each case determined experimentally, such that the interval [F n (−ε n ) − ε n , F n (ε n ) + ε n ] is small enough and declare that an estimate of F (0) belongs to this interval. In practice, sometimes theoretical values might lie outside this interval. Such values of ε n are pointed out below in each Table below. It is important to compare a numerical method with theoretical results. Fortunately such results are available for p = 1 2 . In this case the explicit formula for P X0 (X T = 0) for any µ, σ is known (see e.g. Dawson [6]) −1] , µ = 0. Our numerical results are obtained for the following values of parameters p, X 0 , µ, σ, T : -p = 1/2, 3/4 -X 0 = 1/10, 1/4, 1 -µ = −1, 0, 1 -T = 3, 6, 9 -σ = 1. We use Monte-Carlo simulations of 10 3 independent copies of the process (X n t ) t∈[0,T ] generated by Euler-Maruyama algorithm with t n k − t n k−1 ≡ 10 −2 . Numerical results are given below for p = 1 2 and p = 3 4 . For p = 1 2 , where theoretical formulae are available, they show good agreement for ε n = 10 −7 , 10 −6 and X 0 = 1 4 and σ 2 = 1. Here "lower bound" and "upper bound" denote respectively P(X n T ≤ −ε) − ε and P(X n T ≤ ε) + ε.  Table 1.  Table 2.
A heuristic explanation why a larger p leads to smaller ruin probabilities follows from the fact that the diffusion parameter x p for x ∈ (0, 1) decreases in p. The same argument clarifies why smaller ε n is needed to achieve the desired accuracy for larger p.

5.
Existence and Uniqueness of solution in the CEV model. The proposition below is given only for reader convenience (see Deelstra and Delbaen [7] for more details). Proof. Define a sequence of processes indexed by integers i, i > 1/X 0 , (X i t ) t≥0 , such that X i t is a strong solutions to the following stochastic differential equation The diffusion coefficient σ(i −1 ∨ |x|) p is Lipschitz continuous therefore X i t is the unique strong solution of (14). Set ϑ i = inf{t : X i t = i −1 }. Note that for t ≤ ϑ i X i+1 t = X i t , and it follows that ϑ i+1 > ϑ i . A strong solution to (1) is constructed by a natural prolongation Finally, Yamada-Watanabe's theorem (see, e.g. [23], p. 265) guaranties the uniqueness of the strong solution of equation (1), because the Hölder parameter p ≥ 1 2 .