Abstract
We are interested in the numerical solution of mean-reverting CEV processes that appear in financial mathematics models and are described as nonnegative solutions of certain stochastic differential equations with sublinear diffusion coefficients of the form where . Our goal is to construct explicit numerical schemes that preserve positivity. We prove convergence of the proposed SD scheme with rate depending on the parameter . Furthermore, we verify our findings through numerical experiments and compare with other positivity preserving schemes. Finally, we show how to treat the two-dimensional stochastic volatility model with instantaneous variance process given by the above mean-reverting CEV process.
1. Introduction
Consider the following stochastic models in Itô form: where represents the underlying financially observable variable, is the instantaneous volatility when or the instantaneous variance when , and the Wiener processes , have correlation .
We assume that is a mean-reverting CEV process of the above form, with the coefficients for and , since the process has to be nonnegative. To be more precise the above restriction on implies that is positive; that is, is unattainable, as well as nonexplosive; that is, is unattainable, as can be verified by Feller’s classification of boundaries [17, Proposition 5.22] (see also Appendix A). The steady-state level of is and the rate of mean-reversion is .
System (1) for is the Heston model. When we get the Brennan-Schwartz model [2, Section II] that despite its simple form cannot provide analytical expressions for
Process for , also known as the CIR process [6, (13)], by the initials of the authors that proposed it for the term structure of interest rates, has received a lot of attention and we just mention the latest two contributions to the study of such processes (see [4, 5] and references therein).
Process for has been also considered for the dynamics of the short-term interest rate [3, (1)]. The stationary distribution of the process has also been derived in [7, Proposition 2.2].
We aim for a positivity preserving scheme for the process The scheme that we propose and denote as semidiscrete (SD) preserves the analytical property of staying positive. The idea of the semidiscrete method is that we discretize a part of the original SDE and then apply Itô’s formula (cf. [8] where the method originally appeared and [5, 9, 23]). The explicit Euler scheme fails to preserve positivity, as well as the standard Milstein scheme. We intend to apply the semidiscrete method for the numerical approximation of in model (1) with and compare with other positivity preserving methods such as the balanced implicit method (BIM) (introduced by [11, (3.2)] with the positivity preserving property [25, Section 5] and its stability properties [13]; see also [14] for an extended balanced method with better stability behavior) and the balanced Milstein method (BMM) [25, Theorem 5.9] (we give in Appendix B the form of all the above schemes for the approximation of ). Finally, we approximate the stochastic volatility model (1) with In [15] a thorough treatment can be found, where also another stochastic volatility model is suggested.
Section 2 provides the setting and the main results, Theorems 1 and 2, concerning the -convergence of the proposed semidiscrete method to the true solution of mean-reverting CEV processes of the form of the stochastic volatility in (1). The rate of mean-square convergence in Theorem 1 is logarithmic and in Theorem 2 is polynomial with magnitude The main ingredient of the approach we adopt, inspired by [16], is a change of the initial Brownian motion to another Brownian motion justified by Lévy’s martingale characterization of Brownian motion.
Section 3 is devoted to the logarithmic rate of convergence of the proposed semidiscrete scheme, while Section 4 concerns the proof of the polynomial rate of convergence. In Section 5 we briefly discuss the case where we do not alter the initial Brownian motion This approach produces reduced convergence rate. Finally, Section 6 presents illustrative figures where the behavior of the proposed scheme, regarding the order of convergence, is shown and a comparison with BIM and BMM schemes is given. In Section 7 we treat the full model (1) for a special case. Concluding remarks are in Section 8 and in Appendix B we briefly present numerical schemes for the integration of the variance-volatility process
2. The Setting and the Main Results
We consider the following SDE: where are positive and Then, Feller’s test implies that there is a unique strong solution such that a.s. when a.s. Let where and
Let the partition with and consider the following process: with a.s. or more explicitly for , where represents the level of implicitness and with
Process (6) is well defined when and this is true when and Furthermore, (6) has jumps at nodes Solving for , we end up with the following explicit scheme: with solution in each step given by [1, (4.39), page 123] which has the pleasant feature
Inspired by [16] we remove the term from (6) by considering the process which is a martingale with quadratic variation and thus a standard Brownian motion with respect to its own filtration, justified by Lévy’s theorem [17, Theorem 3.16, page 157]. Therefore, the compact form of (6) becomes for where Consider also the process The process of (2) and the process of (14) have the same distribution. We show in the following that as ; thus the same holds for the unique solution of (2); that is, as To simplify notation we write as We end up with the following explicit scheme: where is as in (8).
Assumption A. Let the parameters be positive such that and consider such that , for Moreover assume a.s. and for some
Theorem 1 (logarithmic rate of convergence). Let Assumption A hold. The semidiscrete scheme (15) converges to the true solution of (2) in the mean-square sense with rate given by where is independent of and given by where
Assumption B. Let Assumption A hold where now and .
Theorem 2 (polynomial rate of convergence). Let Assumption B hold. Then the semidiscrete scheme (15) converges to the true solution of (2) in the mean-square sense with rate given by where and is the constant described in (83) and is an appropriately chosen positive parameter which satisfies (84) and always exists, , and
In the following sections we write for simplicity or for .
3. Logarithmic Rate of Convergence
3.1. Moment Bounds
Lemma 3 (moment bound for SD approximation). It holds that for any , where .
Proof of Lemma 3. We first observe that is bounded in the following way: a.s., where the lower bound comes from the construction of and the upper bound follows from a comparison theorem. We will bound and therefore , since a.s. Set the stopping time , for with the convention Application of Itô’s formula on implies where in the second step we have used the fact that , in the third step the inequality , valid for and with , and in the final step the fact that and Taking expectations in the above inequality and using that is a local martingale vanishing at , we get where we have applied the Gronwall inequality [18, (7)]. We have that Thus taking expectations in the above inequality and using the estimated upper bound for we arrive at and, taking the limit as , we get Let us fix The sequence of stopping times is increasing in and as , and thus the sequence is nondecreasing in and as Application of the monotone convergence theorem implies for any Using again Itô’s formula on , taking the supremum, and then using Doob’s martingale inequality on the diffusion term we bound and thus
Lemma 4 (error bound for SD scheme). Let be an integer such that Then for any , where the positive quantities do not depend on
Proof of Lemma 4. First we take a We get that where we have used the Cauchy-Schwarz inequality. Taking expectations in the above inequality and using Lemma 3 and Doob’s martingale inequality on the diffusion term we conclude where the positive quantity , except on , depends also on the parameters , but not on Now, for , we get where we have used Jensen’s inequality for the concave function Following the same lines, we can show that for any , where the positive quantity , except on , depends also on the parameters , but not on
For the rest of this section we rewrite again the compact form of (12) in the following way: where is given by (3) and the auxiliary process is close to as shown in the next result.
Lemma 5 (moment bounds involving the auxiliary process). For any it holds that and for one has that for any , where the positive quantities do not depend on
Proof of Lemma 5. We have that for any , where we have used (33). Using Lemma 3 we get the left part of (34). Now for and noting that we get the right part of (34), where we have used Lemma 3. The case follows by Jensen’s inequality as in Lemma 4.
Furthermore, for and , we derive that where we have used (30) and in the same manner The case follows by Jensen’s inequality.
3.2. Convergence of the Auxiliary Process to in
We will use the representation (33) and write
Proposition 6. Let Assumption A hold. Then one has for any , where and
Proof of Proposition 6. Let the nonincreasing sequence with and We introduce the following sequence of smooth approximations of (method of Yamada and Watanabe [19]): where the existence of the continuous function with and support in is justified by The following relations hold for with ,We have thatMoreover we find that where we have used properties of Hölder continuous functions and, namely, the fact that is -Hölder continuous for , that is, , and that is -Hölder continuous since Application of Itô’s formula to the sequence implies where in the second step we have used (46) and (47) and the properties of and Taking expectations in the above inequality yields where we have used Lemma 5 in the second step and the Hölder inequality and Lemmas 3 and 4 in the third step and the fact that (the function belongs to the space of real-valued measurable -adapted processes such that ; thus [20, Theorem 1.5.8] implies ). Thus (45) becomes where in the second step we have used the asymptotic relations as for any as for any as ; in the last step we have used the Gronwall inequality and is as defined in Proposition 6 while Taking the supremum over all gives (41).
3.3. Convergence of the Auxiliary Process to in
Proposition 7. Let Assumption A hold. Then one has where is independent of and given by , where
Proof of Proposition 7. We estimate the difference It holds that where in the second step we have used the Cauchy-Schwarz inequality and (46) and Taking the supremum over all and then expectations we have where in the second step we have used Lemma 4 and Doob’s martingale inequality with , since is an -valued martingale that belongs to We find that where we have used (47). Now, Lemmas 3, 4, and 5 imply where we have used the asymptotic relations for all as and the quantity is given by
Relation (56) becomes where we have used Proposition 6 in the second step with the sequence as defined there and Gronwall’s inequality in the last step and the asymptotic relation as , for any , and is independent of and given by
We take , with , to be specified soon and note that as , since as Moreover we have that Now, since there is an small enough such that We take and conclude that as which in turn implies the asymptotic relation as , with the logarithmic rate stated before. In the same way we can show as , by taking We finally arrive at by taking , which implies (53).
3.4. Proof of Theorem 1
In order to finish the proof of Theorem 1 we just use the triangle inequality, Lemma 5, and Proposition 7 to getwhere , given in the statement of Theorem 1.
4. Polynomial Rate of Convergence
We work with the stochastic time change inspired by [21]. We define the process and the stopping time The process is well defined since a.s. and (see Section 2).
The difference is estimated as in Section 3 and we get, as in (56), that where is a stopping time and independent of is as in the proof of Proposition 7. The main difference here will be the estimation of the last term in (66). The approach in Section 3 resulted in the estimation where we used the Yamada-Watanabe approach. Now, we use the Berkaoui approach. Relation (47) becomes where we have used the inequality valid for all , and Consequently, we get the upper bound where we used the Hölder inequality; independent of is as in proof of Proposition 7. Relation (66) becomes where we have used Lemma 5 in the second step. At this point we want to estimate the inverse moments of and to do so we consider the transformation and apply Itô’s formula to get for , where Denote the drift coefficient of the process by and consider the function where Some elementary calculations show that this function attains its minimum at and , thus Consider the process defined through for with Process (74) is a square root diffusion process and when or the process is a CIR process which remains positive if By a comparison theorem [17, Proposition 5.2.18] we obtain that a.s. or a.s. or equivalently a.s. The inverse moment bounds of follow [22, (3.1)]:by choosing big enough and particularly such that (75) holds strictly. Therefore, Relation (77) for implies where in the last step we have used Gronwall’s inequality. Using again relation (77) for and under the change of variables we get where in the last steps we have used (78). We proceed by showing that Markov’s inequality implies for any The following bound holds: and thus where It remains to bound the exponential inverse moments of defined through the stochastic integral equation (2). Exponential inverse moments for the CIR process are known [10, Theorem 3.1] and are given by for where the positive constant is explicitly given in [10, (10)] and depends on the parameters , but is independent of Thus the other condition that we require for parameter is When (84) is satisfied then (75) is satisfied too; thus there is actually no restriction on the coefficient in (83) since we can always choose appropriately a such that (84) holds. Relation (82) becomesWe therefore require that and can always find a , such that the above relation holds by choosing appropriately as discussed before.
Relation (85) becomes and therefore where is chosen such that (86) holds with We conclude by choosing where , is as given in statement of Theorem 2.
5. Alternative Approach with Reduced Rate of Convergence
In this section we briefly discuss the case where instead of (12) we use directly (6). Then, Lemmas 3, 4, and 5 still hold; that is, the moment bounds and error bounds of , as well as the moment bounds involving the auxiliary process , are true. The proof of the convergence results follows the same lines as in Sections 3 and 4. The main difference is in the estimation (47) that now becomes The first term on the right-hand side of the above inequality containing the will contribute in a negative way to the rate of convergence. We do not give all the details but just mention that in order to bound the expectation of that term, which can be done in the following way, we need to estimate the probability of being negative when at the same time , for
Lemma 8. For every it holds that where and and Relation (92) implies that as
Proof of Lemma 8. By the definition (7) of for and for , we have that where The following inclusion relations hold for the event : when and , where We obtain for every standard normal random variable , where in the last step we have used [17, (9.20), page 112] valid for . Using the fact that is a standard normal r.v. and ignoring the exponential term in (96), since its exponent is negative, we get that The following inclusion relations hold for the event : when and Using again (96) we have that Taking probabilities in the inclusion relation (93) and using (97) and (99) we get since as Finally, note that as which justifies the notation (see, e.g., [24]).
6. Numerical Experiments
We discretize the interval with a number of steps in power of The semidiscrete (SD) scheme is given by for , where are the increments of the Brownian motion which are Gaussian random variables with
The ALF (Alfonsi) scheme [4, Section 3] is an implicit scheme which requires solving the nonlinear equation and then computing The estimation of in (102) can be done, for example, with Newton’s method but requires a small enough (in the CIR case, that is, when , (102) simplifies to a solution of a quadratic equation). We also consider a scheme recently proposed in [16] using again the SD method, but in a different way, for Note the similarity in the expressions of (103) and the SD scheme (101) proposed here. This is not strange, because they both rely on the same way of splitting the drift coefficient. In particular, in the explicit HAL scheme, the following process is considered: for with a.s. where now A comparison with (3) and (4) shows that and , for We write (104) again as and the process (106) is well defined when The reader can compare again with (6) for Solving for , we end up with The main result in [16] is when (107) holds, implying a rate of convergence at least which is bigger than the rate of convergence of the SD scheme proposed here which is at least (see Theorem 2).
We also consider two more linear-implicit schemes that were stated in the Introduction and discussed in Appendix B. Namely, we compare with the balanced implicit method (BIM) with appropriate weight functions to guarantee positivity ([25, Theorem 5.9]), which reads and the balanced Milstein method (BMM) with the suggested weight functions [25, Theorem 5.9] that is given by We take the relaxation parameter to be as recommended in [25, (5.10)].
We aim to show experimentally the order of convergence for the above positivity preserving methods for the estimation of the true solution of the CEV model (2), that is, the semidiscrete methods SD (101) and the HAL scheme (103), as well as the implicit ALF scheme (102) and the linear-implicit schemes BIM and BMM. The choice of the parameters is the same as in [15, Figure 6] with In particular
Furthermore, we would also like to reveal the dependence of the order of the semidiscrete methods on ; that is, we want to verify our theoretical results and in particular the order shown in Theorem 2. We take the level of implicitness of the SD method (101) to be ; that is, we consider the fully implicit scheme. We also discuss the fully explicit scheme, that is, when , but also an intermediate scheme , in Section 7.
We want to estimate the endpoint -norm of the difference between the numerical scheme evaluated at step size and the exact solution of (2). For that purpose, we compute batches of simulation paths, where each batch is estimated by and the Monte Carlo estimator of the error is and requires Monte Carlo sample paths. The reference solution is evaluated at step size of the numerical scheme. For the SD case, we have shown in Theorems 1 and 2 that it strongly converges to the exact solution. We simulate paths, where the choice for is as in [28, page 118]. The choice of the number of trajectories is also considered in [26, Section 5] where a fundamental mean-square theorem is proved for SDEs with superlinear growing coefficients satisfying a one-side Lipschitz condition, but unfortunately it is not positivity preserving. Of course, the number of Monte Carlo paths has to be sufficiently large, so as not to significantly hinder the mean-square errors.
We plot in a - scale and error bars represent confidence intervals. The results are shown in Table 1 and Figure 1. Table 1 does not present the computed Monte Carlo errors with confidence, since they were at least times smaller that the mean-square errors.
In Table 2 we present the computational times (we simulate with GHz Intel Pentium, GB of RAM in MATLAB Software. The random number generator is Mersenne Twister. The evaluated times do not include the random number generation time, since all the methods we compare involve the same amount of random numbers) of fully implicit SD, HAL, ALF, BIM, and BMM, for the same problem. Figure 2 shows the relation between the error and computer time consumption. As one can see from Table 2 the CPU times for ALF are at least times bigger than the other schemes; thus we choose in Figure 2 to restrict our attention to the rest of the methods.
We show, in Table 3, the -distance between our proposed method and the other methods for the numerical approximation of (2). We work as before and estimate the distance, between method and , by considering sufficient small , and in particular for
Finally, we examine the behavior of all the methods for a value of the parameter close to The results are shown in Table 4.
The following points of discussion are worth mentioning:(i)The performance of all methods, as shown in Table 1 and Figure 1, implies, in terms of error estimates, that the implicit ALF scheme performs better, for values of discretization steps Actually for these step sizes the ALF method starts to converge, and the same is true for the HAL and BIM methods. All the methods except ALF, that is, the semidiscrete SD and HAL, the BIM, and the BMM have a similar behavior for all values of in the sense of error estimation as Figure 1 shows. The similarity of SD, HAL, BIM, BMM, and especially between SD and BMM is also indicated in Table 3, where we see how close they are with respect to the -norm. Nevertheless, Table 3 also shows that in order to get an accuracy to at least two decimal digits, which in practice may be adequate concerning that we want, for example, to evaluate an option and thus our results are in euros, we are free to use any of the above available methods. We may then choose the fastest one, as will be discussed later on.(ii)The experimental strong order of convergence of implicit SD for problem (2) is (at least as shown theoretically and presented in Table 1). We also see that all methods converge with similar orders and the theoretically rate of the ALF method [4] does not hold for these values of . Thus, again we see that the rate in practical situations does not necessarily matter, if one has to consider very small values of to achieve it. Moreover, we present in Table 5 the performance of the explicit SD method and see that it is very close to the implicit, which is of course natural to happen.(iii)Table 4 concerns the case where the parameter is We do not present the ALF method since it required smaller values of All the methods again behave quite the same, with the BIM performing better with respect to error estimation.(iv)In practice, the computer time consumed to provide a desired level of accuracy is of great importance. In particular, in financial applications, a scheme is considered better when, except of its accuracy, it is implemented faster. As mentioned before, the SD method as well as the HAL method performs well in that aspect, compared to the implicit ALF method, which requires the estimation of a root of a nonlinear equation in each step and is therefore time consuming. This is presented in Table 2 and Figure 2 which illustrates the advantage of the semidiscrete method SD, performing slightly better than HAL and BMM, better than BIM, and of course a lot better compared with ALF (over times quicker to achieve an accuracy of almost two decimal digits). Moreover, the explicit SD performs slightly better in that aspect.(v)A negative step of a numerical method appears when the computer-generated random variable exceeds a certain threshold, which tends to increase as the step size decreases. Thus, the undesirable effect of negative values that are produced by some numerical schemes (such as the explicit Euler (EM) and standard Milstein (M)) tends to disappear, since, after a certain small step size, the threshold exceeds the maximum standard normal random number attainable by the computer system.
7. Approximation of Stochastic Model (1)
So far we have focused on the process , which is one part of system (1). Nevertheless, it can be treated independently, since the only way that it interacts with the process is through the correlation of the Wiener processes. First we apply Itô’s formula on to get
Then, we consider two different schemes for the integration of (113) (the reason for not considering other schemes such as the two-dimensional Milstein is that they generally are time consuming, since they involve additional random number generation for the approximation of double Wiener integrals). The first is the EM scheme which reads which has strong convergence order and is easy to implement. The second scheme, which is based on an interpolation of the drift term and an interpolation of the diffusion term, considering decorrelation of the diffusion term, including a higher order Milstein term [15, Section 4.2], is denoted IJK and is given by [15, (137)]
We therefore consider the EM scheme (114) combined with SD (101) and the IJK scheme (115) combined with SD (101) and compare with the case where the stochastic variance is integrated with BMM scheme (110), for three different correlation parameters, , and with , as in [15, Section 5]. We present in Tables 6, 7, and 8 the errors, in the sense of distance (112), for all the above considered ways of numerical integration of process , for different step sizes, as well as the average computational time (in seconds) consumed for each discretization. We also give an illustrative representation just for one case in Figure 3.
Tables 6, 7, and 8 indicate that in all cases the favorable choice is to integrate using IJK method combined with the SD scheme for in model (1). The IJK-SD approximation of system (1) seems to be the better one, with respect to CPU time, for every correlation coefficient considered.
8. Conclusion
In this paper, we exploit further the semidiscrete method (SD), which originally appeared in [8], to numerically approximate stochastic processes that appear in financial mathematics and are meant to be nonnegative. In [23] we examined the Heston -model, that is, a mean-reverting process with superlinear diffusion, described by a SDE of the form (2) with Now, we deal with SDEs with sublinear diffusion coefficients of the type with These kinds of SDEs, called mean-reverting CEV processes, appear in stochastic models, where they represent the instantaneous volatility-variance of an underlying financially observable variable. We prove theoretically the strong convergence of our proposed SD scheme, revealing the order of convergence. The resulting polynomial rate is shown in Theorem 1. We present a comparative study between various positivity preserving schemes and the SD method seems to be the best with respect to CPU time consumption. The advantage of the SD method here is that although implicit, it has an explicit formula and thus requires fewer arithmetic operations and consequently less computational time. Moreover, our method can cover cases where (2) has time varying coefficients, that is,
We also treat the two-dimensional stochastic volatility model (1). In order to do that, we actually integrate the process which satisfies a SDE of the form (113) and in the end transform back for We only consider two different schemes for the integration of , namely, the Euler Maruyama (EM) scheme, which is easy to implement and the IJK scheme [15, (137)] which is shown to be the most efficient method, robust and simple as EM [15]. We do not apply other two-dimensional schemes, such as, for example, the Milstein scheme, since they are in general time consuming, as they involve approximations of double Wiener integrals which require additional random number generation. We therefore combine the EM scheme with SD ((114) and (101)) and the IJK scheme with SD ((115) and (101)) and compare with the case where the stochastic variance is integrated with BMM scheme (110), for three different correlation parameters, , and with , as in [15, Section 5]. The combination IJK with SD seems to be the most favorable with respect to CPU time, for all the cases.
Appendices
A. Boundary Classification of One-Dimensional Time-Homogeneous SDEs
Let us now recall some results [17, Section 5.5] concerning the boundary behavior of SDEs of the form Let be an interval with and define the exit time from to be Let also the coefficients of (A.1) satisfy the following conditions: Nondegeneracy (ND) Local Integrability (LI)
Then for , we can define the scale function whose behavior at the endpoints of determines the boundary behavior of [17, Proposition 5.22]. In particular, we get for the dynamics of the mean-reverting CEV process of (1) a boundary behavior which is determined by the scale function for any , where Let and take We compute when ; thus by [17, Proposition 5.22c] we have that
B. Some Numerical Schemes for the Integration of the Variance-Volatility Process
We consider a partition of the time interval with and discretization steps for Moreover, we denote by the increments of the Brownian motion. We show in the following subsections some numerical schemes for the approximation of and make some brief comments on them. We also denote
B.1. Standard Euler-Maruyama Scheme
The Euler method, applied to the SDE setting, already appeared in the s through Maruyama [27] and thereafter there has been an extensive study on numerical approximations of solutions of SDEs (we just mention [12] for a recent review on numerical methods for SDEs with applications in finance and references therein).
The explicit Euler-Maruyama (EM) scheme for the process is given by for Clearly ; thus the EM scheme can produce negative values with positive probability, or in the notion of [29] we say that (B.2) has a finite life time.
B.2. Standard Milstein Scheme
The standard one-dimensional Milstein (M) scheme contains some extra terms derived by Itô-Taylor expansion [1, Section 5] and applied to reads for where we have retained terms of order Again (M) scheme has a finite life time.
B.3. Balanced Implicit Method
The balanced implicit method (BIM) [11, (3.2)] was the first attempt to treat the problem of invariance-preserving of specific domains of the underlying process and reads for where and are appropriate weight functions. The choice and preserves positivity [25, Section 5]. Rearranging the above equation, we get the expression
B.4. Balanced Milstein Method
The balanced Milstein method (BMM) was proposed in [25], for an improvement of the BIM in the stability behavior as well as in the rate of convergence. It is given by the following linear-implicit relation: for where and are appropriate weight functions. The choice , where and implies an eternal life time for the scheme [25, Theorem 5.9], in the sense that The step sizes have to be such that The relaxation parameter resembles the implicitness parameter ( in our notation). For there is no restriction in the step size, but it is recommended when possible [25, Remark 5.10] to take Rearranging with the above specifications leads to
Finally, the proposed semidiscrete (SD) scheme reads Increasing the time horizon results in an increase of the percentage of negative paths of EM and M. On the other hand BIM, BMM, and of course SD are not affected by that, since they preserve their positivity on any interval .
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors thank an anonymous referee for helpful comments.