Stochastic differential equations (SDEs) are fundamental in mathematical finance. Particularly, they serve as models for describing the evolution of certain financial variables, such as the stock price, interest rates or volatility of an asset. Although there are extensive studies on SDEs, explicit solutions of SDEs are rarely known or very difficult to obtain, so numerical approximations are relied on. Actually a wealth of numeric methods have been proposed and tested (see for example     ). In our work, we mainly dedicate to applying rather different numerical algorithms to evaluate financial derivatives with comparison to methods based on Euler discretization. We concentrate on the Cox-Ingersoll-Ross (CIR) model  and the Heston model  as they are two of the fundamental models in mathematical finance.
1.1. Problem Statement and Motivation
In some cases, the Itô-type SDEs of the form
are well-defined only with certain boundary conditions. For example, the Cox-Ingersoll-Ross (CIR)  model and the Heston model  which were proposed to describe the short rate of interest and stock price dynamics respectively are well-defined in the domain , where implies an absorbing or reflecting state. Precisely, the CIR process can be expressed in the form
where is a Wiener process and are positive constants. The Heston model is a two-factor model of the form:
where are two correlated Wiener processes with , and are positive parameters. The component S describes the evolution of a financial variable such as stock index or exchange rate, and V describes the stochastic variance of its returns. One can notice that the component V in the Heston model evolves according to the CIR process (1.2).
First let’s mention that the SDE (1.2) for CIR process is not explicitly solvable but its transition probability is known; it can be represented by a non-central chi-square density. Depending on the number of degrees of freedom , there are significant differences in the boundary behavior of CIR process. According to Feller’s classification , if , the boundary is unattainable and the strong solution is strictly positive; if , the boundary is attracting and attainable but strongly reflecting, i.e., each sample path will reflect instantaneously into the positive domain once it reaches 0 (see e.g. Karlin & Taylor  ). The reflecting behavior in the latter case is particularly difficult to capture numerically. During our work, we only focus on the case that the boundary is attainable, i.e. .
Second, as the square-root term in CIR process avoids the possibility of negative values, Euler-Maruyama  and higher order Taylor-type methods   cannot be applied in this case as they lead to negative values of in practice. The numerical methods which preserve non-negativity are preferred. Actually various numerical schemes have been devised to meet this requirement, among which the schemes based on Euler discretization normally gives high efficiency, but their strong convergence rate and discretization errors are difficult to establish. In contrast to the Euler-type methods giving strong convergence of
order , Moro and Schurz  provide the splitting-step methods of strong
convergence at least of order 1. Although they are more costly with respect to computational time, they give higher accuracy and convergence rate. With this motivation, we would like to evaluate some important financial instruments using splitting-step methods with comparison to Euler-type numerical schemes. We apply them to option pricing including a path-dependent option valuation comparing both accuracy and computational cost.
Third, exact simulation methods exist for both the CIR process (see Glasserman  ) and the Heston model (see Broadie and Kaya   ). The drawback of these exact simulation methods is that they are very time-consuming. They are competitive when one just simulates the process at one time (or few times), for example to price a European option with a Monte-Carlo algorithm. On the contrary, they are too slow when one has to simulate the process along a time-grid, for example to compute a path-dependent European option. Thus, the splitting-step methods may have advantage on pricing path-dependent options as they are more efficient.
1.2. Review of Research Status
For the CIR process, a number of numerical schemes based on implicit time-stepping integrators have been devised for the case of unattainable boundary condition, see for example Alfonsi , Kahl and Schurz  and Dereich et al., 2012  . There also are some splitting methods that serve for both boundary conditions without any restriction of parameters, such as Alfonsi , Moro and Schurz  . The latter paper gives a splitting-step method which highly relied on the knowledge of explicit solution of SDEs. This method is what we mainly use throughout our project.
Other direct approaches that can be applied to both attainable and unattainable zero boundary cases are based on modification of Euler-Maruyama approximations. See for example Deelstra and Delbaen , Bossy and Diop , Berkaoui, Bossy and Diop , Higham and Mao , Lord et al., 2008  . These methods all are proved to converge to the exact solution as the time step tends to 0, but the strong convergence rate and discretization errors are difficult to establish. The full truncation method in Lord et al., 2008  has been shown in practice to be the leading method in this class.
The exact simulation method for CIR model exists, see Glasserman , but it requires more computational time than discretization schemes (up to a factor 10). This was analyzed in Alfonsi  . For the Heston model, the exact simulation scheme (see Broadie and Kaya   ) has the same drawback of highly time-consuming and were analyzed in Broadie and Kaya   and Lord et al., 2008  . After comparing to various numerical schemes, Kahl and Jäckel  conclude that combining the balanced Milstein method (BMM) for the variance process with their bespoke IJK method for the logarithm of the stock price gives the best results with respect to strong convergence measure. The BMM method actually preserves positivity for the variance process if : But this restriction is rarely satisfied in practice, and one typically finds that the sampling scheme for V produces negative values with substantial probability.
Anderson  considered two new discretization schemes: the truncated Gaussian scheme (TG) and the quadratic-exponential scheme (QE) with additional martingale-corrected versions labeled by “TG-M” and “QE-M” respectively. Both TG and QE schemes were attested to outperform the full truncation method proposed by Lord et al., 2008  with respect to biases and only cost marginally more computational time than the Euler scheme. The TG scheme is more natural to use than the QE scheme in multi-asset applications that involve several correlated variance processes, but generally performs somewhat worse than QE at practically relevant time steps. Thus they concluded that the QE scheme should be the default choice of the schemes considered due to its simplicity and strong performance; martingale correction (the QE-M scheme) is optional.
Further, for pricing financial derivatives under Heston model, closed form semi-analytical formulae for plain vanilla option prices have been derived in  . However, these formulae require the evaluation of logarithms with complex arguments during the involved inverse Fourier integration step. Carr & Madan  pioneered the use of the fast Fourier transform (FFT) technique by mapping the Fourier transform directly to option prices via the characteristic function. Lewis  considered the same type of approach, except that there the transform was taken with respect to the log-spot price. An important step in Lewis  was made by considering the resulting integral as a contour integral in the complex plane. Lee  generalizes Carr & Madan  and unifies it with extensions of some relevant elements and proved an error analysis for these FFT methods. Kahl & Jächel  propose a new approach which enables the use of Heston’s analytics for practically all levels of parameters and even maturities of many decades since most implementations of Heston’s formulae are not robust for moderate to long dated maturities or strong mean reversion. Lord & Kahl  present the optimal contour of the Fourier integral, taking into account numerical issues such as cancellation and explosion of the characteristic function. Fang & Oosterlee  proposed the COS method which is based on the Fourier and Fourier-cosine expansion.
1.3. Summary of the Paper
In the remaining paper, we introduce the general structure and properties of splitting-step algorithm proposed by Moro and Shurtz in  in Section 2, including an efficient sample manner for generating a non-central chi-square distribution random number.
In Section 3, we present numerical results of applications to CIR model, evaluating a European plain vanilla call option and a path-dependent option with comparisons to the methods based on Euler discretization. The comparisons are with respect to both accuracy and computational time. All of our numerical results are generated on a MacBook Pro with an Intel Core i7 2.3 GHz processor, 16 GB 1600 MHz DDR3 memory, using Xcode 6.4 and R 3.1.0 in a Mac OS X Yosemite environment.
In Section 4 we study the option valuation under the Heston model. The results show that the splitting-step method gives the best convergence rate for pricing a European plain vanilla call option among the five algorithms utilized in our project.
Finally we conclude and point out the future work in Section 5.
2. General Structure of Splitting-Step Algorithms
2.1. Construction of the Algorithms
The splitting-step algorithm  has general structure as followed. Suppose that the more general equation to be integrated is of the form
where is a standard Wiener process? Then decompose the above equation into two subsystems
where it’s required that we know the exact strong solution for or the conditional probability . Afterwards, approximate the solution of (2.1) along the time interval using the following two-step algorithm for each .
Step 1: Knowing the value , taking it as an initial data of (2.2), i.e. , we obtain an intermediate value through the exact integration of (2.2), or alternatively, through the conditional transition probability .
Step 2: Using Xt obtained in step 1 as the initial condition for (2.3), i.e. , integrate (2.3) by any converging deterministic numerical algorithm to get (at least of deterministic order 1). Then
Easily speaking, the procedure is as follow:
We adapt the example in  to see how to obtain through the conditional transition probability . Suppose that , , then the conditional probability distribution is given by
where is the Dirac delta function and I is the modified Bessel function of the first kind with index n which is given by
with gamma function specially for . Then we can sample the conditional probability distribution function using the rejection or inverse methods but it’s computational expensive. There is a more simple way to obtain . Noticing that the variable
follows a non-central c2-distribution, that is
where and is c2-pdf with 2j degrees of freedom. Then
where and K is chosen from a Poisson distribution with mean and
are independent Gaussian random numbers with zero mean and unit variance.
Generally, a non-central chi-square distribution random variable with d degree of freedom and non-centrality parameter λ whose probability density function is given by
see Johnson et al., 1994  . This distribution is properly defined for d positive, and was extended to the case by Siegel , afterwards was extended to the case for by Moro and Schurz  . They used the fact that the distribution of (0.1) can be expressed as a mixture of central c2 variables with Poisson weights, then it can be sample as follow: Choose K from a Poisson distribution with mean = 2, then
where can be sampled using any standard random number generator of the distribution. This sampling is used specially when is small. For the case of large , other approximations (see e.g. Johnson et al., 1994  ) might be taken into account for improving the speed of the splitting-step algorithms.
2.2. Strong and Weak Convergence of Splitting-Step Methods
Let denote the vector space of continuous functions which are i times continuously differentiable with respect to the space
Coordinate , and j times continuously differentiable with respect to the time coordinate . Let
be any random partition of the given time interval with sufficiently small maximum step size . Then the time discretized approximation of a continuous-time process X, is said to be of general strong order of convergence g to X at time T if there exists a positive constant C, which does not depend on D, and a such that the following strong error satisfies
for each .
Along with the strong convergence, the weak convergence can be defined. A discrete-time approximation is said to converge with weak order to X at time T as if for each smooth function g of polynomial growth there exists a constant , which does not depend on D and such that the following weak error satisfies the estimate
for each .
The splitting-step algorithm has been shown that has strong and weak order 1.0 of convergence under certain assumption of the coefficient functions in . We quote the convergence theorem here:
Recall that the original SDE is
We refer to the splitting
Theorem Assume that the coefficient functions and with exclusively uniformly bounded derivatives are such that
for a fixed finite, nonrandom terminal time . Then the splitting-step algorithm with step 1 and 2 (see section §2.1) has (global) strong and weak order 1.0 of convergence on the interval (in the worst case).
The proof and a general theorem on L2-convergence based on variation-of-constants formula (VOP) see Moro and Schurz  . We check the weak convergence with numerical experiment by taking g the identity function in section 3, where we present numerical studies of various strategies for integrating SDEs.
3. Numerical Studies with the CIR Model
3.1. CIR Model and Its Properties
The famous Cox-Ingersoll-Ross (CIR) model was introduced in 1985 by John C. Cox, Jonathan E. Ingersoll and Stephen A. Ross in  as an extension of the Vasicek model. It serves to describe the evolution of interest rates. This model specifies that the interest rate dynamics follows the stochastic differential equation, which is also called CIR Process:
with k, θ and σ strictly positive parameters and a Wiener process. The parameter k determines the speed of adjustment, θ is the long-run mean and σ is the so-called volatility to volatility. The drift factor, is exactly the same as in the Vasicek model. It ensures mean reversion of the interest rate towards the long-term value θ, with speed of adjustment governed by the strictly positive parameter k. According to Feller’s boundary criteria, see Feller , different k, θ and σ values may produce distinct behavior at the boundary
・ If , the upward drift is sufficiently large to make the boundary unattainable, i.e., the solution is always positive if
・ If , there are infinite many values of for which . The boundary becomes attainable, but it is strongly reflecting. That is, when a sample path reaches 0, then it returns immediately to the positive domain in a reflecting manner.
Since for all positive values of k and θ, the standard deviation factor rejects the possibility of negative interest rates, integration strategies must preserve non-negativity. Otherwise, it not only lacks any possible interpretation in the context of finance but also could induce severe errors in option valuation.
3.2. Simulation Schemes for the CIR model
We now turn to the simulation of CIR model (12). The exact simulation method for CIR model exists, see Glasserman , but it’s very time-consuming. A simple Euler discretization of the CIR process may produce negative values of . Practitioners often choose either a absorption or reflection fix whenever the process attains a negative value. That is to say, they deal with it by either setting or setting whenever (see e.g.  ). Other direct approaches base on forced Euler-Maruyama approximations by slight modification of the model are also prevalent, see for example Higham and Mao , Deelstra and Delbaen , Lord et al., 2008 , which the full truncation method in Lord et al., 2008  in practice has been proved to be the leading method of this class. Apart from Euler-type schemes, the splitting-step methods we introduced in the foregoing section can be applied to CIR model without any modification of the SDE. Moreover, the non-negativity is preserved during all the time interval.
3.2.1. Numerical Schemes Based on Euler Discretization
Here we outline the numerical schemes based on Euler discretization which we mentioned above. Let
be any random partition of the given time interval .
・ Higham and Mao  considered an Euler discretization of the CIR model with a novel fix, for which they prove strong convergence. Their scheme deals with the squart-root term in (12) by taking absolute value of the inside component, i.e.
Although it can be applied to the simulation, it leads to negative values of in practice as we can see from Figure 1
・ In Lord et al., 2008 , when they referred the method Higham and Mao, they had one more step which was taking the absolute value of the right side of (3.2), we call this method Higham and Mao complemented, which is as follow
Figure 1. Trajectories of CIR model with . The splitting-step algorithm preserves non-negativity for positive initial data while Higham-Mao produces negative values occasionally.
・ Deelstra and Delbaen  proposed an Euler-type scheme by taking positive part of the component inside the squared-root, this scheme is called partial truncation in Lord et al., 2008  and has the form as follow
・ Full truncation was devised by Lord, Koekkoek & Van Dijk in  which is the leading method of these classes. The difference between full truncation and partial truncation is the treatment of the drift term, where the full truncation method has one more x and can be expressed as follow
We refer the above four numerical schemes as Higham and Mao, Higham and Mao complemented, partial truncation and full truncation respectively.
3.2.2. Splitting-Step Algorithms Applied to CIR Process
According to the general structure of splitting-step algorithm in Section 2, we split the SDE (3.1) into two equations as follows.
We note that the process defined by (3.13) is a kθ-dimensional squared Bessel
process (BESQ) with index of the process = (see Makarov and
Glew  ). The boundary is entrance if , regular if , or exit if . For the regular diffusion on the transition probability density function (PDF) is given by
in the case of a regular boundary, where is the modified Bessel function of the first kind with index n. Then we notice that the transition density can be written in terms of a non-central distribution. Recalling the probability density function (PDF) of a non-central chi-square distribution random variable with d degree of freedom and non-centrality parameter λ:
Comparing the equation above and (3.15), we have
as . The Equation (3.16) shows that the random process can be represented by a non-central chi-square distribution with degree of freedom and non-centrality parameter . Hence we have
Then along the time interval , given , first we take it as an initial data of (3.13) and integrate the SDE though the transition conditional probability (3.15) to obtain which can be done by sampling a random number in (3.17). Second we regard has the initial data of (3.14) and integrate it by any deterministic numerical method of at least order 1. In our project, we practically sample the random number by the aid of the function rnchisq() inside the Rmath library from R and integrate (3.14) with deterministic Euler method.
In Figure 1, three trajectories for the CIR process are shown with a comparison to Higham-Mao which is simply taking absolute value of in the squart-root term. Here we choose , so that the boundary is attainable.
The figure depicts that the splitting-step algorithm preserves non-negativity for positive initial data while Higham-Mao produces negative values occasionally. Thus splitting method is preferable as the CIR process is not defined for negative values.
3.2.3. Weak Convergence
In this section we verify that the splitting-step algorithm has weak convergence of order 1.0 on the interval [0,T] for a fixed finite, nonrandom terminal time .
We take g the identity function and use the model
Then split it in this way:
Similarly to (3.17) (3.18), we have
We use (3.22) (3.23) to obtain which is corresponding to (3.20) and deterministic Euler method to integrate (3.21).
Now, we calculate the expected value of exact solution of (3.19) with the initial value .
Taking the expected value of both sides of (3.19) reads:
Notice that the second term on the right side as is a standard Wiener process. We define
Then is a deterministic function of t. Equation (3.24) becomes
with . Solving the above equation we have
For , we have
Then we apply the splitting-step method to obtain the numerical approximation . Figure 2 depicts the error
Versus decreasing uniform step size Dt, which shows that the method has weak order 1.0 while using constant step sizes Dt. Moreover, we compare to Higham-Mao complemented, partial truncation and full truncation for the same equation. Figure 2 shows that splitting methods has the highest precision among these four algorithms.
3.3. Option Valuations
Among the variety of financial derivatives, the option is one of the most important financial instruments. In current financial markets, there are mainly four kinds of options: American option, European option, Asian option, and Barrier option. In our work, we only focus on pricing European options, which can only be exercised on the maturity date whereas an American option can be exercised at any time before expiration.
Figure 2. Value of weak convergence error (3.27) as a function of Dt for and , , using the splitting-step, Higham-Mao complemented, partial truncation and full truncation algorithms in log-log scale for 107 realizations. The dashed line is proportional to Dt.
3.4. A European Call Option
A European call option on an asset is a contract that allows the buyer to buy (but not obligate) this asset at the price K at time T. The number is called the exercise (or strike) price and the maturity (or expiration) time. We pay K to buy the asset at time T if and sell it immediately at making a profit of . The option is worthless if as we can buy it cheaper than K.
Thus at time T, the call option has the value
Suppose that CT is the numerical approximation of CT, then the error
reveals the precision of an algorithm. Since the exact solution CT no explicitly known, we compute a numerical reference solution with a very small t and regard it as CT. Here we use the numerical approximation of CT with computed by splitting method as reference. Figure 3(a) depicts the comparison of the four algorithms with respect to the precision of calculating call option for each step size t and Figure 3(b) shows the computational cost versus precision among them. We observe that Higham-Mao complemented, partial truncation and full truncation work better on this call option pricing as they have more precision and higher speed.
3.4.1. A Path-Dependent Option
Apart from the typical European call option in the previous section, more exotic options exist in the market (see, for example   ). Here we consider one of path-dependent options: up-and-out call option (hereafter referred as U&O option). An up-and-out call option pays off the usual at expiry unless at any time during the life of the option the underlying asset reaches certain level (from below, obviously) then it is said to knock out and become worthless. The Figure 4 explains the payoff of each underlying asset. We assume that a underlying will be knocked-out once it reaches L from below at any time prior to expiry time T. Then the underlying asset corresponds to the blue one has payoff 0 as it’s knocked-out.
As we say that an U&O option expires worthless if the asset price touches some barrier L from below, say, at any time prior to expiry, where L is larger than the present asset value , the price of U&O call option at expiry time T is given by
where is one if x is true and zero otherwise.
The exact option price is not available since the exact solution of U&O call option under CIR model is not explicitly known. Then instead of evaluating biases we work on the convergence properties of the algorithms.
Figure 3. (a) The call option valuation error (3.29) as a function of Dt for , , , , using the splitting-step, Higham-Mao complemented, partial truncation and full truncation algorithms in log-log scale for 107 realizations. The dashed line is proportional to Dt. (b) Computational time in seconds as a function of the call option valuation error for the same parameter setting.
Figure 4. Payoff of each sample path. The underlying asset corresponds to the blue one has payoff 0 as it’s knocked-out.
Let which depends on the time step Dt be a numerical approximation of an exact value u. The numerical method said to be of order p means that there exists a number C independent of Dt such that
At least for sufficiently small Dt. It’s also said that the convergence rate of the method is . Normally the error depends smoothly on Dt. Then
An approach of p is to check the relative differences between computed for different Dt. In most cases we compare solutions where Dt is halved successively. Then we get
with Hence, we can get an estimate of the order of accuracy p after computing and for different Dt.
The relative differences indicate the convergence rate of an algorithm. The faster converges a scheme (has higher p value), the faster relative differences reduce to 0 as Dt tends to 0. Now we define the relative difference of U&O call option value as:
where is the numerical approach of using time step Dt. Thus , with C a constant independent of Dt. We show the relation between and the time step Dt to get an estimate of p for the four algorithms. We also analyze the computational cost with respect to the relative difference . For sufficient small Dt, the relative differences computed by splitting method are much less than standard error of the mean (SEM) which suggests to increase substantially the number of sample paths for
the sake of validity of these . But roughly where N is size of
the sample, which means that the SEM will reduce to 10% if we increase N to 100N. As we have already used sample paths, it’s impossible in practice to increase it to at least to make larger than SEM. Then for sufficient small Dt, we estimate computed by splitting by the fit of existing data. This has implied that splitting convergent the fastest as it’s the first one to reach the SEM level. Figure 5(a) depicts that the splitting-step method has the highest order of accuracy p since p corresponds to the slope of these lines. The parameter what we used are , , , , , and 107 sample paths are generated. Figure 5(b) portrays the computational time in seconds with respect to the relative difference . At the beginning, splitting performs the worst as it’s the most time-consuming one, see for example when . But along the trend of smaller relative differences, splitting performs better and better and costs the least since certain value of relative difference. For example, with the same parameter setting, we can deduce that splitting will cost the least for the cases of from the tendency of these four algorithms.
4. Option Pricing with Heston’s Stochastic Volatility Model
4.1. Heston Model Basics
The Heston model  is defined by the coupled two-dimensional stochastic differential equation (SDE):
Figure 5. (a) The relative difference (3.31) of values of up-and-out option at time T as a function of Dt for , , , , , , using the splitting-step, Higham-Mao complemented, partial truncation and full truncation algorithms in log-log scale for 107 realizations. Dash lines are p power laws of Dt and . (b) Computational time in seconds as a function of the relative difference (3.31) for the same parameter setting.
where are positive constants, and are Wiener processes with correlation, i.e., . The parameters μ is the rate of return of the asset. θ is the long variance, which means as t tends to infinity, the expected value of tends to θ. k is the rate at which reverts to θ σ is the so-called vol of vol, which determines the variance of . We note that the variance (4.2) follows a CIR process.
For it has the same boundary behavior as we mentioned in §3.1.
1) 0 is an attainable boundary when . The boundary is strongly reflecting;
2) ∞ is an unattainable boundary.
There is an additional condition for which is iii) has an absorbing barrier at 0.
4.2. Path Simulation
There are plenty of methods can be used to simulate Heston model. Broadie and Kaya   derived an exact simulation method without bias, but highly time-consuming. This is analyzed in Broadie and Kaya , and Lord et al., 2008  . The exact scheme is competitive when one simulates the process just at one time (or few times), for example to price European options with a Monte-Carlo algorithm. On the contrary, they are drastically too slow for simulating the process along a time-grid, which occurs when computing path-dependent options prices. Kahl and Jäckel  state that combining the balanced Milstein method (BMM) for the variance process with their bespoke IJK method for the logarithm of the stock price gives the best results with respect to strong convergence measure. The BMM method actually preserves positivity for the variance process if : But this restriction is rarely satisfied in practice, and one typically finds that the sampling scheme for V produces negative values with substantial probability. The TG and QE schemes proposed by Anderson  were attested to outperform the full truncation method proposed by Lord et al., 2008  with respect to biases and only cost marginally more computational time than the Euler scheme.
Comparing to the methods above, Euler discretization and splitting-step are very simple to implement. For simulating the variance process (4.2) we use Splitting-step, Higham-Mao, partial truncation and full truncation methods introduced in §3. Then we switch to logarithms for the asset price , as in Lord et al., 2008 , i.e.
with two independent Wiener processes. For Higham-Mao, as it may produce negative values of with substantial probability which we have seen in §3.2.1, we follow their spirit to take the absolute value of in (4.3), i.e.,
4.3. A European Call Option with Heston Model
Closed form semi-analytical formulae for plain vanilla option prices have been derived in  . However, these formulae require the evaluation of logarithms with complex arguments during the involved inverse Fourier integration step. Carr & Madan , Lewis  considered the use of the fast Fourier transform (FFT) technique by mapping the Fourier transform directly to option prices via the characteristic function. Lee  generalizes Carr & Madan  and unifies it with extensions of some relevant elements and proved an error analysis for these FFT methods. Kahl & Jächel  propose a new approach which enables the use of Hestons analytics for practically all levels of parameters and even maturities of many decades since most implementations of Hestons formulae are not robust for moderate to long dated maturities or strong mean reversion. Lord & Kahl  present the optimal contour of the Fourier integral, Fang & Oosterlee  proposed the COS method which is based on the Fourier and Fourier-cosine expansion.
For simplicity and efficiency, in this section we mainly work on the convergence performance of the algorithms we mentioned in §4.2 and demonstrate that besides the well performance in convergence rates, splitting has superiority in computational efficiency.
We define to be the exact European call option price at maturity time T based on Heston model, i.e.
and is an numerical approximation of with time step Dt.
Then the relative difference is defined as
The same technique with §3.3.2 gives that , with C a constant independent of Dt, where p is the order of accuracy defined in §3.3.2. We estimate p by computing for different Dt. The slope of the plots of with respect to Dt gives the approximation of p.
Figure 6 are drawn with parameters , , , , , , , , and 107 sample paths. Obviously the splitting-step method has the highest order of accuracy p. Although the splitting method is more time-consuming than the others for each Dt, and with respect to relative differences, it costs the most at the beginning as for the case illustrated in Figure 6(b). But it costs the least computational time if one restricts that the relative differences should be less than a certain value. For example, if one requires that the relative difference should be less than 0.10 with the same parameter setting, splitting practically works much faster than the others. The Figure 6(b) gives the evidence.
More precisely, we analyze the computational cost by taking the case as an example. As partial truncation and full truncation give fluctuant points (see Figure 6(b)), we choose Higham-Mao to compare with. The fitted line of Higham-Mao is extended to intersect with the vertical line . From Figure 7, one can find that for reducing to 0.065, splitting needs around 104 seconds » 2.8 hours while Higham-Mao needs much more than 106 seconds » 277.8 hours which is more than 100 times than splitting.
A number of existing numerical algorithms for integrating SDEs cannot be used to simulate certain financial models which require non-negativity, such as the CIR model, the Heston model that we have seen. A slight modification of Euler-Maruyama  method conducts to various Euler-type numerical schemes. They are simple to implement and efficient enough, but have lower accuracy compared to exact simulation schemes. Sometimes they also give (4.5). The vertical dash line corresponds to .
Less precision is compared to semi-analytical numerical schemes such as the splitting-step methods which we are analyzing in our project. Since the exact simulation schemes are thoroughly too slow for simulating a process along a time-grid, splitting-step methods preponderate them with respect to efficiency.
Figure 6. (a) Relative differences of call option (4.5) as a function of Dt in log-log scale with , , , , , , , , , and 107 sample paths. Dash lines are power laws of Dt and . (b) Computational time in seconds as a function of the call option valuation relative differences for the same parameter setting.
Figure 7. Computaional cost in seconds as a function of relative difference .
The splitting algorithms heavily rely on the exploitation of the specific structure of original system (2.1). One can decompose the original system into two subsystems appropriately for which either one knows the explicit solution or the conditional transition probability of one of the subsystems and integrate the remaining one by deterministic numerical methods of at least order 1. In this way, it preserves the non-negativity and a maximum of convergence order 1.0 both in strong and weak sense. For the analytical solution of SDEs, an extensive list of known one can be found in textbooks such as Kloeden and Platen  .
Among the numerical schemes in §3, splitting-step method gives the best convergence rate and normally highest accuracy. We observe that it doesn’t perform better than other four schemes on call option valuation with CIR model. Despite it costs more computational time to generate sample paths, it may cost even less for certain × errors. For example, if one requires that the relative difference of A European call option under the Heston model should be less than 0.10 with the parameter setting in ×4.3, the splitting method costs the least computational time than others.
6. Future Work
Since derivatives are one of the three main categories of financial instruments where the other two being stocks (i.e., equities or shares) and debt (i.e., bonds and mortgages), and as we said that the splitting methods are much more efficient but has less accuracy than exact simulation methods which are not suitable for pathwise simulations, then one direction of our interest is applying splitting-step method to other path-dependent options, probably with comparison to other efficient algorithms like what we have done in §3.3.2.
Another direction is the application of splitting-step models to portfolio problems which are aiming to find the optimal investment strategy of an investor. Following the optimal portfolio strategy leads to the maximum expected utility of the terminal wealth. For example the classical Merton’s problem  and mixed stock & bond portfolio problem like in  are under consideration.
As Moro and Schurz  have mentioned that the splitting-step scheme is applicable to multi-dimensional problems, we intend to encounter some suitable splitting manners to deal with them which in particular can be SDEs with linear, decoupled, commutative or diagonal noise terms. Furthermore, non-linear partial differential equations, for instance, Super-Brownian motion interest us to verify its properties by implement splitting-step methods.
I would take the opportunity to thank my tutor Esteban Moro Egido who gives the support during our project.
 Heston, S. (1993) A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Reviews of Financial Studies, 6, 327-343.
 Moro, E. and Schurz, H. (2007) Boundary Preserving Semi-Analytic Numerical Algorithms for Stochastic Differential Equations. SIAM Journal on Scientific Computing, 29, 1525-1549.
 Broadie, M. and Kaya, O. (2004) Exact Simulation of Option Greeks under Stochastic Volatility and Jump Diffusion Models. In: Ingalls, R.G., Rossetti, M.D., Smith, J.S. and Peters, B.A., Eds., Proceedings of the 2004 Winter Simulation Conference, INFORMS, Washington, 1607-1615.
 Dereich, S., Neuenkirch, A. and Szpruch, L. (2012) An Euler-Type Method for the Strong Approximation of the Cox-Ingersoll-Ross Process. Proceedings of the Royal Society of London A, 468, 1105-1115.
 Alfonsi, A. (2010) High Order Discretization Schemes for the CIR Process: Application to Affine Term Structure and Heston Models. Mathematics of Computation, American Mathematical Society, 79, 209-237.
 Deelstra, G. and Delbaen, F. (1998) Convergence of Discretized Stochastic (Interest Rate) Processes with Stochastic Drift Term. Applied Stochastic Models and Data Analysis, 14, 77-84.
 Lewis, A.L. (2001) A Simple Option Formula for General Jump-Diffusion and Other Exponential Levy Processes.
 Lord, R. and Kahl, C. (2007) Optimal Fourier Inversion in Semi-Analytical Option Pricing.
 Fang, F. and Oosterlee, C.W. (2008) A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions. SIAM Journal on Scientific Computing, 31, 826-848.
 Shreve, S.E. (1996) Stochastic Calculus and Finance.
 Korn, R. and Kraft, H. (2002) A Stochastic Control Approach to Portfolio Problems with Stochastic Interest Rates. SIAM Journal on Control and Optimization, 40, 1250-1269.