# Efficient Simulation for Pricing Barrier Options with Two-Factor Stochastic Volatility and Stochastic Interest Rate

Authors: Zhang Sumei; Zhao Jieqiong
Journal: Mathematical Problems in Engineering
 (2017)

---

## 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.


---

## Body

## 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 [3–9] 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 [14–17], 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. Section2 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{Ω,F,{Ft},P} be a complete probability space, where P is a risk-neutral bank account measure. Assume that WjS(t), WjV(t) (j=1,2), and Wr(t) are all standard Brownian motions which are Ft-adapted. The double Heston model proposed by Christoffersen et al. [2] is defined by the following stochastic partial differential equations: (1)dSt=rStdt+∑j=12VjtStdWjSt,dV1t=b1θ1-V1tdt+σ1V1tdW1Vt,dV2t=b2θ2-V2tdt+σ2V2tdW2Vt,where r is constant interest rate, WjS(t) and WjV(t) (j=1,2) are correlated with parameter ρj. bj, θj, σj are the mean reversion speed, long-term volatility level, and volatility of process Vj(t) (j=1,2), respectively. 2bjθj≥σj2 are supposed to make the processes Vj(t) remain strictly positive. Conditional on Vj(t), Vj(t+Δ) obey σj2(1-e-bjΔ)/4bj times a noncentral χ2 distribution with dj degree of freedom and noncentrality parameter ξj; that is, (2)PVjt+Δ≤xj∣Vjt=Fχdj2ξj4bjσj21-e-bjΔxj,where dj=4bjθj/σj2, ξj=4bje-bjΔVj(t)/σj2(1-e-bjΔ), and Fχdj2(ξj)(zj)=∑i=1∞‍ξjie-ξj/2/2dj/2+2ii!Γ(i+dj/2)∫0zj‍zjdj/2+i-1e-v/2dv.We consider the Hull-White stochastic interest rate [22] which is driven by the following mean-reverting process: (3)drt=θt-brrtdt+σrdWrt,where positive constants br, σr are the mean reversion speed and volatility of process r(t), respectively. θ(t)>0 is used to fit the initial term structure of process r(t). Following Brigo and Mercurio [23], the price of a zero-coupon bond maturing at time T can be formulated as follows: (4)Pt,T=At,Texp⁡-Bt,Trt,where (5)At,T=P0,TP0,texp⁡Bt,Tf0,t-σr24br1-e-2brtBt,T2,Bt,T=1-e-brT-tbr,f0,t=-∂ln⁡P0,t∂t.We replace constantr 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: (6)dSt=rtStdt+∑j=12VjtStdWjSt,(7)dV1t=b1θ1-V1tdt+σ1V1tdW1Vt,(8)dV2t=b2θ2-V2tdt+σ2V2tdW2Vt,(9)drt=θt-brrtdt+σrdWrt.Suppose V1(0)=V1, V2(0)=V2, r(0)=r, and S(0)=S. Assume that any two random processes are uncorrelated with each other except cov(dW1S(t),dW1V(t))=ρ1dt, cov(dW2S(t),dW2V(t))=ρ2dt.


## 3. Simulation Scheme for the DHHW Model

To reduce the dimension of the Monte Carlo simulation, we change from the measureP to the T-forward measure Q by using P(0,T) as numeraire. Set (10)zt,T=ln⁡StPt,T.By the Itô formula, we rewrite (6) as follows: (11)dzt,T=-12V1t+V2t+σr2B2t,T+∑j=12VjtWjSt+σrBt,TWrt.


### 3.1. Variance Simulation

Based on the fact that a noncentralχ2 distribution with high noncentrality parameter can be well approximated by a normal distribution, we use quadratic exponential scheme [24] for discrete variance process Vj(t).By (2) and simple calculation, we have (12)mj=EVjt+Δ∣Vjt=θj+Vjt-θje-bjΔ,sj2=VVjt+Δ∣Vjt=σj2e-bjΔVjtbj1-e-bjΔ+θjσj22bj1-e-bjΔ2.Set ψj=sj2/mj2. Provided that ψj≤2 and ξj→∞, we approximate Vj(t+Δ) (the segment of high value) by (13)Vjt+Δ≈αjβj+ZVj2,where αj=mj/1+bj2, βj2=2ψj-1-1+22ψj-12ψj-1-1, and ZVj are independent standard normal random variables.Provided thatψj≥1 and ξj→0, we approximate Vj(t+Δ) (the segment of low value) by(14)Vjt+Δ≈mjψj+12log⁡21-Uj1+ψj,ψj-1ψj+1<Uj≤10,0≤Uj≤ψj-1ψj+1,where Uj are independent uniform random numbers.


