An iterative splitting method for pricing European options under the Heston model

An option is a financial instrument that provides the investor with the opportunity to buy or sell an asset at a specified price. Options make it possible for investors to benefit from the return of the underlying asset at a low cost and play an important role in hedging risk as well.

Research on option pricing can be traced back to Bachelier’s work in 1900 [1], since then an increasing number of mathematical tools have been introduced into the financial field [29]. In 1973, the seminal work of Black and Scholes opened a new era of the option pricing theory [2]. After that, the option market expanded rapidly and option pricing attracted considerable attention of researchers.

Due to the volatility smile and skew reflected by the market data, the assumption of constant volatility in the Black-Scholes model need relaxing and different models have been proposed to better reflect the market features based on the Black-Scholes model. Assume that the volatility of the underlying is a function of the price and time, then we get to the well-accepted local volatility model [13] [11]. By calibrating it to the market prices of European options, we can obtain the implied volatility surface which can be used to price other options. But it always leads to the behavior of the volatility smile opposite to the market as the price changes[19]. In jump diffusion models, Poisson processes are added into the dynamic of the underlying to capture the phenomena of the volatility [28, 12, 25]. However, they are not applicable to riskless hedging. Lévy processes are continuous-time stochastic processes with stationary independent increments and are used to model the return of the underlying [9, 18, 7]. The Black-Scholes model and jump diffusion models can be seen as special cases of models based on Lévy processes. Because of the inherent property of Lévy processes, i.e., independent increments, they are not able to capture the dependence structure. Stochastic volatility models are constructed with an additional stochastic process to describe the dynamic of the volatility [24, 31, 21, 30, 20, 17]. Among them, Heston’s stochastic volatility model is the most popular and widely used in both academia and industry. The volatility in the Heston model follows the Cox-Ingersoll-Ross (CIR) process, which ensures that it cannot be negative. The dependence structure can be incorporated and it has a closed-form solution when pricing European options.

As we have shown above, different continuous-time stochastic processes are introduced to model the price or other state variables such as the volatility of the underlying and always appear in literature as stochastic differential equations (SDEs). Under the risk neutral probability measure, the price of the option can be represented as the conditional expectation of the discounted payoff determined by the contract. By the Feynman-Kac theorem, it can be considered as a solution of a partial differential equation (PDE) corresponding to the dynamics or SDEs we assume. Closed-form solutions are difficult to derived or calculated for some of the models and options so that a variety of numerical methods have been developed. By lattice methods, the continuous distribution is transformed into a discrete distribution, which is much easier to understand and implement [10, 5]. However, lattice methods are not suitable for complicated models. To deal with the continuous-time stochastic processes or SDEs directly, we have the well-known Monte Carlo simulation method [4, 24, 3]. It is flexible and applies to high-dimensional problems, but it gives solutions with probable errors and requires a large number of simulation paths to guarantee the accuracy. Researchers have tried to use numerical integration methods to compute the price of the option based on the Fourier Transform with the characteristic functions of the distributions of the underlying [8, 14]. These methods are efficient and can be useful for calibration, while they always need careful derivation of the formulae for different models and options. Another class of methods is to price options by solving the PDEs and a wide range of numerical techniques are then available including finite difference methods [6, 22, 33, 23, 16], finite element methods [15, 32] and so on. PDE methods especially finite difference methods are intuitive and interpretable, but memory space and computational time they cost may increase rapidly as the dimension of the problem grows.

In this paper, we take the Heston stochastic volatility model as an example and try to give a new numerical method for the option pricing problem, which can maintain the advantages of finite difference methods and reduce the influence of the dimension of the PDE. Compared with the Black-Scholes model, the Heston model involved the CIR process to describe the volatility of the underlying, which is a mean-reverting process and can be seen as adding fluctuations to a given mean volatility to some extent. Accordingly, the PDE derived from the Heston model can be seen as adding terms related to the volatility to the PDE corresponding to the Black-Scholes model, and the dimension is then increased from one to two. On this basis, we propose an iterative splitting method to split the operators and solve the quasi one-dimensional PDEs iteratively. We make use of the analytical solution of the Black-Scholes model, lower the computational cost and maintain the simplicity, flexibility and interpretability of the PDE methods. A mixed method is also used to guarantee the convergence of the iterative process. We conduct numerical experiments to price the European options. Boundary conditions from [20] and artificial boundary conditions obtained with the methods in [27] are used for a comparison. Numerical results show that the iterative splitting method together with an artificial boundary condition based on the method from [27] gives the most accurate solutions and Greeks. By contrast with a classic two-dimensional finite difference scheme, its solutions have smaller errors as well. Besides, the iterative splitting method can be easily extended to other models and options.

