Non-Negativity Preserving Numerical Algorithms for Problems in Mathematical Finance

Yuan Yuan

Hydroforce (Beijing) Technology Co., Ltd., Beijing, China.

DOI: 10.4236/am.2018.93024                

Abstract

We give a study result to analyze a rather different, semi-analytical numerical algorithms based on splitting-step methods with their applications to mathematical finance. As certain subsistent numerical schemes may fail due to producing negative values for financial variables which require non-negativity preserving. These algorithms which we are analyzing preserve not only the non-negativity, but also the character of boundaries (natural, reflecting, absorbing, etc.). The derivatives of the CIR process and the Heston model are being extensively studied. Beyond plain vanilla European options, we creatively apply our splitting-step methods to a path-dependent option valuation. We compare our algorithms to a class of numerical schemes based on Euler discretization which are prevalent currently. The comparisons are given with respect to both accuracy and computational time for the European call option under the CIR model whereas with respect to convergence rate for the path-dependent option under the CIR model and the European call option under the Heston model.

Keywords

Stochastic Differential Equations, Numerical Methods, Splitting-Step Algorithm, Non-Negativity, Convergence, Simulation

1. Introduction

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 [1] [2] [3] [4] ). 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 [5] and the Heston model [6] 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

$\text{d}X\left(t\right)=f\left(X\left(t\right)\right)\text{d}t=\sigma \left(X\left(t\right)\right)\text{d}W\left(t\right),X\left(0\right)=X_0\in D$ (1.1)

are well-defined only with certain boundary conditions. For example, the Cox-Ingersoll-Ross (CIR) [5] model and the Heston model [6] which were proposed to describe the short rate of interest and stock price dynamics respectively are well-defined in the domain $D=\left[0,\infty \right)$ , where $X\left(t\right)=0$ implies an absorbing or reflecting state. Precisely, the CIR process can be expressed in the form

$\text{d}V\left(t\right)=\kappa \left(\theta -V\left(t\right)\right)\text{d}t-\sigma \sqrt{V\left(t\right)}\text{d}{W}^V\left(t\right)$ (1.2)

where $W^V\left(t\right)$ is a Wiener process and $\kappa ,\theta ,\sigma$ are positive constants. The Heston model is a two-factor model of the form:

$\text{d}S\left(t\right)=\mu S\left(t\right)\text{d}t+\sqrt{V\left(t\right)}S\left(t\right)\text{d}{W}^S\left(t\right)$ (1.3)

$\text{d}V\left(t\right)=\kappa \left(\theta -V\left(t\right)\right)\text{d}t+\sigma \sqrt{V\left(t\right)}\text{d}{W}^V\left(t\right)$ (1.4)

where $W^S\left(t\right),W^V\left(t\right)$ are two correlated Wiener processes with $\text{d}{W}^S\left(t\right)\text{d}{W}^V\left(t\right)=\rho \text{d}t,\rho \in \left(-1,1\right)$ , and $\mu ,\kappa ,\theta ,\sigma$ 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 $d:=4\kappa \theta /{\sigma }^2$ , there are significant differences in the boundary behavior of CIR process. According to Feller’s classification [7], if $d\ge 2$ , the boundary $V\left(t\right)=0$ is unattainable and the strong solution is strictly positive; if $d<2$ , 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 [8] ). 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. $d<2$ .

Second, as the square-root term in CIR process avoids the possibility of negative values, Euler-Maruyama [1] and higher order Taylor-type methods [1] [2] cannot be applied in this case as they lead to negative values of $X\left(t\right)$ 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 $\frac{1}{2}-\epsilon$ , Moro and Schurz [9] 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 [10] ) and the Heston model (see Broadie and Kaya [11] [12] ). 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 [13], Kahl and Schurz [14] and Dereich et al., 2012 [15] . There also are some splitting methods that serve for both boundary conditions without any restriction of parameters, such as Alfonsi [16], Moro and Schurz [9] . 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 [17], Bossy and Diop [18], Berkaoui, Bossy and Diop [19], Higham and Mao [20], Lord et al., 2008 [21] . 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 [21] has been shown in practice to be the leading method in this class.

The exact simulation method for CIR model exists, see Glasserman [10], but it requires more computational time than discretization schemes (up to a factor 10). This was analyzed in Alfonsi [13] . For the Heston model, the exact simulation scheme (see Broadie and Kaya [11] [12] ) has the same drawback of highly time-consuming and were analyzed in Broadie and Kaya [11] [12] and Lord et al., 2008 [21] . After comparing to various numerical schemes, Kahl and Jäckel [22] 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 $4\kappa \theta >{\sigma }^2$ : 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 [23] 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 [21] 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 [6] . However, these formulae require the evaluation of logarithms with complex arguments during the involved inverse Fourier integration step. Carr & Madan [24] pioneered the use of the fast Fourier transform (FFT) technique by mapping the Fourier transform directly to option prices via the characteristic function. Lewis [25] 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 [25] was made by considering the resulting integral as a contour integral in the complex plane. Lee [26] generalizes Carr & Madan [24] and unifies it with extensions of some relevant elements and proved an error analysis for these FFT methods. Kahl & Jächel [27] 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 [28] present the optimal contour of the Fourier integral, taking into account numerical issues such as cancellation and explosion of the characteristic function. Fang & Oosterlee [29] 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 [9] 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 [9] has general structure as followed. Suppose that the more general equation to be integrated is of the form

$\text{d}X\left(t\right)=\left[\alpha \left(X\left(t\right),t\right)+\beta \left(X\left(t\right),t\right)\right]\text{d}t+\sigma \left(X\left(t\right),t\right)\text{d}W\left(t\right),$ (2.1)

