Exponential integrability properties of Euler discretization schemes for the Cox-Ingersoll-Ross process

We analyze exponential integrability properties of the Cox-Ingersoll-Ross (CIR) process and its Euler discretizations with various types of truncation and reflection at 0. These properties play a key role in establishing the finiteness of moments and the strong convergence of numerical approximations for a class of stochastic differential equations arising in finance. We prove that both implicit and explicit Euler-Maruyama discretizations for the CIR process preserve the exponential integrability of the exact solution for a wide range of parameters, and find lower bounds on the explosion time.


Introduction
The Cox-Ingersoll-Ross process was originally proposed in Cox et al. (1985) for short-term interest rate modeling, and is the solution to the following stochastic differential equation (SDE): dy t = k y (θ y − y t )dt + ξ y √ y t dW t , (1.1) where W = (W t ) t≥0 is a one-dimensional Brownian motion, whereas y 0 , k y , θ y and ξ y are strictly positive real numbers. According to Karatzas and Shreve (1991), (1.1) admits a unique strong solution, which is strictly positive when the Feller condition is satisfied, i.e., when 2k y θ y > ξ 2 y . The desirable features of the CIR process under consideration, such as non-negativity and mean-reversion, make it very popular when modeling interest rates or variances, e.g., in Heston's stochastic volatility model (Heston 1993). The Feller condition is typically satisfied in practice in the former case, but often fails to hold in the latter case.
The conditional distribution of the CIR process is noncentral chi-squared and hence its increments can be simulated exactly. On the other hand, discretization schemes are usually preferred when the entire sample path of the CIR process has to be simulated, or when the process is part of a system of SDEs. For instance, when pricing path-dependent financial derivatives written on an underlying process S = (S t ) t∈[0,T ] modeled by a d-dimensional SDE, with CIR dynamics in one or more dimensions, we need to evaluate where f : ([0, T ], R d ) → R is the discounted payoff. In particular, this class of SDEs contains the popular Heston model and extensions thereof, such as stochastic interest rates (Grzelak and Oosterlee 2011;Ahlip and Rutkowski 2013) or stochastic-local volatility (van der Stoep et al. 2014). However, one can rarely find an explicit formula for the quantity in (1.2), in which case we approximate the solution to the SDE via a discretization scheme and employ Monte Carlo simulation methods (see Glasserman 2003). Since we cannot use the standard Euler-Maruyama scheme to approximate (y t ) t∈[0,T ] defined in (1.1) because of the non-zero probability of the approximation process becoming negative, we set it equal to zero when it turns negative (absorption fix) or reflect it in the origin (reflection fix). An overview of the explicit Euler schemes considered thus far in the literature can be found in Lord et al. (2010). Alternatively, we can use an implicit scheme to discretize the CIR process. Although weak convergence is important when estimating expectations of payoffs, strong convergence may be required for complex path-dependent derivatives and plays a crucial role in multilevel Monte Carlo methods (Giles 2008). An important step in deriving strong convergence is proving the finiteness of moments of order higher than one of the process and its approximation (Higham et al. 2002;Cozma and Reisinger 2015). In addition, for a number of stochastic volatility models, moments of order higher than one can explode in finite time (Andersen and Piterbarg 2007). This, however, can cause serious problems in practice when valuing securities whose payoffs have super-linear growth, as is the case with some commonly traded fixed income contracts. For instance, the risk-neutral valuation of CMS swaps and caps or Eurodollar futures contracts involves the evaluation of the second moment (Andersen and Piterbarg 2007). Therefore, moment explosions may lead to infinite prices of derivatives. The same issue can also be observed for Euler approximations of SDEs with super-linearly growing drift or diffusion coefficients (Hutzenthaler et al. 2011;Hutzenthaler and Jentzen 2015).
Hence, we need to examine the stability of moments of the actual and the approximated processes, a problem directly related to the exponential integrability of the CIR process and its discretization (Cozma and Reisinger 2015). Hutzenthaler et al. (2014) proved that a class of stopped increment-tamed Euler approximations for nonlinear systems of SDEs with locally Lipschitz drift and diffusion coefficients retain the exponential integrability of the exact solution under some mild assumptions. However, the diffusion coefficient in (1.1) is not locally Lipschitz, so their analysis does not apply to the present work. Cozma and Reisinger (2015) derived the exponential integrability of full truncation Euler approximations (Lord et al. 2010) for the CIR process up to a critical time. In the present work, we first extend the aforementioned result to more general exponential functionals of the CIR process, and then, we prove that the drift-implicit and a number of explicit Euler discretizations for the CIR process preserve exponential integrability properties.
The paper is structured as follows. In Section 2, we discuss the discretization schemes and their strong convergence. In Section 3, we deduce the uniform exponential integrability of functionals of the CIR process and its explicit and implicit Euler discretizations. Section 4 examines moment stability for a particular model (the Heston model) in more detail. Finally, Section 5 summarizes the results and outlines possible future work.