### 3.2. Asset Price Simulation

Under theT-forward measure Q, by integrating (11) and applying the Cholesky decomposition, we rewrite asset price process (6) as follows: (15)zt+Δ,T=zt,T-12∫tt+ΔV1u+V2u+σr2B2u,Tdu+∑j=12ρj∫tt+ΔVjudWjVu+∑j=121-ρj2∫tt+ΔVjudWju+σr∫tt+ΔBu,TdWru,where Wj(u) are Brownian motions independent of WjV(u) and Wr(u).By the drift interpolation method [25], we approximate the integral of the variance process Vj(t) (j=1,2) by (16)∫tt+ΔVjudu≈12Vjt+Vjt+ΔΔ.By (7) and (8), we have (17)∫tt+ΔVjudWjVu=1σjVjt+Δ-Vjt-bjθjΔ+bj∫tt+ΔVjudu.The integral∫tt+Δ‍Vj(u)dWj(u) (j=1,2) obeyed normal distribution with mean zero and variance ∫tt+Δ‍Vj(u)du. By (16), we approximate ∫tt+Δ‍Vj(u)dWj(u) by (18)∫tt+ΔVjudWju≈12Vjt+Vjt+ΔΔZj,where Zj are independent standard normal random variables.By direct calculation, we have(19)∫tt+ΔB2u,Tdu=1br2Δ+e-2brT-t-Δ-e-2brT-t2br+2e-brT-t-2e-brT-t-Δbr=Ct,t+Δ.Accordingly, we have the integral ∫tt+Δ‍B(u,T)dWr(u) obeying normal distribution with mean zero and variance ∫tt+Δ‍B2(u,T)du; that is, (20)∫tt+ΔBu,TdWru=Ct,t+ΔZr,where Zr is a standard normal random variable.Substituting (16)–(20) into (15), we have the simulation scheme for the asset price process as follows: (21)zt+Δ,T=zt,T+K0+∑j=12K1jVjt+K2jVjt+Δ+K3jZj+K4Zr,where (22)K0=-σr22br2Δ+e-2brT-t-Δ2br-2e-brT-t-Δbr-e-2brT-t2br+2e-brT-tbr-bjθjρjΔσj,K1j=-Δ4+bjρjΔ2σj-ρjσj,K2j=-Δ4+bjρjΔ2σj+ρjσj,K3j=Δ21-ρj2Vjt+Vjt+Δ,K4=σr2br2Δ+e-2brT-t-Δ2br-2e-brT-t-Δbr-e-2brT-t2br+2e-brT-tbr.


## 3.1. Variance Simulation

Based on the fact that a noncentralχ2 distribution with high noncentrality parameter can be well approximated by a normal distribution, we use quadratic exponential scheme [24] for discrete variance process Vj(t).By (2) and simple calculation, we have (12)mj=EVjt+Δ∣Vjt=θj+Vjt-θje-bjΔ,sj2=VVjt+Δ∣Vjt=σj2e-bjΔVjtbj1-e-bjΔ+θjσj22bj1-e-bjΔ2.Set ψj=sj2/mj2. Provided that ψj≤2 and ξj→∞, we approximate Vj(t+Δ) (the segment of high value) by (13)Vjt+Δ≈αjβj+ZVj2,where αj=mj/1+bj2, βj2=2ψj-1-1+22ψj-12ψj-1-1, and ZVj are independent standard normal random variables.Provided thatψj≥1 and ξj→0, we approximate Vj(t+Δ) (the segment of low value) by(14)Vjt+Δ≈mjψj+12log⁡21-Uj1+ψj,ψj-1ψj+1<Uj≤10,0≤Uj≤ψj-1ψj+1,where Uj are independent uniform random numbers.


## 3.2. Asset Price Simulation

