Abstract

This paper presents an extension of the double Heston stochastic volatility model by combining Hull-White stochastic interest rates. By the change of numeraire and quadratic exponential scheme, this paper develops a new simulation scheme for the extended model. By combining control variates and antithetic variates, this paper provides an efficient Monte Carlo simulation algorithm for pricing barrier options. Based on the differential evolution algorithm the extended model is calibrated to S&P 500 index options to obtain the model parameter values. Numerical results show that the proposed simulation scheme outperforms the Euler scheme, the proposed simulation algorithm is efficient for pricing barrier options, and the extended model is flexible to fit the implied volatility surface.

1. Introduction

A barrier option is a path-dependent option which is exterminated (knocked out) or initiated (knocked in) if the underlying spot price hits the specified barrier level during the life of the option. Because of this supplementary risk, barrier options are cheaper than plain vanilla options and thus are widely traded in exchanges worldwide. One-factor stochastic volatility models can generate “smile," leverage effects, and term structure effects which cannot be explained by the Black-Scholes model [1, 2]. Consequently many papers [39] evaluate barrier options under one-factor stochastic volatility models.

However, one-factor stochastic volatility models including the Heston model [10] present the poor performance when fitting the stiff volatility skews [11]. One extension by using multiple stochastic volatility factors has been presented in some literatures [2, 12, 13]. Christoffersen et al. [2] confirm that the double Heston model significantly improves the flexibility of the one-factor stochastic volatility model in capturing the volatility term structure. In addition, stochastic interest rate is crucial for option pricing because it ensures proper discounting of future payoffs. In recent literatures [1417], Hull-White stochastic interest rate which is analytically tractable has been incorporated into one-factor stochastic volatility model for pricing path-dependent options. Therefore, the model which incorporates multifactor stochastic volatility and stochastic interest rate may be more reasonable for pricing barrier options.

Barrier options with less stochastic factors can be efficiently evaluated by partial differential equation (PDE) methods [3, 4, 18, 19]. However, the evaluation of barrier options with multiple stochastic factors is to solve a high-dimensional PDE which makes PDE methods quite complex and potentially prone to accuracy and stability problems [1]. A more efficient method for pricing barrier options in this case is the Monte Carlo method; see [20, 21]. The main purpose of this paper is to provide a Monte Carlo method for pricing barrier options under a two-factor stochastic volatility and stochastic interest rate model.

The main contributions of this paper are threefold. Firstly, this paper extends the double Heston model to stochastic interest rate. Secondly, this paper provides a new simulation scheme for the extended model. Thirdly, the paper develops an efficient Monte Carlo algorithm for pricing barrier options. The rest of the paper is organized as follows. Section 2 presents the extended model. Section 3 details the simulation scheme for the proposed model. Section 4 develops the simulation algorithm for pricing barrier options. Section 5 provides numerical experiments. Section 6 concludes.

2. The Pricing Model

Let be a complete probability space, where is a risk-neutral bank account measure. Assume that , (), and are all standard Brownian motions which are -adapted. The double Heston model proposed by Christoffersen et al. [2] is defined by the following stochastic partial differential equations: where is constant interest rate, and () are correlated with parameter . , , are the mean reversion speed, long-term volatility level, and volatility of process (), respectively. are supposed to make the processes remain strictly positive. Conditional on , obey times a noncentral distribution with degree of freedom and noncentrality parameter ; that is, where , , and .

We consider the Hull-White stochastic interest rate [22] which is driven by the following mean-reverting process: where positive constants , are the mean reversion speed and volatility of process , respectively. is used to fit the initial term structure of process . Following Brigo and Mercurio [23], the price of a zero-coupon bond maturing at time can be formulated as follows: where

We replace constant in (1) with Hull-White stochastic interest rate (3) and define the double Heston Hull-White (DHHW) model by a four-dimensional system of stochastic differential equations: Suppose , , , and . Assume that any two random processes are uncorrelated with each other except , .

3. Simulation Scheme for the DHHW Model

To reduce the dimension of the Monte Carlo simulation, we change from the measure to the -forward measure by using as numeraire. Set By the Itô formula, we rewrite (6) as follows:

3.1. Variance Simulation

Based on the fact that a noncentral distribution with high noncentrality parameter can be well approximated by a normal distribution, we use quadratic exponential scheme [24] for discrete variance process .

By (2) and simple calculation, we have Set . Provided that and , we approximate (the segment of high value) by where , , and are independent standard normal random variables.

Provided that and , we approximate (the segment of low value) bywhere are independent uniform random numbers.

3.2. Asset Price Simulation

Under the -forward measure , by integrating (11) and applying the Cholesky decomposition, we rewrite asset price process (6) as follows: where are Brownian motions independent of and .

By the drift interpolation method [25], we approximate the integral of the variance process () by

By (7) and (8), we have

The integral () obeyed normal distribution with mean zero and variance . By (16), we approximate by where are independent standard normal random variables.

By direct calculation, we have Accordingly, we have the integral obeying normal distribution with mean zero and variance ; that is, where is a standard normal random variable.

Substituting (16)–(20) into (15), we have the simulation scheme for the asset price process as follows: where

4. Simulation Algorithm for Pricing Barrier Options

Under the -forward measure , we evaluate barrier options by the following formula: where is a payoff function, is barrier level, and is the first time when barrier is hit. For a down and out call (DOC) option is defined as follows: is given by where is exercise price at maturity time .

4.1. The Basic Algorithm