The discretization schemes
The classical Euler-Maruyama scheme does not preserve the non-negativity of the process and hence it is not well-defined when applied to (1.1) directly due to the square-root diffusion coefficient. A number of corrections have been proposed in the literature, by either setting the process equal to zero when it turns negative, or by reflecting it in the origin. Consider a uniform grid: δt = T/N , t n = nδt, ∀n ∈ {0, 1, ..., N }. The partial truncation Euler (PTE) schemeỹ where y + = max (0, y) and δW tn = W t n+1 − W tn , was proposed in Deelstra and Delbaen (1998), whereas the full truncation Euler (FTE) schemẽ y t n+1 =ỹ tn + k y (θ y −ỹ + tn )δt + ξ y ỹ + tn δW tn (2.2) was studied in Lord et al. (2010). The absorption (ABS) scheme reads as For the schemes (2.1) -(2.3), the piecewise constant time-continuous interpolation is defined as Y t =ỹ + tn , whenever t ∈ [t n , t n+1 ). The reflection (REF) schemẽ was introduced in Higham and Mao (2005), and we define Y t = |ỹ tn |, whenever t ∈ [t n , t n+1 ). The symmetrized Euler (SYM) schemẽ y t n+1 = ỹ tn + k y (θ y −ỹ tn )δt + ξ y ỹ tn δW tn (2.5) was studied in Bossy and Diop (2007), and we let Y t =ỹ tn , whenever t ∈ [t n , t n+1 ). Finally, assuming that the Feller condition holds and applying Itô's formula to x t = √ y t leads to where α = 4k y θ y − ξ 2 y 8 , β = − k y 2 and γ = ξ y 2 . (2.7) The drift-implicit (square-root) Euler schemẽ was proposed in Alfonsi (2005) and later studied in Dereich et al. (2012). Because α, γ > 0 and β < 0, (2.8) has the unique positive solutioñ This method is also called the backward Euler-Maruyama (BEM) scheme (Neuenkirch and Szpruch 2014). We let Y t =x 2 tn , whenever t ∈ [t n , t n+1 ). The classical convergence theory (Kloeden and Platen 1999;Higham et al. 2002) does not apply to the CIR process because the square-root diffusion coefficient is not Lipschitz. Consequently, alternative approaches have been employed by the authors to prove the strong or weak convergence of their particular discretization. Strong convergence, either without a rate or with a logarithmic rate, of the partial truncation, the full truncation, the reflection and the symmetrized Euler schemes was established in Deelstra and Delbaen (1998), Alfonsi (2005), Higham and Mao (2005) and Lord et al. (2010). Strong convergence of order 1/2 of the symmetrized and the drift-implicit Euler schemes was proved in Berkaoui et al. (2008) and Dereich et al. (2012), respectively, albeit under some very restrictive assumptions for the former. Recently, Alfonsi (2013) and Neuenkirch and Szpruch (2014) improved the rate of strong convergence of the drift-implicit Euler scheme to 1. To the best of our knowledge, convergence properties of the absorption scheme are not treated in the literature.

Exponential integrability
The goal of this paper is to establish the exponential integrability of functionals of the CIR process and its approximations. To this end, let (Y t ) t∈[0,T ] be the piecewise constant timecontinuous approximation of (y t ) t∈[0,T ] corresponding to one of the discretization schemes from (2.1) -(2.5) or (2.9) and, for some λ, µ ∈ R, define Assuming that t ∈ [t n , t n+1 ] and conditioning on the σ-algebra G tn , we find that Since Y is always non-negative, if ∆ ≤ 0, then E tn Θ t ≤ E tn Θ tn . From the law of iterated expectations, which concludes the proof.