where $W\left(t\right)$ is a standard Wiener process? Then decompose the above equation into two subsystems

$\text{d}{X}_1\left(t\right)=\beta \left(X_1\left(t\right),t\right)\text{d}t+\sigma \left(X_1\left(t\right),t\right)\text{d}W\left(t\right),$ (2.2)

$\text{d}{X}_2\left(t\right)=\alpha \left(X_2\left(t\right),t\right)\text{d}t$ (2.3)

where it’s required that we know the exact strong solution for $X_1\left(t\right)$ or the conditional probability $P\left[X_1\left(t\right)|X_1\left(0\right)\right]$ . Afterwards, approximate the solution of (2.1) along the time interval $\left[t,\Delta t\right]$ using the following two-step algorithm for each $\Delta t$ .

Ÿ Step 1: Knowing the value $X_t$ , taking it as an initial data of (2.2), i.e. $X_t=X_1\left(t\right)$ , we obtain an intermediate value ${\tilde{X}}_t=X_1\left(t+\Delta t\right)$ through the exact integration of (2.2), or alternatively, through the conditional transition probability $P\left[X_1\left(t+\Delta t\right)|X_1\left(t\right)\right]$ .

Ÿ Step 2: Using X_t obtained in step 1 as the initial condition for (2.3), i.e. ${\tilde{X}}_t=X_2\left(t\right)$ , integrate (2.3) by any converging deterministic numerical algorithm to get $X_2\left(t+\Delta t\right)$ (at least of deterministic order 1). Then $X_{t+\Delta t}=X_2\left(t+\Delta t\right).$

Easily speaking, the procedure is as follow:

We adapt the example in [9] to see how to obtain $X_1\left(t\right)$ through the conditional transition probability $P\left[X_1\left(t\right)|X_1\left(0\right)\right]$ . Suppose that $\beta \left(X_1,t\right)=0$ , $\sigma \left(X_1,t\right)=\sqrt{X_1}$ , then the conditional probability distribution is given by

$\begin{array}{l}P\left[X_1\left(t\right)|X_1\left(0\right)\right]=\frac{2}{t}{\left(\frac{X_1\left(0\right)}{X_1\left(t\right)}\right)}^{\frac{1}{2}}{I}_1\left(\frac{4}{t}\sqrt{X_1\left(t\right)X_1\left(0\right)}\right){\text{e}}^{-\frac{2}{t}\left[X_1\left(t\right)+X_1\left(0\right)\right]}\\ +{\text{e}}^{-\frac{4}{t}{X}_1\left(0\right)}\delta (X_1(t))\end{array}$

where $\delta \left(x\right)$ is the Dirac delta function and I is the modified Bessel function of the first kind with index n which is given by

$I_{\nu }\left(x\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\frac{1}{n!\left(\nu +n+1\right)}{\left(\frac{x}{2}\right)}^{\nu +2n}},x>0$

with gamma function $\Gamma \left(t\right):={\displaystyle {\int }_0^{\infty }{x}^{t-1}{\text{e}}^{-x}\text{d}x},t>0$ specially $\Gamma \left(n\right)=\left(n-1\right)!$ for $n=0,1,2,\cdots$ . 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 $X_1\left(t\right)$ . Noticing that the variable

$Z\left(t\right)=\frac{4}{t}{X}_1\left(t\right)$ follows a non-central c²-distribution, that is

$P\left[Z|Z_0\right]={\displaystyle \underset{j=1}{\overset{\infty }{\sum }}\frac{{\left(\lambda /2\right)}^j{\text{e}}^{-\lambda /2}}{j!}{P}_{{\chi }_{2j}^2}\left(Z\right)}+{\text{e}}^{-\frac{t}{4}\lambda }\delta (Z)$

where $\lambda =\frac{4}{t}{Z}_0$ and $P_{{\chi }_{2j}^2}\left(x\right)$ is c²-pdf with 2j degrees of freedom. Then

