The Euler-Maruyama approximation for the absorption time of the CEV diffusion

A standard convergence analysis of the simulation schemes for the hitting times of diffusions typically requires non-degeneracy of their coefficients on the boundary, which excludes the possibility of absorption. In this paper we consider the CEV diffusion from the mathematical finance and show how a weakly consistent approximation for the absorption time can be constructed, using the Euler-Maruyama scheme.


Introduction
In scientific computations, expectations with respect to probabilities, induced by continuous time processes, are often replaced by Monte Carlo averages over independent trajectories. For diffusions, generated by stochastic differential equations (SDEs), the trajectories are usually approximated numerically (see e.g. [16]). The accuracy assessment of such numerical procedures is a well studied topic and the available theory establishes and quantifies the convergence of the approximations to the actual solution in a variety of modes, depending on the properties of the SDE coefficients. This in turn typically suffices to claim convergence of expectations for path functionals, continuous in an appropriate topology, but unfortunately, may not apply to discontinuous functionals, some of which arise naturally in applications.
One important example of such a functional is the hitting time of a domain boundary. Let X = (X t ) t≥0 be the diffusion process on R, generated by the Itô SDE dX t = b(X t )dt + a(X t )dB t , X 0 = x ∈ R, (1.1) recursion (2.2) and (2.3) below) and suppose that X δ approximates the diffusion X in the sense of weak convergence. More precisely, let C([0, ∞), R) be the space of real valued continuous functions on [0, ∞), endowed with the metric Then X δ converges weakly to X, if for any bounded and continuous functional h : Eh(X δ ) = Eh(X).
In particular, if φ : C([0, ∞), R) →R + is a functional, almost surely continuous with respect to the measure induced by X, then (1.3) implies lim δ→0 Ef φ(X δ ) = Ef φ(X) (1.4) for any continuous and bounded function f :R + → R. In other words, the weak convergence of the processes X δ → X implies the weak convergence of the random variables φ(X δ ) → φ(X) or, equivalently, the convergence of the probability distribution functions lim δ→0 P(φ(X δ ) ≤ x) = P(φ(X) ≤ x) for any x ∈R + , at which the distribution function P(φ(X) ≤ x) is continuous.
On the other hand, τ ℓ (·) is continuous at any u, which either upcrosses ℓ: If this type of paths is typical for the diffusion (1.1), i.e. if the set of all such paths is of full measure, induced by the process X, then τ ℓ (·) is essentially continuous and the weak convergence X δ → X still implies the weak convergence τ ℓ (X δ ) → τ ℓ (X). However, if ℓ is an accessible absorbing boundary of X, the paths which hit ℓ cannot leave it at any further time. For such diffusions τ ℓ (·) is discontinuous on a set of positive probability and the weak convergence τ ℓ (X δ ) → τ ℓ (X) cannot be directly deduced from X δ → X.
Hitting times play an important role in applied sciences, such as physics or finance, and since their exact probability distribution cannot be found in a closed form beyond special cases, practical approximations are of considerable interest. There are two principle approaches to compute such approximations. One is based on the fact that the expectation of a given function of the hitting time solves the Dirichlet boundary problem for an appropriate PDE, and thus the approximations can be computed using the generic tools from the PDE numerics. Sometimes, the particular structure of the emerging PDE can be exploited to calculate expectations of special functions of hitting times, such as moments (as e.g. in the linear programming approach of [10]).
The probabilistic alternative is to use the Monte Carlo simulations, in which the diffusion paths are approximated by numerical solvers. Typically the diffusions are simulated on a discrete grid of points and the evaluation of the hitting times requires construction of the continuous paths through an interpolation. The naive approach is to use the general purpose interpolation techniques, such as the one used in our paper (see (2.3) below). Better results are obtained if the possibility of having a hit between the grid points is taken into account as in e.g. [17], [6], [7], [12].
The convergence analysis of the approximations of the hitting times, based on the various numerical schemes, appeared in [19], [21,20], [9], [8]. The results obtained in these works, assume ellipticity or hypoellipticity of the diffusion processes under consideration, which corresponds to the case of non-absorbing boundary in the preceding discussion. The analysis beyond these non-degeneracy conditions appears to be a much harder problem.
In this paper we consider a particular diffusion on R + , with an absorbing boundary at {0}. As explained above, the absorption time in this case is a genuinely discontinuous functional of the diffusion paths, which makes the convergence analysis of the approximations a delicate matter. We propose a simple approximation procedure, based on the Euler-Maruyama scheme, and prove its weak consistency.
In the next section we formulate the precise setting of the problem and state the main result, whose prove is given in Section 3. The results of numerical simulations are gathered in Section 4 and some supplementary calculations are moved to the appendices.