The Cox-Ingersoll-Ross process
The exponential integrability result below is an extension of Proposition 3.1 in Andersen and Piterbarg (2007) to the case ∆ ≤ 0. When ∆ > 0, we can derive Proposition 3.2 directly from Andersen and Piterbarg (2007), by writing the moments of the Heston model in terms of E Θ T . Our proof takes a fairly different approach and is included for completeness.
Proposition 3.2. The first moment of the exponential functional of the CIR process defined in (3.1) is uniformly bounded up to a critical time T * , i.e.,

8)
and Since the first exponential on the right-hand side of (3.15) is a continuous function of time and since (y t ) t∈[0,T ] is a stationary Markov process, it suffices to prove the finiteness of the supremum over τ of is a martingale. Applying Itô's formula to the right-hand side and setting the resulting drift term equal to zero, we find a PDE for F (τ, y): Suppose that the solution to the PDE is of the form The conditions under which G(τ ) -and so F (τ, y) -blows up are connected to the position of the roots of the polynomial f (x) = ax 2 + bx + c, i.e., First of all, if ∆ = 0, then zero is a root of the polynomial. Since f (x) is locally Lipschitz, we conclude that the ODE (3.23) with the initial condition G(0) = 0 has a unique solution, namely G(τ ) = 0, ∀τ ≥ 0. From (3.24), we find that H(τ ) = 0, ∀τ ≥ 0. Substituting back into (3.21) and making use of (3.18), we deduce that E Θ t = 1, ∀t ≥ 0. Alternatively, note that (Θ t ) t∈[0,T ] is a Doléans exponential and a true martingale (see, for instance, Cheridito et al. 2007). Therefore, when (2011), one can easily show the finiteness of G(τ ), and hence of H(τ ) and F (τ, y), for all τ < T * , and the explosion of G(τ ), for all τ ≥ T * . Therefore, F (τ, y) is continuous and finite on [0, T ], for all T < T * , so its supremum over the time interval is finite by the boundedness theorem, which concludes the proof.

The drift-implicit Euler (BEM) scheme
Suppose that 2k y θ y > ξ 2 y and let Y t =x 2 tn , ∀t ∈ [t n , t n+1 ), withx defined in (2.9). In order to derive the exponential integrability of the BEM scheme, we first prove an auxiliary result.
Lemma 3.3. Suppose that ∆ > 0. Then we can find η > 1 such that for all ω ∈ [0, 1], if and only if T < T * , where T * is given below: 1. When µ < 0 and λ < 3 2 µ 2 , (3.28) 2. When µ < 0 and λ ≥ 3 2 µ 2 , or when µ ≥ 0, Proof. Fix any η > 1 and define the polynomial with two distinct real roots However, (3.32) holds for some η > 1 if and only if the second-order polynomial in η on the left-hand side has two distinct real roots, one of which is greater than one. Substituting back into (3.32) with γ = 0.5ξ y from (2.7), we find the necessary and sufficient conditions: Some straightforward calculations lead to an equivalent set of conditions: Henceforth, the conclusion follows relatively easily.

Explicit Euler schemes with absorption fixes
First of all, consider the FTE scheme and let Y t =ỹ + tn , ∀t ∈ [t n , t n+1 ), withỹ from (2.2).
The following result is an extension of Proposition 3.3 in Cozma and Reisinger (2015).
Proposition 3.7. If ∆ ≤ 0 and T ≥ 0 or otherwise, if ∆ > 0 and T ≤ T * , with T * from (3.62) -(3.63), then there exists δ T > 0 such that for all δt ∈ (0, δ T ), the first moments of the exponential functionals from (3.2) of the partial truncation and absorption schemes are uniformly bounded, i.e., sup Proof. We follow the argument of Proposition 3.6 closely and note that (3.74) holds for the partial truncation scheme if δ T ≤ k −1 y , while for the absorption scheme we have equality.