The paper is organized as follows. Section 2 presents the PDE derived from the Heston stochastic volatility model. The iterative splitting method is proposed in Section 3 and the corresponding algorithms are also given. In Section 4, we derive and show different artificial boundary conditions with the methods in [27]. Section 5 gives the results of numerical experiments with different parameters, boundary conditions and numerical methods. We summarize in Section 6.

Denote the asset spot price and its variance by $S(t)$ and $v(t)$. According to [20], assume that $S(t)$ and $v(t)$ satisfy the following two stochastic processes:

The parameter $\mu$ is the expected return of the underlying, $\kappa$ is the speed at which the variance $v(t)$ reverts to its long run mean $\theta$ and $\sigma$ is the so-called ”volatility of volatility”. In addition, $W^{(1)}(t)$ and $W^{(2)}(t)$ are standard Wiener processes with a constant correlation coefficient $\rho$.

The price of a financial derivative $C(S,v,t)$ can then be considered as the solution of a Garman’s type PDE:

where $r$ is the risk-free interest rate and the price of volatility risk $\lambda(S,v,t)$ can be set to be 0.

We take the European call option as an example and denote the maturity and the strike by $T$ and $K$. Let $\tau=T-t$ represent the time to maturity and

then the PDE (3) is equivalent to

Besides, to give the price of the derivative, an initial value and proper boundary conditions are also needed. Referring to [20], the conditions can be rewritten as follows with the transforms mentioned above:

First, we split the operators on the left-hand side of the PDE (4) into two parts. Let

then the PDE (4) can be rewritten as

If $\mathscr{L}_{2}U$ is known and can be denoted by $Q(\tilde{S},v,\tau)$, the two-dimensional problem (4)-(9) with a mixed spatial-derivative term can be transformed into two quasi one-dimensional problems with the same operator. Let

we get to the following two problems:

and

We can easily give the solution to the problem (13) based on the Black-Scholes formula:

Here, $N(x)$ represents the cumulative distribution function of the standard normal distribution.

The problem (14) can be solved only if $Q(\tilde{S},v,\tau)=\mathscr{L}_{2}U$ is known, so we try to use an iterative process to find the solution. Since the operators in the original problem have been split as above, we call it an iterative splitting method.

Choose a time step size $\Delta\tau=T/N$, and let $\tau_{n}=n\cdot\Delta\tau,\ n=0,...,N$. Denote $U^{n}(\tilde{S},v,\tau)=U(\tilde{S},v,\tau_{n})$ and define $U_{1}^{n}$, $U_{2}^{n}$and $Q^{n}$ similarly. For the $n$-th time step, we first calculate $U_{1}^{n}$ by (15). With the known $U^{n-1}$, we compute $Q^{n-1}=\mathscr{L}_{2}U^{n-1}$ and use it to approximate $Q^{n}$. Then we can solve the problem (14) numerically to get $U_{2}^{n}$. Since $U^{n}=U^{n}_{1}+U^{n}_{2}$, $U^{n}$ can be updated, with which we obtain a new $Q^{n}$. The problem (14) can be solved again, and we can recalculate $U^{n}_{2}$, $U^{n}$ and $Q^{n}$ again and again until some termination condition is satisfied. The algorithm is shown in Algorithm 1.

The interval $[0,T]$ has been divided into $N$ pieces. Choose $\Delta\tilde{S}=\tilde{S}_{max}/I$ and $\Delta v=v_{max}/J$, and let $\tilde{S}_{i}=i\cdot\Delta\tilde{S},\ i=0,...,I,\ v_{j}=j\cdot\Delta v,\ j=0,...,J$, then we get a uniform grid. Denote $U_{i,j}^{n}=U(\tilde{S}_{i},v_{j},\tau_{n})$ and $U_{1,i,j}^{n}$, $U_{2,i,j}^{n}$ and $Q_{i,j}^{n}$ are similar.