$X_1\left(t\right)=\frac{1}{2k}\{\begin{array}{l}0\text{if}K=0,\\ {\displaystyle \underset{i=1}{\overset{2K}{\sum }}{z}_i^2}\text{if}K\ne 0,\end{array}$

where $k=\frac{2}{t}$ and K is chosen from a Poisson distribution with mean $\lambda /2$ and

$z_i$ are independent Gaussian random numbers with zero mean and unit variance.

Generally, a non-central chi-square distribution random variable ${{\chi }^{\prime }}_d^2\left(\lambda \right)$ with d degree of freedom and non-centrality parameter λ whose probability density function is given by

$P\left[{{\chi }^{\prime }}_d^2\left(\lambda \right)=x\right]=p\left(x;d,\lambda \right)=\frac{{\text{e}}^{-\left(\lambda +x\right)/2}}{2}{\left(\frac{x}{\lambda }\right)}^{\left(d-2\right)/4}{I}_{\left(d-2\right)/2}\left(\sqrt{\lambda x}\right)$ (2.4)

see Johnson et al., 1994 [30] . This distribution is properly defined for d positive, and was extended to the case $d=0$ by Siegel [31], afterwards was extended to the case for $d=-2,-4,\cdots$ by Moro and Schurz [9] . They used the fact that the distribution of (0.1) can be expressed as a mixture of central c² variables with Poisson weights, then it can be sample as follow: Choose K from a Poisson distribution with mean = 2, then

${{\chi }^{\prime }}_d^2\left(\lambda \right)={\chi }_{d+2K}^2$

And

${{\chi }^{\prime }}_d^2\left(\lambda \right)=\{\begin{array}{l}0\text{if}d+2K\le 0\\ {{\chi }^{\prime }}_{d+2K}^2\text{if}d+2K>0\end{array}d=0,-2,-4,\cdots$

where ${\chi }_{d+2K}^2$ can be sampled using any standard random number generator of the ${\chi }_d^2$ distribution. This sampling is used specially when $\lambda$ is small. For the case of large $\lambda$ , other approximations (see e.g. Johnson et al., 1994 [30] ) 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 $C^{i,j}\left({ℝ}^d\times \left[0,T\right]\right)$ denote the vector space of continuous functions $f=f\left(x,t\right)$ which are i times continuously differentiable with respect to the space

Coordinate $x_k\in ℝ$ , $\left(k=1,2,\cdots ,d\right)$ and j times continuously differentiable with respect to the time coordinate $t\in \left[0,T\right]$ . Let

$0<t_0<t_1<\cdots <t_n<t_{n+1}<\cdots <t_N=T$

be any random partition of the given time interval $\left[0,T\right]$ with sufficiently small maximum step size $\Delta ={\max }_{i=1,2,\cdots ,N}\left|t_i-t_{i-1}\right|\le 1$ . Then the time discretized approximation ${\tilde{X}}_{\Delta }$ 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 ${\delta }_0>0$ such that the following strong error $\epsilon \left(\Delta \right)$ satisfies

$\epsilon \left(\Delta \right)=\text{E}\left(\left|X\left(T\right)-{\tilde{X}}_{\Delta }\left(T\right)\right|\right)\le C{\Delta }^{\gamma }$

for each $\Delta \in \left(0,{\delta }_0\right)$ .

Along with the strong convergence, the weak convergence can be defined. A discrete-time approximation $Y_{\Delta }$ is said to converge with weak order $\beta >0$ to X at time T as $\Delta \to 0$ if for each smooth function g of polynomial growth there exists a constant $C_g$ , which does not depend on D and ${\Delta }_0\in \left[0,1\right]$ such that the following weak error $ϵ\left(\Delta \right)$ satisfies the estimate

$ϵ\left(\Delta \right)=\left|E_g\left(X\left(T\right)\right)-E_g\left(Y_{\Delta }\left(T\right)\right)\right|\le C_g{\Delta }^{\beta }$

for each $\Delta \in \left(0,{\Delta }_0\right)$ .

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 [9]. We quote the convergence theorem here:

Recall that the original SDE is

$\text{d}X\left(t\right)=\left[\alpha \left(X\left(t\right),t\right)+\beta \left(X\left(t\right),t\right)\right]\text{d}t+\sigma \left(X\left(t\right),t\right)\text{d}W\left(t\right).$ (2.5)

We refer to the splitting

$\text{d}{X}_1\left(t\right)=\beta \left(X\left(t\right),t\right)\text{d}t+\sigma \left(X\left(t\right),t\right)\text{d}W\left(t\right)$ (2.6)

$\text{d}{X}_2\left(t\right)=\alpha \left(X\left(t\right),t\right)\text{d}t$ (2.7)

Theorem Assume that the coefficient functions $\alpha ,\beta \in C^{2,1}\left({ℝ}^d\times \left[0,T\right]\right)$ and $\sigma \in C^{3,2}\left({ℝ}^d\times \left[0,T\right]\right)$ with exclusively uniformly bounded derivatives are such that

$\underset{0\le t\le T}{\sup }E\left[\left|X{\left(t\right)}^2\right|+{\left|\alpha \left(X\left(t\right),t\right)\right|}^2+{\left|\beta \left(X\left(t\right),t\right)\right|}^2+{\left|\sigma \left(X\left(t\right),t\right)\right|}^2\right]<+\infty$

for a fixed finite, nonrandom terminal time $T>0$ . 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 $\left[0,T\right]$ (in the worst case).

The proof and a general theorem on L²-convergence based on variation-of-constants formula (VOP) see Moro and Schurz [9] . 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 [5] 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:

$\text{d}V\left(t\right)=-\kappa \left(V\left(t\right)-\theta \right)\text{d}t+\sigma \sqrt{V\left(t\right)}\text{d}W\left(t\right),\forall t\in {ℝ}^{+}$ (3.1)

with k, θ and σ strictly positive parameters and $W\left(t\right)$ 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, $\kappa \left(\theta -V\left(t\right)\right)$ 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 [32], different k, θ and σ values may produce distinct behavior at the boundary $V\left(t\right)=0$

・ If $\kappa \theta \ge {\sigma }^2/2$ , the upward drift is sufficiently large to make the boundary unattainable, i.e., the solution is always positive $V\left(t\right)>0$ if $V_0>0.$

・ If $\kappa \theta <{\sigma }^2/2$ , there are infinite many values of $t>0$ for which $V\left(t\right)=0$ . 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 $\sigma \sqrt{V\left(t\right)}$ 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 [10], but it’s very time-consuming. A simple Euler discretization of the CIR process may produce negative values of $V\left(t\right)$ . 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 $V\left(t\right)=0$ or setting $V\left(t\right)=-V\left(t\right)$ whenever $V\left(t\right)<0$ (see e.g. [33] ). Other direct approaches base on forced Euler-Maruyama approximations by slight modification of the model are also prevalent, see for example Higham and Mao [20], Deelstra and Delbaen [17], Lord et al., 2008 [21], which the full truncation method in Lord et al., 2008 [21] 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

$0<t_0<t_1<\cdots <t_n<t_{n+1}<\cdots <t_N=T$

be any random partition of the given time interval $\left[0,T\right]$ .

・ Higham and Mao [20] 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.

$V\left(t_{n+1}\right)=V\left(t_n\right)-\kappa \Delta t\left(V\left(t_n\right)-\theta \right)+\sigma \sqrt{\left|V\left(t_n\right)\right|}\Delta W_n,$ (3.2)

$V\left(t_0\right)=V\left(0\right),$ (3.3)

where $\Delta W_n=W\left(t_{n+1}\right)-W\left(t_n\right)$ .

Although it can be applied to the simulation, it leads to negative values of $V\left(t\right)$ in practice as we can see from Figure 1

・ In Lord et al., 2008 [21], 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 $V_0=1,T=1,\kappa =1,\theta =1,\sigma =2$ . The splitting-step algorithm preserves non-negativity for positive initial data while Higham-Mao produces negative values occasionally.

$\tilde{V}\left(t_{n+1}\right)=\tilde{V}\left(t_n\right)-\kappa \Delta t\left(\tilde{V}\left(t_n\right)-\theta \right)+\sigma \sqrt{\left|\tilde{V}\left(t_n\right)\right|}\Delta W_n,$ (3.4)

$V\left(t_{n+1}\right)=\left|\tilde{V}\left(t_{n+1}\right)\right|,$ (3.5)

$V\left(t_0\right)=V\left(0\right).$ (3.6)

・ Deelstra and Delbaen [17] 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 [21] and has the form as follow

$\tilde{V}\left(t_{n+1}\right)=\tilde{V}\left(t_n\right)-\kappa \Delta t\left(\tilde{V}\left(t_n\right)-\theta \right)+\sigma \sqrt{{\left(\tilde{V}\left(t_n\right)\right)}^{+}}\Delta W_n$ (3.7)

$\tilde{V}\left(t_{n+1}\right)={\left(\tilde{V}\left(t_{n+1}\right)\right)}^{+}$ (3.8)

$V\left(t_0\right)=V\left(0\right)$ (3.9)

・ Full truncation was devised by Lord, Koekkoek & Van Dijk in [21] 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

$\tilde{V}\left(t_{n+1}\right)=\tilde{V}\left(t_n\right)-\kappa \Delta t\left(\tilde{V}{\left(t_n\right)}^{+}-\theta \right)+\sigma \sqrt{{\left(\tilde{V}\left(t_n\right)\right)}^{+}}\Delta W_n$ (3.10)

$\tilde{V}\left(t_{n+1}\right)={\left(\tilde{V}\left(t_{n+1}\right)\right)}^{+}$ (3.11)

$V\left(t_0\right)=V\left(0\right)$ (3.12)

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.

$\text{d}{V}_1\left(t\right)=\kappa \theta \text{d}t+\sigma \sqrt{V_1\left(t\right)}\text{d}W\left(t\right)$ (3.13)

$\text{d}{V}_2\left(t\right)=-\kappa V_2\left(t\right)\text{d}t$ (3.14)

We note that the process defined by (3.13) is a kθ-dimensional squared Bessel

process (BESQ) with index of the process = $\nu =\frac{2\kappa \theta }{{\sigma }^2}-1$ (see Makarov and

Glew [34] ). The boundary $V_1\left(t\right)=0$ is entrance if $\nu \ge 0$ , regular if $-1<\nu <0$ , or exit if $\nu \le -1$ . For the regular diffusion on $\left(0,\infty \right)$ the transition probability density function (PDF) is given by

$P\left(V_1\left(t+\Delta t\right)|V_1\left(t\right)\right)={\left(\frac{V_1\left(t+\Delta t\right)}{V_1\left(t\right)}\right)}^{\frac{\nu }{2}}\frac{{\text{e}}^{-2\left(V_1\left(t+\Delta t\right)+V_1\left(t\right)\right)/{\sigma }^2\Delta t}}{{\sigma }^2\Delta t/2}{I}_{\nu }\left(\frac{4\sqrt{V_1\left(t+\Delta t\right)V_1\left(t\right)}}{{\sigma }^{2}t}\right)$ (3.15)

in the case of a regular boundary, where $I_{\nu }$ is the modified Bessel function of the first kind with index n. Then we notice that the transition density $P\left[V_1\left(t+\Delta t\right)|V_1\left(t\right)\right]$ can be written in terms of a non-central ${\chi }^2$ 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 λ:

$P\left[{{\chi }^{\prime }}_d^2\left(\lambda \right)=x\right]=p\left(x;d,\lambda \right)=\frac{{\text{e}}^{-\left(\lambda +x\right)/2}}{2}{\left(\frac{x}{\lambda }\right)}^{\left(d-2\right)/4}{I}_{\left(d-2\right)/2}\left(\sqrt{\lambda x}\right),x>0$

Comparing the equation above and (3.15), we have

$P\left(V_1\left(t+\Delta t\right)|V_1\left(t\right)\right)=\frac{4}{{\sigma }^2\Delta t}p\left(\frac{4V_1\left(t+\Delta t\right)}{{\sigma }^2\Delta t};\frac{4\kappa \theta }{{\sigma }^2},\frac{4V_1\left(t\right)}{{\sigma }^2\Delta t}\right),$ (3.16)

as $\nu =\frac{d}{2}-1$ . The Equation (3.16) shows that the random process $\frac{4V_1\left(t+\Delta t\right)}{{\sigma }^2\Delta t}$ can be represented by a non-central chi-square distribution with $d=\frac{4\kappa \theta }{{\sigma }^2}$ degree of freedom and non-centrality parameter $\lambda =\frac{4V_1\left(t\right)}{{\sigma }^2\Delta t}$ . Hence we have

$V_1\left(t+\Delta t\right)=\frac{{\sigma }^2\Delta t}{4}{{\chi }^{\prime }}_d^2\left(\lambda \right)$ (3.17)

where

$\lambda =\frac{4V_1\left(t\right)}{{\sigma }^2\Delta t},d=\frac{4\kappa \theta }{{\sigma }^2}$ (3.18)

Then along the time interval $\left[t,t+\Delta t\right]$ , given $V\left(t\right)$ , first we take it as an initial data of (3.13) and integrate the SDE though the transition conditional probability (3.15) to obtain $V_1\left(t+\Delta t\right)$ which can be done by sampling a ${{\chi }^{\prime }}_d^2\left(\lambda \right)$ random number in (3.17). Second we regard $V_1\left(t+\Delta t\right)$ 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 ${{\chi }^{\prime }}_d^2\left(\lambda \right)$ 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 $V\left(t\right)$ in the squart-root term. Here we choose $\kappa =1,\theta =1,\sigma =2$ , 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 $T>0$ .

We take g the identity function and use the model

$\text{d}V\left(t\right)=\left(a+bV\left(t\right)\right)\text{d}t+\sigma \sqrt{V\left(t\right)}\text{d}W\left(t\right)$ (3.19)

Then split it in this way:

$\text{d}{V}_1\left(t\right)=a\text{d}t+\sigma \sqrt{V_1\left(t\right)}\text{d}W\left(t\right)$ (3.20)

$\text{d}{V}_2\left(t\right)=b{V}_2\left(t\right)\text{d}t$ (3.21)

Similarly to (3.17) (3.18), we have

$V_1\left(t+\Delta t\right)=\frac{{\sigma }^2\Delta t}{4}{{\chi }^{\prime }}_d^2\left(\lambda \right)$ (3.22)

where

$\lambda =\frac{4V_1\left(t\right)}{{\sigma }^2\Delta t},d=\frac{4a}{{\sigma }^2}$ (3.23)

We use (3.22) (3.23) to obtain $V_1\left(t+\Delta t\right)$ 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 $V\left(0\right)=V_0\ge 0$ .

Taking the expected value of both sides of (3.19) reads:

$\begin{array}{l}\text{d}\mathrm{E}\left[V\left(t\right)|V\left(0\right)=V_0\right]\\ =\left(a+b\mathrm{E}\left[V\left(t\right)|V\left(0\right)=V_0\right]\right)\text{d}t+\mathrm{E}\left[\sigma \sqrt{V\left(t\right)}\text{d}W\left(t\right)|V\left(0\right)=V_0\right]\end{array}$ (3.24)

Notice that the second term on the right side $\mathrm{E}\left[\sigma \sqrt{V\left(t\right)}\text{d}W\left(t\right)|V\left(0\right)=V_0\right]=0$ as $W\left(t\right)$ is a standard Wiener process. We define

$y\left(t\right)=\mathrm{E}\left[V\left(t\right)|V\left(0\right)=V_0\right]$

Then $y\left(t\right)$ is a deterministic function of t. Equation (3.24) becomes

$\text{d}y\left(t\right)=\left(a+by\left(t\right)\right)\text{d}t$

with $y\left(0\right)=V_0$ . Solving the above equation we have

$y\left(t\right)=y\left(0\right){\text{e}}^{bt}+\frac{a}{b}\left({\text{e}}^{bt}-1\right)$

Therefore,

$\mathrm{E}\left[V\left(t\right)|V\left(0\right)=V_0\right]=y\left(t\right)=V_0{\text{e}}^{bt}+\frac{a}{b}\left({\text{e}}^{bt}-1\right)$ (3.25)

For $a=1,b=1$ , we have

$\mathrm{E}\left[V\left(t\right)|V\left(0\right)=V_0\right]={\text{e}}^t\left(1+V_0\right)-1$ (3.26)

Then we apply the splitting-step method to obtain the numerical approximation $\tilde{V}\left(t\right)$ . Figure 2 depicts the error

${ϵ}_{\Delta t}=\left|\mathrm{E}\left[V\left(t\right)\right]-\mathrm{E}\left[\tilde{V}\left(t\right)\right]\right|$ (3.27)

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 $T=1.0$ and $S\left(0\right)=1.0$ , $a+b=1$ , $\sigma =2$ using the splitting-step, Higham-Mao complemented, partial truncation and full truncation algorithms in log-log scale for 10⁷ 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 $K>0$ is called the exercise (or strike) price and $T>0$ the maturity (or expiration) time. We pay K to buy the asset at time T if $V\left(T\right)\ge K$ and sell it immediately at $V\left(T\right)$ making a profit of $V\left(T\right)-K$ . The option is worthless if $V\left(T\right)<K$ as we can buy it cheaper than K.

Thus at time T, the call option has the value

$C_T=\mathrm{E}\left[{\left(V\left(T\right)-K\right)}^{+}\right]=\mathrm{E}\left[\max \left(V\left(T\right)-K,0\right)\right]$ (3.28)

Suppose that C_T is the numerical approximation of C_T, then the error

${ϵ}_{call}=\left|{\tilde{C}}_T-C_T\right|$ (3.29)

reveals the precision of an algorithm. Since the exact solution C_T no explicitly known, we compute a numerical reference solution with a very small t and regard it as C_T. Here we use the numerical approximation of C_T with $t={10}^4$ 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 [35] [36] ). 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 $\max \left(V\left(T\right)-K,0\right)$ 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 $V_0$ , the price of U&O call option at expiry time T is given by

$C^{U\&O}\left(T\right)=\mathrm{E}[{\left(V\left(T\right)-K\right)}^{+}{\text{1}}_{\left\{0<V\left(t\right)<L,0\le t\le T\right\}}]$ (3.30)

where $1_{\left\{x\right\}}$ is one if x is true and zero otherwise.

The exact option price $C^{U\&O}\left(T\right)$ 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.

(a) (b)

Figure 3. (a) The call option valuation error (3.29) as a function of Dt for $V\left(0\right)=1.0$ , $\kappa =\theta =1$ , $\sigma =2$ , $K=2$ , $T=10$ using the splitting-step, Higham-Mao complemented, partial truncation and full truncation algorithms in log-log scale for 10⁷ 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 ${\tilde{u}}_{\Delta t}$ 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

$\left|{\tilde{u}}_{\Delta t}-u\right|\le C\Delta t^p$

At least for sufficiently small Dt. It’s also said that the convergence rate of the method is $\Delta t^p$ . Normally the error ${\tilde{u}}_{\Delta t}-u$ depends smoothly on Dt. Then

${\tilde{u}}_{\Delta t}-u=C\Delta t^p+O\left(\Delta t^{p+1}\right)$

i.e.,

${\tilde{u}}_{\Delta t}=u+C\Delta t^p+O\left(\Delta t^{p+1}\right)$

An approach of p is to check the relative differences between ${\tilde{u}}_{\Delta t}$ computed for different Dt. In most cases we compare solutions where Dt is halved successively. Then we get

$\begin{array}{c}{\tilde{u}}_{\Delta t}-{\tilde{u}}_{\Delta t/2}=u+C\Delta t^p-u-C{\left(\Delta t/2\right)}^p+O\left(\Delta t^{p+1}\right)\\ =C\Delta t^p\left(1-\frac{1}{2^p}\right)+O\left(\Delta t^{p+1}\right)\\ =C^{\prime }\Delta t^p+O\left(\Delta t^{p+1}\right)\end{array}$

with $C^{\prime }=C\left(1-\frac{1}{2^p}\right)$ Hence, we can get an estimate of the order of accuracy p after computing ${\tilde{u}}_{\Delta t}$ and ${\tilde{u}}_{\Delta t/2}$ for different Dt.

The relative differences ${\tilde{u}}_{\Delta t}-{\tilde{u}}_{\Delta t/2}$ 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:

${ϵ}_{\Delta t}^{U\&O}=\left|{\tilde{C}}_{\Delta t}^{U\&O}\left(T\right)-{\tilde{C}}_{\Delta t/2}^{U\&O}\left(T\right)\right|$ (3.31)

where ${\tilde{C}}_{\Delta t}^{U\&O}\left(T\right)$ is the numerical approach of $C_{\Delta t}^{U\&O}\left(T\right)$ using time step Dt. Thus ${ϵ}_{\Delta t}^{U\&O}=C\Delta t^p+O\left(\Delta t^{p+1}\right)$ , with C a constant independent of Dt. We show the relation between ${ϵ}_{\Delta t}^{U\&O}$ 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 ${ϵ}_{\Delta t}^{U\&O}$ . For sufficient small Dt, the relative differences ${ϵ}_{\Delta t}^{U\&O}$ 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 ${ϵ}_{\Delta t}^{U\&O}$ . But roughly $\mathrm{SEM}\sim \frac{1}{\sqrt{N}}$ 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 $N={10}^7$ sample paths, it’s impossible in practice to increase it to at least $N={10}^{11}$ to make ${ϵ}_{\Delta t}^{U\&O}$ larger than SEM. Then for sufficient small Dt, we estimate ${ϵ}_{\Delta t}^{U\&O}$ 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 $V\left(0\right)=1.0$ , $\kappa =\theta =1$ , $\sigma =2$ , $T=10$ , $K=2$ , $L=10$ and 10⁷ sample paths are generated. Figure 5(b) portrays the computational time in seconds with respect to the relative difference ${ϵ}_{\Delta t}^{U\&O}$ . At the beginning, splitting performs the worst as it’s the most time-consuming one, see for example when ${ϵ}_{\Delta t}^{U\&O}=\text{0}\text{.02}$ . 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 ${ϵ}_{\Delta t}^{U\&O}<\text{0}\text{.001}$ from the tendency of these four algorithms.

4. Option Pricing with Heston’s Stochastic Volatility Model

4.1. Heston Model Basics

The Heston model [6] is defined by the coupled two-dimensional stochastic differential equation (SDE):

(a) (b)

Figure 5. (a) The relative difference (3.31) of values of up-and-out option at time T as a function of Dt for $T=10$ , $T=10$ , $T=10$ , $\sigma =2$ , $K=2$ , $L=10$ , using the splitting-step, Higham-Mao complemented, partial truncation and full truncation algorithms in log-log scale for 10⁷ realizations. Dash lines are p power laws of Dt and $\sqrt{\Delta t}$ . (b) Computational time in seconds as a function of the relative difference (3.31) for the same parameter setting.

$\text{d}S\left(t\right)=\mu S\left(t\right)\text{d}t+\sqrt{V\left(t\right)}S\left(t\right)\text{d}{W}^S\left(t\right)$ (4.1)

$\text{d}V\left(t\right)=\kappa \left(\theta -V\left(t\right)\right)\text{d}t+\sigma \sqrt{V\left(t\right)}\text{d}{W}^V\left(t\right)$ (4.2)

where $\kappa ,\theta ,\sigma$ are positive constants, $W^S\left(t\right)$ and $W^V\left(t\right)$ are Wiener processes with correlation, i.e., $\text{d}{W}^S\left(t\right)\text{d}{W}^V\left(t\right)=\text{d}t,\rho \in \left(-1,1\right)$ . 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 $V\left(t\right)$ tends to θ. k is the rate at which $V\left(t\right)$ reverts to θ σ is the so-called vol of vol, which determines the variance of $V\left(t\right)$ . We note that the variance (4.2) follows a CIR process.

For $V\left(t\right)$ it has the same boundary behavior as we mentioned in §3.1.

1) 0 is an attainable boundary when ${\sigma }^2>2\kappa \theta$ . The boundary is strongly reflecting;