Under theT-forward measure Q, by integrating (11) and applying the Cholesky decomposition, we rewrite asset price process (6) as follows: (15)zt+Δ,T=zt,T-12∫tt+ΔV1u+V2u+σr2B2u,Tdu+∑j=12ρj∫tt+ΔVjudWjVu+∑j=121-ρj2∫tt+ΔVjudWju+σr∫tt+ΔBu,TdWru,where Wj(u) are Brownian motions independent of WjV(u) and Wr(u).By the drift interpolation method [25], we approximate the integral of the variance process Vj(t) (j=1,2) by (16)∫tt+ΔVjudu≈12Vjt+Vjt+ΔΔ.By (7) and (8), we have (17)∫tt+ΔVjudWjVu=1σjVjt+Δ-Vjt-bjθjΔ+bj∫tt+ΔVjudu.The integral∫tt+Δ‍Vj(u)dWj(u) (j=1,2) obeyed normal distribution with mean zero and variance ∫tt+Δ‍Vj(u)du. By (16), we approximate ∫tt+Δ‍Vj(u)dWj(u) by (18)∫tt+ΔVjudWju≈12Vjt+Vjt+ΔΔZj,where Zj are independent standard normal random variables.By direct calculation, we have(19)∫tt+ΔB2u,Tdu=1br2Δ+e-2brT-t-Δ-e-2brT-t2br+2e-brT-t-2e-brT-t-Δbr=Ct,t+Δ.Accordingly, we have the integral ∫tt+Δ‍B(u,T)dWr(u) obeying normal distribution with mean zero and variance ∫tt+Δ‍B2(u,T)du; that is, (20)∫tt+ΔBu,TdWru=Ct,t+ΔZr,where Zr is a standard normal random variable.Substituting (16)–(20) into (15), we have the simulation scheme for the asset price process as follows: (21)zt+Δ,T=zt,T+K0+∑j=12K1jVjt+K2jVjt+Δ+K3jZj+K4Zr,where (22)K0=-σr22br2Δ+e-2brT-t-Δ2br-2e-brT-t-Δbr-e-2brT-t2br+2e-brT-tbr-bjθjρjΔσj,K1j=-Δ4+bjρjΔ2σj-ρjσj,K2j=-Δ4+bjρjΔ2σj+ρjσj,K3j=Δ21-ρj2Vjt+Vjt+Δ,K4=σr2br2Δ+e-2brT-t-Δ2br-2e-brT-t-Δbr-e-2brT-t2br+2e-brT-tbr.


## 4. Simulation Algorithm for Pricing Barrier Options

Under theT-forward measure Q, we evaluate barrier options by the following formula: (23)Vbt,T=Pt,TEQΛSτB,τB,where Λ(S(τB),τB) is a payoff function, B is barrier level, and τB is the first time when barrier is hit. For a down and out call (DOC) option τB is defined as follows: (24)τB=inf⁡0≤t≤T:St>B.Λ(S(τB),τB) is given by (25)ΛSτB,τB=max⁡ST-K,0,τB=T0,τB<T,where K is exercise price at maturity time T.


### 4.1. The Basic Algorithm

Based on the simulation scheme for the DHHW model, we evaluate a DOC option by Algorithm1.Algorithm 1:The basic algorithm for pricing a DOC option.
(
1
)Initial K, S, B, T, r0, br, σr, bj, θj, σj, ρj, Vj, ψc;
(
2
)Compute P(0,T) using (4);
(
3
)
s
u
m
_
p
=
0;
(
4
)    repeatfor i=1 to M
(
5
)
Compute mj,sj2,ψj,αj,βj for j=1,2;
(
6
)
Draw two uniform random numbers Uj;
(
7
)
Draw two standard normal random numbers ZVj by Φ-1(Uj);
(
8
)
if
ψ
j
≤
ψ
c
then compute Vj(t) by αj(βj+ZVj)2;
(
9
)
else compute Vj(t) by (14);
(
10
)
end if
(
11
)
Compute K0, K1j, K2j, K3j, K4;
(
12
)
Compute z(t) by (21);
(
13
)
Set S(t)=P(t,T)ez(t);
(
14
)
if
S
≤
B, then  payoff = 0;
(
15
)
else  payoff = max⁡(S-K,0);
(
16
)
end if
(
17
)
s
u
m
_
p
=
s
u
m
_
p
+
payoff;
(
18
)
until the last simulation
(
19
)Compute VDOC(0) by VDOC(0)=P(0,T)sum_p/M.