We approximate $Q^{n}=\mathscr{L}_{2}U^{n}$ with an upwind scheme as follows:

To compute $U^{n}_{2}$, we give an implicit central difference scheme:

It is noteworthy that the iterative process is not always convergent with the algorithm and the finite difference scheme above. We calculate $Q^{n}$ based on the solution $U^{n}$ totally from the last iteration and it may lead to divergence when the coefficients in the expression of $Q^{n}$ (16) get large. A mixed method should be used to calculate $Q^{n}$ to guarantee convergence.

Let

and define $Q_{1}^{n}$, $Q_{2}^{n}$, $Q^{n}_{1,i,j}$, $Q^{n}_{2,i,j}$ for $i=0,...,I,\ j=0,...,J,\ n=0,...,N$ similarly to $U^{n}$. To get $U_{2}^{n}$, a central difference scheme is still used to approximate $\mathscr{L}_{1}$ as in (17). $Q^{n}$ can be divided into two parts including $Q^{n}_{1}$ and $Q^{n}_{2}$. We compute $Q^{n}_{1}=\mathscr{L}_{2}U^{n}_{1}$ with the same discretization of $\mathscr{L}_{2}$ in (16). In fact, $Q^{n}_{1}$ can also be calculated exactly since the expression of $U^{n}_{1}$ is known. In terms of $Q^{n}_{2}$, we still discretize $\mathscr{L}_{2}U_{2}$ as in (16), but instead of computing with $U_{2}^{n}$ totally from the last iteration, we make the terms on the current grid point $(i,j)$ in the discretization come from the current iteration step and keep others unchanged from the last step. We present the improved algorithm in Algorithm 2.

To make it clearer, the following numerical scheme is used to get $U_{2}^{n}$:

Here, $\hat{Q}_{2,i,j}^{n}$ is computed by the following expression but note that $U^{n}$ from the last iteration step should be used:

Boundary conditions (6)-(9) from [20] have been widely used in both academia and industry. For the conditions imposed at infinity, it is common to use them on a chosen boundary to find the solution and then truncate the domain into a smaller one to guarantee accuracy, which may cause large computational cost.

For the boundary condition in the $\tilde{S}$-direction, we use the original boundary condition for $U$ in [20] to give the boundary condition for $U_{2}$ as is shown in (14) (we denote it as OriginalBC). Besides, we derive another two boundary conditions for $U_{2}$ based on the artificial boundary methods in [27] which have been proven to be able to improve the accuracy while lower the cost.

We introduce $\tilde{S}_{max}>1$ as an artificial boundary. Then with the artificial boundary methods in [27], two artificial boundary conditions on $\{\tilde{S}=\tilde{S}_{max}\}$ called MApABC1 and MApABC2 can be derived and presented in the same form:

In MApABC1,

In MApABC2,

In addition, we replace the boundary condition (9) for $U$ on $v\to+\infty$ by the one in [26] for accuracy and then for $U_{2}$, the boundary condition should be of the same form:

We also choose $v_{max}$ to restrict the option pricing problem on $[0,\tilde{S}_{max}]\times[0,v_{max}]\times[0,T]$.

The discrete boundary condition on $\tilde{S}=0$ can be given easily:

For OriginalBC, we have

For MApABC1 and MApABC2, we have

where

For MApABC1, $H_{j}^{n}$ can be calculated approximately with the trapezoidal rule applied to (21) with a discrete $Q$ the same as (16). Also, for a better algorithm, we are not supposed to calculate $Q$ in (21) by $U$ totally from the last iteration step. Instead, we need to divide $Q$ into $Q_{1}$ and $Q_{2}$ and use a mixed method for $Q_{2}$ as we have done to the equation before.

For MApABC2, We use a family of curves

where $\gamma_{0},\gamma_{1},\mu,\sigma$ are all parameters to fit $Q_{j}^{n}(\tilde{S})$ with $Q_{j}^{n}(\tilde{S_{i}})=Q_{i,j}^{n},\ i=1,...,I-1,\ j=1,...,J-1$ approximated in the same way of (16). When the parameters of the curve are determined, the integral in (22) in the $\tilde{S}$-direction can be easily calculated and a trapezoidal rule is also used in the $\tau$-direction to obtain $H_{j}^{n}$. The convergence of the iterative process should be taken into account and dealt with as before.