2) ∞ is an unattainable boundary.

There is an additional condition for $S\left(t\right)$ which is iii) $S\left(t\right)$ 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 [11] [12] derived an exact simulation method without bias, but highly time-consuming. This is analyzed in Broadie and Kaya [21], and Lord et al., 2008 [21] . 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 [22] 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 $4\kappa \theta >{\sigma }^2$ : 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 [23] were attested to outperform the full truncation method proposed by Lord et al., 2008 [21] 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 $S\left(t\right)$ , as in Lord et al., 2008 [21], i.e.

$\ln S\left(t+\Delta t\right)=\ln S\left(t\right)+\left(\mu -\frac{1}{2}V\left(t\right)\right)\Delta t+\sqrt{V\left(t\right)}\left(\rho \Delta W_V\left(t\right)+\sqrt{1-{\rho }^2}\Delta W_S\left(t\right)\right)$ (4.3)

with $W_S\left(t\right),W_V\left(t\right)$ two independent Wiener processes. For Higham-Mao, as it may produce negative values of $V\left(t\right)$ with substantial probability which we have seen in §3.2.1, we follow their spirit to take the absolute value of $V\left(t\right)$ in (4.3), i.e.,

$\ln S\left(t+\Delta t\right)=\ln S\left(t\right)+\left(\mu -\frac{1}{2}\left|V\left(t\right)\right|\right)\Delta t+\sqrt{\left|V\left(t\right)\right|}\left(\rho \Delta W_{V\left(t\right)}+\sqrt{1-{\rho }^2}\Delta W_S\left(t\right)\right)$ (4.4)