Moment stability in the Heston model
In this section, we consider a filtered probability space (Ω, F, {F t } t≥0 , Q) and suppose that the dynamics of the underlying process are governed by the Heston model under measure Q: (4.1) where W s and W v are correlated Brownian motions with constant correlation ρ ∈ (−1, 1), r is an arbitrary non-negative number and k, θ, ξ > 0 as before. We decompose W s and write it as a linear combination of independent Brownian motionsW s and W v . An application of Itô's formula leads to (4.2) In particular, we are interested in the evaluation of moments E [S ω T ] for ω > 1. Andersen and Piterbarg (2007) give several examples of fixed income securities with super-linear payoffs, whose risk-neutral valuation involves the calculation of the second moment. Hence, moment explosions may lead to infinite prices of derivatives. Moreover, establishing the existence of moments of order higher than one of the process and its approximation is also important for the convergence analysis.
Conditioning on the σ-algebra Define the explosion time of the moment of order ω to be the first time beyond which the moment E [S ω T ] will cease to exist, i.e., If the moment does not explode in finite time, then T * (ω) = ∞.
For convenience, we fix S 0 = 1 and r = 0. Proposition 3.2 derives sharp conditions on the finiteness of moments of the process, whereas Propositions 3.4 and 3.6 to 3.9 give lower bounds on the explosion times of moments of the different discretizations. For illustration, we plot in Figure 1 the explosion time and the corresponding lower bounds with different schemes against the model parameters. Since ∆ = 1 2 ω(ω − 1), (4.5) Propositions 3.2, 3.4 and 3.6 to 3.9 ensure the finiteness of the first moment of the process and its discretizations for all T , i.e., T * (1) = ∞. On the other hand, we infer from Figure  1a that both the explosion time for the exact process and the lower bounds with the explicit Euler discretizations approach infinity as ω approaches one, i.e., that lim ω→1 + T * (ω) = ∞. This ensures the uniform boundedness of moments, for ω sufficiently close to one, of the explicit schemes even for very long maturities, an important ingredient in proving the strong convergence of the approximation process (see Cozma and Reisinger 2015). Note that the green and the yellow curves in Figure 1 overlap when ρ ≥ 0.  The data in Figure 1b suggest that there exists a critical correlation level ρ * such that E [S ω T ] < ∞ for all T , provided ρ ≤ ρ * , and E [S ω T ] = ∞ for some T , provided ρ > ρ * . When k = 0.4, ξ = 0.3 and ω = 2, we find ρ * = −0.04. Moreover, we also infer from the data in Figure 1b that decreasing the correlation has a damping effect on the second moment of the process, and that for strongly negative correlations between the underlying process and the variance, as is usually the case in equity markets (the so-called leverage effect), the lower bounds on the explosion time of the second moment -with all but the reflection schemeare above the typical maturity range of equity derivatives.
The data in Figures 1c and 1d indicate that increasing the speed of mean reversion and decreasing the volatility of volatility have a damping effect on the second moment of the process and its explicit Euler discretizations, which is to be expected considering that larger values of k and smaller values of ξ lead to smaller fluctuations in the variance over time. Next, we assign the following values to the underlying model parameters: k = 0.4, θ = 0.12, ξ = 0.3, v 0 = 0.12 and ρ = 1, and note that the Feller condition is satisfied. Henceforth, we estimate the second moment by a standard Monte Carlo estimator.
From (3.12), the explosion time of the second moment of the process is: T * = 5.77. On the other hand, we infer from the data in Figure 2 that for sufficiently small values of the time step -for instance, when δt = 0.02 -the second moment with the BEM scheme will cease to exist after some time T BEM close to T * . The first sign of moment explosion can be observed when T BEM = 5.98, however this phenomenon is more pronounced when T = 6.14, where the moment jumps to 9.7 × 10 4 .  By close inspection of the data in Figure 3, we infer that for sufficiently small values of the time step δt, the approximation to the second moment with either the implicit or one of the explicit schemes considered in this paper explodes after some critical time close to T * . Combined with the previous observation, this suggests that the explosion times of the second moments of the process and its discretizations become close as we increase the number of time steps. Therefore, we deduce from the data in Figure 1 that the lower bounds for the partial truncation, the full truncation, the absorption and the symmetrized Euler schemes are sharper than the ones for the drift-implicit and the reflection schemes.

Conclusions
In this paper, we have established the uniform exponential integrability of functionals of the Cox-Ingersoll-Ross process and a number of its discretization schemes often encountered in the finance literature. One consequence of this result with obvious practical implications is the stability of moments of numerical approximations for a large class of SDEs arising in finance, which in turn is used to prove strong convergence (Higham et al. 2002;Cozma and Reisinger 2015). An open question is whether we can find sharp conditions on the exponential integrability of Euler approximations for the CIR process.