The setting and the main result
Consider the diffusion process X = (X t ) t≥0 , generated by the Itô SDE where µ ∈ R, σ > 0 and p ∈ [1/2, 1) are constants and B = (B t ) t≥0 is the Brownian motion, defined on a filtered probability space (Ω, F, (F t ), P), satisfying the usual conditions. This SDE has the unique strong solution (Proposition 1 in [1]) and is known in mathematical finance as the Constant Elasticity of Variance (CEV) model (see e.g. [3]). For p = 1/2, it is also Feller's branching diffusion, being the weak limit of the Galton -Watson branching processes under appropriate scaling.
The process X is a regular diffusion on R + ∪ {0} and a standard calculation reveals that {0} is an absorbing (or Feller's exit) boundary (see §6, Ch. 15 in [13]). We will denote by (P x ) x∈R + the corresponding family of Markov probabilities with P x (X 0 = x) = 1, induced by X on the measurable space C([0, ∞), R + ) with the metric (1.2).
Consider the continuous time process X δ = (X δ t ), t ≥ 0, which satisfies the Euler-Maruyama recursion at the grid points t j ∈ δZ + and is piecewise linear otherwise: where x + := max(0, x), δ > 0 is a small time step parameter and (ξ j ) j∈Z + is a sequence of i.i.d. N (0, 1) random variables.
Since the diffusion coefficient of (2.1) degenerates and is not Lipschitz on the boundary {0}, this SDE does not quite fit the standard numerical frameworks such as [15] or [22]. Nevertheless the scheme (2.2) does approximate the solution of (2.1) in the sense of the weak convergence of measures, as was recently shown in [1] (see also [27], [26], [2]). Consequently, for any P x -a.s. continuous functional φ : where w → stands for the the weak convergence, defined in (1.4). Since a typical trajectory of X oscillates around the level a > 0 after crossing it, the functional τ a (·) is P x -a.s. continuous for a > 0 (see Lemma 3.3 below) and hence τ a (X δ ) This argument, however, does not apply to τ 0 (·), since it is essentially discontinuous, as discussed in the Introduction. Leaving the question of convergence τ 0 (X δ ) w → τ 0 (X) open, we shall prove the following result Note that τ δ β (·) is a continuous functional for any fixed δ > 0 and hence can be seen as a mollified version of the discontinuous τ 0 (·). The parameter β controls the mollification, relatively to the step-size parameter of the Euler-Maruyama algorithm.
In practical terms, the convergence (2.5) provides theoretical justification for the procedure, in which the approximate trajectory X δ , generated by (2.2) and (2.3), is stopped not at 0, but at δ β > 0, which only approaches zero as δ → 0. Our method does not quantify the convergence in terms of e.g. rates, but the numerical experiments in Section 4 indicate that this stopping rule produces practically adequate results.
Remark 2.2. As will become clear from the proof, our approach exploits the local behavior of the SDE coefficients near the boundary and can be applied to the more general one-dimensional diffusions of the form (1.1), for which a weakly convergent numerical scheme is available. For example, all the arguments in the proof of Theorem 2.1 directly apply to the diffusion, whose coefficients b(x) and a(x) have a similar local asymptotic at 0 as the CEV model. The SDE (2.1) is a study case, which seems to capture the essential difficulties of the problem, related to the degeneracy of the SDE coefficients on the absorbing boundary. It is a convenient choice, since the weak convergence (2.4), being the starting point of our approach, has been already established for the CEV model in [1], and the explicit formulas for the probability density of τ 0 (X) are available, allowing to carry out the numerical experiments.

The proof of Theorem 2.1
The proof, inspired by the approach of S.Ethier [4], is based on the following observation (a variation of Billingsley's lemma, see Proposition 6.1 [4]) Proposition 3.1. Let Q n , n ∈ N and Q be Borel probability measures on a metric space S such that Q n w − → Q as n → ∞. Let S ′ be a separable metric space with metric ρ, and suppose that φ k , k ∈ N and φ are Borel measurable mappings of S into S ′ such that for an increasing real sequence (α n ).
Proof. Let h be a continuous with respect to ρ bounded real valued function on S ′ , then where E and E n denote the expectations with respect to Q and Q n respectively. Since h is continuous and φ k → φ, Q-a.s., the last term vanishes as k → ∞ by the dominated convergence. Moreover, since for any fixed k, φ k is Q-a.s. continuous and Consequently, where the latter equality holds by (iii). The claim follows by arbitrariness of h.
Let us now outline the plan of the proof. In our context, the probability measures, induced by the family X δ , play the role of Q δ and by Proposition 3.2 they converge weakly to the law Q := P x of the diffusion X. Since the diffusion coefficient of (2.1) is positive, away from the boundary point 0, τ ε (·) is a P x -a.s. continuous functional (Lemma 3.3) and hence (i) of Proposition 3.1 holds. On the other hand, lim ε→0 τ ε (u) = τ 0 (u) for any u ∈ C([0, ∞), R) (Lemma 3.4), which implies (ii) of Proposition 3.1. The more intricate part is the convergence (iii), which is verified in Lemma 3.5, using the particular structure of the SDE (2.1). The statement of Theorem 2.1 then follows from Proposition 3.1.
Proof. For the process, obtained by piecewise constant interpolation of the points generated by the recursion (2.2), the claim is established in Theorem 1.1 in [1]. The extension to the piecewise linear interpolation (2.3) is straightforward.
It is left to show that (3.2) holds. The diffusion X satisfies the strong Markov property and thus (we write τ ε for τ ε (X)), and hence the required claim holds, if P ε (τ ε− = 0) = 1.
Take now f to be a scale function of the diffusion X, i.e. a solution to the equation Lf = 0, where L is the generator of X: It is well known, e.g. [14], that we can take it to be positive and increasing, specifically, for (3. 3), 2µz (σz p ) 2 dz dy can be taken. The process f (X t ) is a nonnegative local martingale and thus a supermartingale (e.g. [14] p. 197). The random variable t ∧ τ ε− is a bounded stopping time and by the optional stopping theorem we have for any t ≥ 0 By the definition of τ ε− and path continuity of X, it follows that X t∧τ ε− ≥ ε, and since f is increasing, we have that f (X t∧τ ε− ) ≥ f (ε). Thus it follows from (3.4) that P ε (X t∧τ ε− = ε) = 1 for all t ≥ 0 and, consequently, [X, X] t∧τ ε− = 0, P ε -a.s. for t ≥ 0. On the other hand, [X, X] τ ε− = τ ε− 0 σ 2 X 2p s ds > 0, on the set τ ε− > 0 , P ε -a.s. The obtained contradiction implies P ε (τ ε− > 0) = 0, as claimed.
Let P be the probability on the space, carrying the sequence (ξ j ) (see (2.2)) and denote by (P x ) x≥0 the Markov family of probabilities, corresponding to the discrete time process (X δ t j ) with P x (X δ 0 = x) = 1. Since the process (X δ t ) is piecewise linear off the grid δZ + , the condition (iii) of Proposition 3.1 follows from Lemma 3.5. For any β ∈ 0, 1/2 1−p and η > 0 Proof. Roughly speaking, (3.5) means that a trajectory of X δ , which approaches the boundary {0}, is very likely to hit it. This seemingly plausible statement is not at all obvious, since the coefficients of our diffusion decrease to zero near this boundary, making it hard to reach. By letting the level δ β decrease to zero at a particular rate allows to approximate expectations of the hitting times τ δ β (X δ ) by those of τ δ β (X), which in turn can be estimated using their relations to the corresponding boundary value problems.
In what follows, C, C 1 ,C 2 etc. denote unspecified constants, independent of ε and δ, which may be different in each appearance. Define the crossing times of level a with inf{∅} = ∞. The sequence (X δ t j ), j ∈ Z + is a strong Markov process and ν a /δ is a stopping time with respect to its natural filtration.
Since X δ is piecewise linear, |τ ε (X δ ) − ν ε | ≤ δ and |τ δ β (X δ ) − ν δ β | ≤ δ on the set τ ε (X δ ) < ∞ and thus by the triangle inequality where η ′ := tan(η/2) (assuming δ is small enough). Further, Proof of (3.6). We shall use the regularity properties of the function near the boundary point 0, summarized in the Appendix A. In particular, ψ is continuous on the interval [0, 1] and is smooth on (0, 1]. Note, however, that the derivatives of ψ explode at the boundary point 0, which is related to the possibility of absorption. We shall extend the domain of ψ to the whole R by continuity, setting Consider the Taylor expansion whereX δ t j−1 is between X δ t j−1 and X δ t j . After rearranging terms, the latter reads where L is the generator defined in (3.3), the second term is the martingale and the last term is the residual Consequently, for an integer k ≥ 1 and the stopping time ν := ν δ β ∧ ν 1 , where r(δ) := ψ X δ ν∧kδ − ψ X δ (ν∧kδ)−δ − M ν/δ∧k − M (ν/δ∧k)−1 , is the residual term, which accommodates the possible overshoot at the terminal crossing time ν.
Remark 3.6. The condition β < 1/2 1−p originates in the estimate (3.14), which is plugged into (3.12). The principle difficulty is that the term X δ cannot be effectively controlled for greater values of β as δ → 0. For example, it is not clear how to bound the right hand side of (3.12) already for β = 1 and p = 1/2.

Numerical experiments
In insurance, one is often interested in calculating the probability of ruin by a particular time t > 0. For the diffusion (2.1), this probability can be found explicitly and for p = 1/2, it has a particularly simple form: Note that τ 0 (X) has an atom at {+∞} when µ > 0: are plotted versus δ (in the log scale), along with the 99% confidence intervals, based on the CLT approximation. The results appear to be practically adequate: for example, the accuracy of 0.1% is obtained already with δ = 10 −3 . The positive bias of the error is not surprising, since the earlier absorption is more probable for the larger threshold δ β . The simulation results also indicate in favor of the convergence which remains a plausible conjecture.
Appendix B. An inequality Lemma B.1. Let (η j ) j∈N be a sequence of random variables and N be an integer valued random variable. Then for constants α > 0 and c > 0 and an integer k E|η j |1 {|η j |>cj α } .
Proof. For a fixed integer k ≥ 1 and a real number α > 0 E|η j |1 {|η j |>cj α } Corollary B.2. Let (ξ j ) j∈N be an i.i.d. N (0, 1) sequence and N be an integer valued random variable. Then for any α > 0, p > 0 and integer k with a constant C α,p depending only on α and p. Moreover, for any r ≥ 0, the sum S n = n j=0 |ξ j | p 1 {|ξ j |≥r} satisfies with a constant C p , depending only on p.