4.3. A European Call Option with Heston Model

Closed form semi-analytical formulae for plain vanilla option prices have been derived in [6] . However, these formulae require the evaluation of logarithms with complex arguments during the involved inverse Fourier integration step. Carr & Madan [24], Lewis [25] considered the use of the fast Fourier transform (FFT) technique by mapping the Fourier transform directly to option prices via the characteristic function. Lee [26] generalizes Carr & Madan [24] and unifies it with extensions of some relevant elements and proved an error analysis for these FFT methods. Kahl & Jächel [27] 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 [28] present the optimal contour of the Fourier integral, Fang & Oosterlee [29] 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 $C^{Hes}\left(T\right)$ to be the exact European call option price at maturity time T based on Heston model, i.e.

$C^{Hes}\left(T\right)=\mathrm{E}\left[{\left(S\left(T\right)-K\right)}^{+}\right]=\mathrm{E}\left[\max \left(S\left(T\right)-K\right),0\right]$

and ${\tilde{C}}_{\Delta t}^{Hes}$ is an numerical approximation of $C^{Hes}\left(T\right)$ with time step Dt.

Then the relative difference is defined as

${ϵ}_{\Delta t}^{Hes}=\left|{\tilde{C}}_{\Delta t}^{Hes}-{\tilde{C}}_{\Delta t/2}^{Hes}\right|$ (4.5)