The boundary conditions in the $v$-direction can be discretized similarly so we ignore the details here.

In this section, we show the results of some numerical experiments for the option pricing problem. We implement the iterative splitting method in Section 3 with the boundary conditions in Section 4.

Let $\Delta\tau=\Delta\tilde{S}=\Delta v:=h$, we focus on the behaviors of the proposed method using different step sizes. We choose $\tilde{S}_{max}=4$ and $v_{max}=4$ as the artificial boundary of the computational domain. The 2nd-order asymptotic solution in [27] is used as the reference solution and we compute the relative errors of the numerical results in $l^{2}$-norm. As for the stopping condition of the iterative process in the iterative splitting method, we set a limitation of $10^{-4}$ on the $l^{2}$-norm of the difference between the $U_{2}$’s calculated in the two adjacent iteration steps.

Table 1 shows the first set of parameters with a large mean reversion rate $\kappa$. The relative errors with respect to the reference solution are presented in Table 2. As the step size gets smaller, the relative error of the solution by the iterative splitting method with MApABC2 decreases and the algorithm achieves almost 1st order convergence. With $h=0.1$ fixed, we show the numerical solutions and the errors with OriginalBC and MApABC2 on $\tilde{S}=3$ and $\tilde{S}=4$ in Figure 1 and Figure 2 respectively. The numerical solutions and the errors on $v=3$ and $v=4$ are presented in Figure 3 and Figure 4 as well. Compared with OriginalBC, MApABC2 gives much better solutions.

The second set of parameters is shown in Table 3. We compare the iterative splitting method proposed in Section 3 with a classic 2-dimensional finite difference method as in [27]. Table 4 shows the relative errors of the solutions corresponding to the two methods and boundary conditions with respect to the reference solution.In this example, the mean reversion rate $\kappa$ and the volatility of the variance $\sigma$ are both relatively small, so that the dynamic of the volatility $v(t)$ should be close to a constant. With the iterative splitting method, the first part of the solution $U_{1}$ accounts for most of the option price $U$ and since we have got the exact solution of $U_{1}$, a small step size is enough to give a solution of appreciable accuracy. As is shown in Table 4, the iterative splitting method has a significant advantage over the 2-dimensional finite difference method.

In Table 5, we present the third set of parameters of moderate sizes. We conduct numerical experiments using both the iterative splitting method and the 2-dimensional finite difference method as we mentioned before and show the relative errors in Table 6. The iterative splitting method with MApABC2 gives the most accurate solution, especially when we set the step size to be $0.05$, the relative error is only $1.6\times 10^{-4}$.

Some of the derivatives of the option price with respect to different variables are called Greeks in finance and are of high importance in industry. For example, the first and second derivatives of the option price with respect to the underlying price are called Delta and Gamma, while Vega corresponds to the first derivative with respect to the variance of the underlying. Greeks are good measures of risk and sometimes determine the trading volume so the computational efficiency and accuracy have received wide attention. With $h=0.1$ fixed, we plot the errors of the option price and Greeks on $\tilde{S}=4$ in Figure 5. The numerical results show that the iterative splitting method with MApABC2 behaves the best in terms of accuracy, and the calculation is simpler and more convenient as well since the iterative splitting method helps to transform a 2-dimensional problem into two quasi 1-dimensional ones.

We propose an iterative splitting method to solve the option pricing problem and concentrate on the Heston stochastic volatility model. We introduce the artificial boundary conditions based on the methods from [27]. Numerical study has shown that for the European options, the iterative splitting method can improve the accuracy by contrast with the classic two-dimensional finite difference scheme. Using the artificial boundary condition called MApABC2, the solution and the Greeks are calculated better in terms of the error. The iterative splitting method transformed a two-dimensional PDE into two quasi one-dimensional PDEs with the idea of operator splitting. The convergence of the iteration is ensured by a mixed method. Our new method is easy to interpret, understand and implement. It can also be extended to other models and options, and then the existed methods for the Black-Scholes model can be available, which will be considered in the future.