### 4.2. The Hybrid Algorithm Based on Variance Reduction Techniques

We combine antithetic variate and control variate techniques to improve the efficiency of the Algorithm1. We use a European option price as control variate and modify step (19) of Algorithm 1 as follows: 19~V^DOC0=P0,Tsum_pM+1M∑j=1MVE0-P0,Tmax⁡ST-K,0,where VE(0) 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, VE(0) can be approximated by the following formula: (26)VE0≈exp⁡-αkπ∑j=1Ne-i2π/Nj-1i-1eibvjψvjη33+-1j-δj-1,where k=log⁡K, α is damping factor, i is an imaginary unit, b=Nh/2, η=2π/Nh (h is a regular spacing), δn is the Kronecker delta function that is unity for n=0 and zero otherwise, and ψ(v)=φ(v-(α+1)i)/α2+α-v2+i(2α+1)v, where φ(u)=EP[e-∫tT‍rudu+iuln⁡S(t)] is the discounted T-forward characteristic function of log⁡S(t).Based on the results of Grzelak et al. [14] and simple calculation, we have the discounted T-forward characteristic function of the DHHW model as follows: (27)φu=eiuln⁡St+Du,τ+∑j=12Eju,τVjt+Fu,τrt,where (28)Du,τ=iu-1-Bt,Tγt+ln⁡P0,tP0,T+12H0,T-H0,t-12u2-1+2iuHt,T+bjθjσj2bj-ρjσjiu-Djτ-2ln⁡1-Gje-Djτ1-Gj,Eju,τ=bj-ρjσjiu-Djσj21-e-Djτ1-Gje-Djτ,Fu,τ=1briu-11-e-brτ,Ht,T=σr2br2T-t+2bre-brT-t-12bre-2brT-t-32br,γt=f0,t+σr21-e-brt22br2,Gj=bj-ρjσjiu-Djbj-ρjσjiu+Dj,Dj=bj-ρjσjiu2+σj2u2+iu,τ=T-t.For computingΦ-1(Uj) in step (7) 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 ZVj, we consider antithetic variable -ZVj and compute Vj(t) by αj(βj-ZVj)2 and thus obtain the second path of the DHHW model. By combining control variate technique we modify the last step of Algorithm 1 as (29)V^DOC0=P0,T0.5sum_p1+sum_p2M+1M∑j=1MVE0-P0,Tmax⁡ST-K,0.


## 4.1. The Basic Algorithm

Based on the simulation scheme for the DHHW model, we evaluate a DOC option by Algorithm1.Algorithm 1:The basic algorithm for pricing a DOC option.
(
1
)Initial K, S, B, T, r0, br, σr, bj, θj, σj, ρj, Vj, ψc;
(
2
)Compute P(0,T) using (4);
(
3
)
s
u
m
_
p
=
0;
(
4
)    repeatfor i=1 to M
(
5
)
Compute mj,sj2,ψj,αj,βj for j=1,2;
(
6
)
Draw two uniform random numbers Uj;
(
7
)
Draw two standard normal random numbers ZVj by Φ-1(Uj);
(
8
)
if
ψ
j
≤
ψ
c
then compute Vj(t) by αj(βj+ZVj)2;
(
9
)
else compute Vj(t) by (14);
(
10
)
end if
(
11
)
Compute K0, K1j, K2j, K3j, K4;
(
12
)
Compute z(t) by (21);
(
13
)
Set S(t)=P(t,T)ez(t);
(
14
)
if
S
≤
B, then  payoff = 0;
(
15
)
else  payoff = max⁡(S-K,0);
(
16
)
end if
(
17
)
s
u
m
_
p
=
s
u
m
_
p
+
payoff;
(
18
)
until the last simulation
(
19
)Compute VDOC(0) by VDOC(0)=P(0,T)sum_p/M.


## 4.2. The Hybrid Algorithm Based on Variance Reduction Techniques