The same technique with §3.3.2 gives that ${ϵ}_{\Delta t}^{Hes}=C\Delta t^p+O\left(\Delta t^{p+1}\right)$ , with C a constant independent of Dt, where p is the order of accuracy defined in §3.3.2. We estimate p by computing ${ϵ}_{\Delta t}^{Hes}$ for different Dt. The slope of the plots of ${ϵ}_{\Delta t}^{Hes}$ with respect to Dt gives the approximation of p.

Figure 6 are drawn with parameters $S\left(0\right)=100$ , $V\left(0\right)=0.04$ , $\mu =0.02$ , $\kappa =1.5$ , $\theta =0.06$ , $\sigma =0.7$ , $\rho =0$ , $K=120$ , $T=10$ and 10⁷ 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 ${ϵ}_{\Delta t}^{Hes}>2$ 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 ${ϵ}_{\Delta t}^{Hes}=0.065$ 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 ${ϵ}_{\Delta t}^{Hes}=0.065$ . From Figure 7, one can find that for reducing ${ϵ}_{\Delta t}^{Hes}$ to 0.065, splitting needs around 10⁴ seconds » 2.8 hours while Higham-Mao needs much more than 10⁶ seconds » 277.8 hours which is more than 100 times than splitting.

5. Conclusions

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 [1] 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 ${ϵ}_{\Delta t}^{Hes}=0.065$ .

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.