Based on the simulation scheme for the DHHW model, we evaluate a DOC option by Algorithm 1.

Initial , , , , , , , , , , , , ;
Compute using (4);
;
    repeatfor to
Compute for ;
Draw two uniform random numbers ;
Draw two standard normal random numbers by ;
ifthen compute by ;
else compute by (14);
end if
Compute , , , , ;
Compute by (21);
Set ;
if, then  payoff = 0;
else  payoff = ;
end if
;
until the last simulation
Compute by .
4.2. The Hybrid Algorithm Based on Variance Reduction Techniques

We combine antithetic variate and control variate techniques to improve the efficiency of the Algorithm 1. We use a European option price as control variate and modify step of Algorithm 1 as follows: where is the exact price of the European call option which can be efficiently calculated by the fast Fourier transform method [26, 27]. Given the characteristic function, can be approximated by the following formula: where , is damping factor, is an imaginary unit, , ( is a regular spacing), is the Kronecker delta function that is unity for and zero otherwise, and , where is the discounted -forward characteristic function of .

Based on the results of Grzelak et al. [14] and simple calculation, we have the discounted -forward characteristic function of the DHHW model as follows: where

For computing in step of Algorithm 1, we use the approximation algorithm of Wichura [28] not by the Matlab function norminv.m. This can greatly decrease the simulation time. With the obtained , we consider antithetic variable and compute by and thus obtain the second path of the DHHW model. By combining control variate technique we modify the last step of Algorithm 1 as

5. Numerical Experiments

5.1. Calibration of the Model

Before applying the DHHW model to option pricing, we need to estimate the model parameters. For this purpose we use financial market data to calibrate the DHHW model. We define the following relative mean squares error (RMSE): where and are the model and market prices, respectively. We calibrate the DHHW model by solving the following nonlinear least-squares optimization problem:

We choose the S&P 500 index call options on March 28, 2017, available online at http://www.cboe.com/. The options have maturities between 366 days and 997 days and strike prices ranging from 2225 to 2500. The closing price of underlying is 2358.57. For the sake of simplicity, we set the dividend yield to zero. We calculate the model prices by (26). We use differential evolution algorithm (global optimum) to seek the optimal parameter set . Table 1 lists the calibrated parameter values and associated RMSE.

5.2. Implementation of the Simulation Algorithm

With the obtained parameter values we first test the performance of the simulation scheme. We apply the simulation scheme to European options and take the fast Fourier transform solutions as benchmark. In (26) we use , , and and set , , and . For our simulation scheme, we specify throughout the paper. We use 1000, 10000, and 100000 simulation trials, respectively. For each simulation trial we use 100 and 500 time steps, respectively. Table 2 compares our simulation scheme with the standard Euler scheme through sample standard errors and relative errors.

Table 2 shows that although there is almost no difference in the sample standard error between the two schemes, significant relative error remains in Euler scheme compared to our simulation scheme. Moreover, it is observed that standard error of our simulation scheme decreases with the increase of the number of simulations and is almost unaffected by time steps. Table 2 implies that our simulation scheme outperforms the Euler scheme.

With the same parameter setting we evaluate the DOC options with barrier by the hybrid algorithm. We still use 1000, 10000, and 100000 simulation trials. Since our simulation scheme is almost unaffected by time steps we only use 100 time steps. Table 3 compares the hybrid algorithm with basic algorithm, control variates, and antithetic variates techniques. The criterion used is sample standard errors produced by the above methods. Table 3 shows that both control variates and antithetic variates can reduce standard errors. The hybrid algorithm which combines the two techniques significantly reduces the standard errors. Table 3 implies that the hybrid algorithm is efficient for pricing barrier options.

Furthermore, based on the hybrid algorithm, we plot the implied volatility surface of the DHHW model. Figure 1 examines the effects of interest rate parameters , and the second variance parameters , , , on the implied volatility surface. To examine the effects of the above parameters we change each parameter by setting three different values while fixing all other parameter values.

Figure 1 summarizes our findings. Both and only affect the long end regarding the maturity of the volatility surface. However, decreasing leads to a significant rise of the slope of the implied volatility smile, while the contrary occurs by decreasing . Increasing results in a fall of volatility, and this effect on the short-term volatility is more significant, which leads to the slight decrease of the curvature of the smile. Increasing results in the rise of overall volatility, but the slope of the smile hardly changed. The curvature of the smile slightly increases by increasing . By changing , the slope of the smile significantly changes and a negative slope is displayed in the implied volatility. Figure 1 implies that interest rates and the second variance process have important effects on the implied volatility surface of the DHHW model.

6. Conclusion

We propose the DHHW model by combining the double Heston stochastic volatility and Hull-White stochastic interest rate. Under the -forward measure, this paper develops a simulation scheme for the DHHW model. Combining the control variates and antithetic variates, this paper provides a hybrid Monte Carlo algorithm for pricing barrier options under the DHHW model. Numerical results show that the proposed scheme outperforms Euler scheme and the hybrid algorithm is efficient and easy to implement in pricing barrier options. Extensive implied volatility experiments show that implied volatility surfaces present different shapes by varying the parameter values of interest rate and the second variance process, which verify that the DHHW model is flexible to fit the implied volatility smile.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors would like to acknowledge the support from the National Natural Science Foundation of China (Grant no. 11601420), the Natural Science Basic Research Plan in Shaanxi Province of China (Grant no. 2017JM1021), and the Scientific Research Program Funded by Shaanxi Provincial Education Department (Grant no. 17JK0714).