We combine antithetic variate and control variate techniques to improve the efficiency of the Algorithm1. We use a European option price as control variate and modify step (19) of Algorithm 1 as follows: 19~V^DOC0=P0,Tsum_pM+1M∑j=1MVE0-P0,Tmax⁡ST-K,0,where VE(0) 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, VE(0) can be approximated by the following formula: (26)VE0≈exp⁡-αkπ∑j=1Ne-i2π/Nj-1i-1eibvjψvjη33+-1j-δj-1,where k=log⁡K, α is damping factor, i is an imaginary unit, b=Nh/2, η=2π/Nh (h is a regular spacing), δn is the Kronecker delta function that is unity for n=0 and zero otherwise, and ψ(v)=φ(v-(α+1)i)/α2+α-v2+i(2α+1)v, where φ(u)=EP[e-∫tT‍rudu+iuln⁡S(t)] is the discounted T-forward characteristic function of log⁡S(t).Based on the results of Grzelak et al. [14] and simple calculation, we have the discounted T-forward characteristic function of the DHHW model as follows: (27)φu=eiuln⁡St+Du,τ+∑j=12Eju,τVjt+Fu,τrt,where (28)Du,τ=iu-1-Bt,Tγt+ln⁡P0,tP0,T+12H0,T-H0,t-12u2-1+2iuHt,T+bjθjσj2bj-ρjσjiu-Djτ-2ln⁡1-Gje-Djτ1-Gj,Eju,τ=bj-ρjσjiu-Djσj21-e-Djτ1-Gje-Djτ,Fu,τ=1briu-11-e-brτ,Ht,T=σr2br2T-t+2bre-brT-t-12bre-2brT-t-32br,γt=f0,t+σr21-e-brt22br2,Gj=bj-ρjσjiu-Djbj-ρjσjiu+Dj,Dj=bj-ρjσjiu2+σj2u2+iu,τ=T-t.For computingΦ-1(Uj) in step (7) 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 ZVj, we consider antithetic variable -ZVj and compute Vj(t) by αj(βj-ZVj)2 and thus obtain the second path of the DHHW model. By combining control variate technique we modify the last step of Algorithm 1 as (29)V^DOC0=P0,T0.5sum_p1+sum_p2M+1M∑j=1MVE0-P0,Tmax⁡ST-K,0.


## 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):(30)RMSE=1NT×NK∑t=1NT∑k=1NKVtkΘ-Vtk2Vtk,where VtkΘ and Vtk are the model and market prices, respectively. We calibrate the DHHW model by solving the following nonlinear least-squares optimization problem: (31)Θ∗=arg⁡min⁡RMSE.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.Table 1
Estimated parameter values for the DHHW model calibrated to S&P 500 index call options on March 28, 2017.
Parameter values
RMSE
br=0.201
σr=0.021
r0=0.012
b1=0.101
θ1=0.019
σ1=0.189
V
1
=
0.011
ρ1=-0.251
0.0679
b2=2.114
θ2=0.012
σ2=0.253
V
2
=
0.021
ρ
2
=
-
0.823


### 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 N=4096, h=π/300=0.01047, and α=1.25 and set S=100, K=100, and T=2. For our simulation scheme, we specify ψc=1.5 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
Comparison of the accuracy of our simulation scheme and Euler scheme for pricing European options. The exact priceV=13.209.
Number of simulations
Number of time steps
Our simulation scheme
Euler scheme
Standard error
Relative error
Standard error
Relative error
1000
100
0.4780
0.0421
0.4681
0.1182
1000
500
0.4559
0.0424
0.4800
0.0472
10000
100
0.1515
0.0166
0.1533
0.0227
10000
500
0.1566
0.0032
0.1561
0.0166
100000
100
0.0480
0.0267
0.0479
0.0523
100000
500
0.0480
0.0257
0.0486
0.0334Table2 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 barrierB=90 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.Table 3
Comparison of the standard errors among the different Monte Carlo techniques for pricing DOC options, where time stepsN=100 and barrier B=90.
K
M
Basic algorithm
Control variates
Antithetic variates
Hybrid algorithm
80
1000
0.6548
0.3167
0.4248
0.1047
10000
0.2079
0.1012
0.1334
0.0342
100000
0.0657
0.0312
0.0428
0.0114
100
1000
0.4671
0.2317
0.3323
0.0112
10000
0.1502
0.0689
0.1098
0.0248
100000
0.0478
0.0219
0.0355
0.0077
120
1000
0.2843
0.1418
0.2454
0.0456
10000
0.0929
0.0437
0.0850
0.0152
100000
0.0298
0.0136
0.0276
0.0049Furthermore, based on the hybrid algorithm, we plot the implied volatility surface of the DHHW model. Figure1 examines the effects of interest rate parameters br, σr and the second variance parameters b2, θ2, σ2, ρ2 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
Implied volatility surface under the DHHW model with different values of parametersbr, σr, b2, θ2, σ2, and ρ2.Figure1 summarizes our findings. Both br and σr only affect the long end regarding the maturity of the volatility surface. However, decreasing br leads to a significant rise of the slope of the implied volatility smile, while the contrary occurs by decreasing σr. Increasing b2 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 θ2 results in the rise of overall volatility, but the slope of the smile hardly changed. The curvature of the smile slightly increases by increasing σ2. By changing ρ2, 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.