(a) (b)

Figure 6. (a) Relative differences of call option (4.5) as a function of Dt in log-log scale with $S\left(0\right)=100$ , $V\left(0\right)=0.04$ , $\mu =0.02$ , $\kappa =1.5$ , $\theta =0.06$ , $\sigma =0.7$ , $\rho =0$ , $K=120$ , $T=10$ , and 10⁷ sample paths. Dash lines are power laws of Dt and $\sqrt{\Delta t}$ . (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 ${ϵ}_{\Delta t}^{Hes}$ .

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 [1] .

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 [37] and mixed stock & bond portfolio problem like in [38] are under consideration.

As Moro and Schurz [9] 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.

Acknowledgements

I would take the opportunity to thank my tutor Esteban Moro Egido who gives the support during our project.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Kloeden, P.E. and Platen, E. (1998) Numerical Solution of Stochastic Differential Equations. Springer, New York.
[2]	Kloeden, P.E., Platen, E. and Schurz, H. (1994) Numerical Solution of SDEs through Computer Experiments. Springer, New York. https://doi.org/10.1007/978-3-642-57913-4
[3]	Schurz, H. (1997) Stability, Stationarity, and Boundedness of Some Implicit Numerical Methods for Stochastic Differential Equations and Applications. Logos-Verlag, Berlin.
[4]	Milstein, G.N. (1995) Numerical Integration of Stochastic Differential Equations. Kluwer, Dordrecht. https://doi.org/10.1007/978-94-015-8455-5
[5]	Cox, J.C., Ingersoll, J.E. and Ross, S.A. (1985) A Theory of Term Structure of Interest Rates. Econometrica, 53, 385-407.
[6]	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. https://doi.org/10.1093/rfs/6.2.327
[7]	Feller, W. (1951) Two Singular Diffusion Problems. Annals of Mathematics, 54, 173-182. https://doi.org/10.2307/1969318
[8]	Karlin, S. and Taylor, H.M. (1981) A Second Course in Stochastic Processes. Academic Press, New York.
[9]	Moro, E. and Schurz, H. (2007) Boundary Preserving Semi-Analytic Numerical Algorithms for Stochastic Differential Equations. SIAM Journal on Scientific Computing, 29, 1525-1549. https://doi.org/10.1137/05063725X
[10]	Glasserman, P. (2003) Monte Carlo Methods in Financial Engineering. Springer Verlag, New York. https://doi.org/10.1007/978-0-387-21617-1
[11]	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. https://doi.org/10.1109/WSC.2004.1371506
[12]	Broadie, M. and Kaya, O. (2006) Exact Simulation of Stochastic Volatility and Other Affine Jump Diffusion Processes. Operations Research, 54, 217-231. https://doi.org/10.1287/opre.1050.0247
[13]	Alfonsi, A. (2005) On the Discretization Schemes for the CIR (and Bessel Squared) Processes. Monte Carlo Methods and Applications, 11, 355-384. https://doi.org/10.1515/156939605777438569
[14]	Kahl, C. and Schurz, H. (2006) Balanced Milstein Methods for Ordinary SDEs. Monte Carlo Methods and Applications, 12, 143-170. https://doi.org/10.1515/156939606777488842
[15]	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. https://doi.org/10.1098/rspa.2011.0505
[16]	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. https://doi.org/10.1090/S0025-5718-09-02252-2
[17]	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. https://doi.org/10.1002/(SICI)1099-0747(199803)14:1<77::AID-ASM338>3.0.CO;2-2
[18]	Bossy, M. and Diop, A. (2004) An Efficient Discretization Scheme for One Dimensional SDEs with a Diffusion Coefficient Function of the Form. INRIA Working Paper No. 5396.
[19]	Berkaoui, A., Bossy, M. and Diop, A. (2008) Euler Scheme for SDEs with Non-Lipchitz Diffusion Coefficient: Strong Convergence. ESAIM Probability and Statistics, 12, 1-11.
[20]	Higham, D.J. and Mao, X. (2005) Convergence of Monte Carlo Simulations Involving the Mean-Reverting Square Root Process. Journal of Computational Finance, 8, 35-62.
[21]	Lord, R., Koekkoek, R. and van Dijk, D. (2008) A Comparison of Biased Simulation Schemes for Stochastic Volatility Models. Journal of Quantitative Finance, 10, 177-194.
[22]	Kahl, C. and Jackel, P. (2006) Fast Strong Approximation Monte-Carlo Schemes for Stochastic Volatility Models. Quantitative Finance, 6, 513-536. https://doi.org/10.1080/14697680600841108
[23]	Andersen, L. (2008) Simple and Efficient Simulation of the Heston Stochastic Volatility Model. Journal of Computational Finance, 11, 1-42. https://doi.org/10.21314/JCF.2008.189
[24]	Carr, P. and Madan, D. (1999) Option Valuation using the Fast Fourier Transform. Journal of Computational Finance, 2, 61-73. https://doi.org/10.21314/JCF.1999.043
[25]	Lewis, A.L. (2001) A Simple Option Formula for General Jump-Diffusion and Other Exponential Levy Processes. http://www.optioncity.net/pubs/ExpLevy.pdf
[26]	Lee, R.W. (2004) Option Pricing by Transform Methods: Extensions, Unification and Error Control. Journal of Computational Finance, 7, 51-86. https://doi.org/10.21314/JCF.2004.121
[27]	Kahl, C. and Jackel, P. (2005) Not-So-Complex Logarithms in the Heston Model. Wilmott Magazine, 19, 94-103.
[28]	Lord, R. and Kahl, C. (2007) Optimal Fourier Inversion in Semi-Analytical Option Pricing. http://papers.tinbergen.nl/06066.pdf
[29]	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. https://doi.org/10.1137/080718061
[30]	Johnson, N.L., Kotz, S. and Balakrishnan, N. (1994) Continuous Univariate Distributions. Wiley, New York, 2.
[31]	Siegel, A.F. (1979) The Noncentral Chi-Squared Distribution with Zero Degrees of Freedom and Testing for Uniformity. Biometrika, 66, 381-386. https://doi.org/10.1093/biomet/66.2.381
[32]	Feller, W. (1971) Probability Theory and Its Applications. 2nd Edition, Wiley, New York.
[33]	Gatheral, J. (2006) The Volatility Surface: A Practitioners Guide. John Wiley and Sons, New York.
[34]	Makarov, R.N. and Glew, D. (2010) Exact Simulation of Bessel Discussions. Monte Carlo Methods and Applications, 16, 283-306. https://doi.org/10.1515/mcma.2010.010
[35]	Shreve, S.E. (1996) Stochastic Calculus and Finance. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.137.6951&rep=rep1&type=pdf
[36]	Wilmott, P., Howison, S. and Dewynne, J. (1995) The Mathematics of Financial Derivatives. Cambridge University Press, New York. https://doi.org/10.1017/CBO9780511812545
[37]	Merton, R.C. (1969) Lifetime Portfolio Selection under Uncertainty: The Continuous-Time Case. Review of Economics and Statistics, 51, 247-257. https://doi.org/10.2307/1926560
[38]	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. https://doi.org/10.1137/S0363012900377791