## 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):(30)RMSE=1NT×NK∑t=1NT∑k=1NKVtkΘ-Vtk2Vtk,where VtkΘ and Vtk are the model and market prices, respectively. We calibrate the DHHW model by solving the following nonlinear least-squares optimization problem: (31)Θ∗=arg⁡min⁡RMSE.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.Table 1
Estimated parameter values for the DHHW model calibrated to S&P 500 index call options on March 28, 2017.
Parameter values
RMSE
br=0.201
σr=0.021
r0=0.012
b1=0.101
θ1=0.019
σ1=0.189
V
1
=
0.011
ρ1=-0.251
0.0679
b2=2.114
θ2=0.012
σ2=0.253
V
2
=
0.021
ρ
2
=
-
0.823


## 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 N=4096, h=π/300=0.01047, and α=1.25 and set S=100, K=100, and T=2. For our simulation scheme, we specify ψc=1.5 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
Comparison of the accuracy of our simulation scheme and Euler scheme for pricing European options. The exact priceV=13.209.
Number of simulations
Number of time steps
Our simulation scheme
Euler scheme
Standard error
Relative error
Standard error
Relative error
1000
100
0.4780
0.0421
0.4681
0.1182
1000
500
0.4559
0.0424
0.4800
0.0472
10000
100
0.1515
0.0166
0.1533
0.0227
10000
500
0.1566
0.0032
0.1561
0.0166
100000
100
0.0480
0.0267
0.0479
0.0523
100000
500
0.0480
0.0257
0.0486
0.0334Table2 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 barrierB=90 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.Table 3
Comparison of the standard errors among the different Monte Carlo techniques for pricing DOC options, where time stepsN=100 and barrier B=90.
K
M
Basic algorithm
Control variates
Antithetic variates
Hybrid algorithm
80
1000
0.6548
0.3167
0.4248
0.1047
10000
0.2079
0.1012
0.1334
0.0342
100000
0.0657
0.0312
0.0428
0.0114
100
1000
0.4671
0.2317
0.3323
0.0112
10000
0.1502
0.0689
0.1098
0.0248
100000
0.0478
0.0219
0.0355
0.0077
120
1000
0.2843
0.1418
0.2454
0.0456
10000
0.0929
0.0437
0.0850
0.0152
100000
0.0298
0.0136
0.0276
0.0049Furthermore, based on the hybrid algorithm, we plot the implied volatility surface of the DHHW model. Figure1 examines the effects of interest rate parameters br, σr and the second variance parameters b2, θ2, σ2, ρ2 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
Implied volatility surface under the DHHW model with different values of parametersbr, σr, b2, θ2, σ2, and ρ2.Figure1 summarizes our findings. Both br and σr only affect the long end regarding the maturity of the volatility surface. However, decreasing br leads to a significant rise of the slope of the implied volatility smile, while the contrary occurs by decreasing σr. Increasing b2 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 θ2 results in the rise of overall volatility, but the slope of the smile hardly changed. The curvature of the smile slightly increases by increasing σ2. By changing ρ2, 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 theT-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.



---

*Source: 3912036-2017-11-14.xml*