Efficient Simulations of Option Pricing and Greeks Under Three-Factor Model by Conditional Monte Carlo Method

Shanshan CHEN, Chenglong XU, Zhaokui SHI

Journal of Systems Science and Information ›› 2023, Vol. 11 ›› Issue (3) : 314-331.

PDF(184 KB)
PDF(184 KB)
Journal of Systems Science and Information ›› 2023, Vol. 11 ›› Issue (3) : 314-331. DOI: 10.21078/JSSI-E2022053
 

Efficient Simulations of Option Pricing and Greeks Under Three-Factor Model by Conditional Monte Carlo Method

Author information +
History +

Abstract

This paper proposes a hybrid Monte Carlo simulation method for pricing European options under the stochastic volatility model and three-factor model. First, the European options are expressed as a conditional expectation formula, which can be used not only for reducing variance of simulations, but also for calculating the value of Greeks easily, due to the elimination of the weak singularity for the payoff of the option. Then, in order to reduce variance further, the authors also construct a new explicit regression based control variate under Heston model and three-factor model respectively. Numerical results of experiments show that the proposed method can greatly reduce the variance of simulation for pricing European option, and is easy to complement for the calculation of Greeks.

Key words

conditional Monte Carlo / control variate / stochastic volatility / stochastic interest rate / option pricing

Cite this article

Download Citations
Shanshan CHEN , Chenglong XU , Zhaokui SHI. Efficient Simulations of Option Pricing and Greeks Under Three-Factor Model by Conditional Monte Carlo Method. Journal of Systems Science and Information, 2023, 11(3): 314-331 https://doi.org/10.21078/JSSI-E2022053

1 Introduction

With the growth of uncertainty in financial markets, derivative pricing models are becoming more and more complicated, and pricing these securities poses a challenging task. In addition to a small part which can get an analytical solution by the probability theory or partial differential equation theory, the vast majority of the products have no closed-form pricing formula, so we need to turn to the numerical method.
There are three kinds of typical numerical method for pricing financial derivatives: Binomial tree method, finite difference method and Monte Carlo method. The methods of Binomial tree and the finite difference are efficient for low-dimensional problems, but don't work for high dimensional problems, because the storage and calculation increase exponentially with the dimensionality. The method of Monte Carlo simulation for pricing option is a popular approach, since efficient algorithms have been developed for optimal stopping problems[1]. The advantage of Monte Carlo simulation method is that it is less sensitive to the dimensionality of the pricing problem and suitable for parallel computation. The main disadvantage is that the rate of convergence is only half order. So some methods for variance reduction, such as the control variables, dual variables, stratified sampling, importance sampling technique and conditional Monte Carlo (CMC) technology in the implementation have been proposed.
The deviation of Monte Carlo simulation consists of two parts: The statistical error and discretization error generated in the simulation process[2]. The statistical error due to the finite sample and different sampling methods can be improved by variance reduction techniques. The discretization error derived from the sample path simulation mainly be reduced by the sample path with higher order convergence.
In the recent work, Kemna and Vorst[3] used the geometric average Asian option as a controlling to price arithmetic average Asian options, achieved great variance reduction effect. Longstaff and Schwartz[4] used the dual variable variance technique to reduce the pricing error of American options. Fouque and Tullie[5] proposed an optimal importance sampling method based on martingale representation theorem to price option with stochastic volatility model. Bolia, Junejia and Glasserman[6] applied the importance sampling method based on function approximation for American option pricing. Broadie and Kaya[7] employed conditional Monte Carlo method to accelerate the option pricing with stochastic volatility model and jump diffusion process in numerical calculation. The variance reduction technique is one of the most important acceleration method, and its popularity rests on the ease of implication, the availability of controls, and the straightforward intuition of the underlying theory. The method of control variate exploits information about the errors in estimates of known quantities, and it is used to reduce the error in an estimate of an unknown quantity.
The aim of this paper is to propose some methods of acceleration simulation for pricing European options: Conditional Monte Carlo method and method of conditional Monte Carlo with control variate, which can been used not only to reduce the variance of simulation, but also can be used for calculation of Greeks, which is almost as important as option pricing.
The innovation of this paper is that we combine conditional Monte Carlo and regression-based control variable methods to accelerate the process of Monte Carlo simulation, and synthesize on the basis of traditional variance reduction techniques to obtain a better variance reduction technique. Based on Taylor expansion, a better random variable is successfully constructed as a control variable, and a better acceleration effect is achieved. We apply the conditional Monte Carlo acceleration algorithm combined with the control variable to the Heston stochastic volatility model and three-factor model, the first- and second-order Greeks calculation of the option price with respect to the initial stock price is carried out. Numerical examples verify the effectiveness of the accelerated method.
The idea of conditional Monte Carlo with control variate may also applicable for pricing options written on several assets, such as the Basket option, the Asian option with discrete observing points and the Barrier option with discrete observing points and other path dependent options.
The rest of the paper is organized as follows. We first give some models of asset in Section 2. In Section 3, we propose the method of conditional Monte Carlo for reducing variance of simulations for pricing European option, under the Heston model[8] and the three-factor model[9] respectively. Then for reducing variance of simulation further, we construct a hybrid acceleration method which combining the method of conditional Monte Carlo and the method of control variate. Numerical results show that the hybrid method has great acceleration effect. In Section 4, we give an application of the calculation of Greeks for conditional Monte Carlo method, which is important for the sensitivity analysis of option. Section 5 gives the numerical results of our methods, which coincide with the theoretical analysis well.

2 Heston Model and Three-Factor Model of Asset

2.1 Heston Stochastic Volatility Model

In this section, we describe how the movement of the underlying asset's price is defined. In the Black-Scholes model, the volatility term is a constant, this assumption is not always satisfied by real-life options as the probability distribution of a stock price has a fatter left tail and a thinner right tail than the log-normal distribution in [10], and the assumption of constant volatility in financial model is incompatible with stock prices observed in the market, verified by volatility smile. The concept of stochastic volatility was introduced by Hull and White[11] and subsequent developments include the work of Scott[12], Stein and Stein[13], Ball and Roma[14], and Heston[8]. They proposed different stochastic volatility models. The models of underlying asset's price is usually given by the Heston model and the three-factor model. We recall that on a filtered probability space (Ω,F,(Ft)t0,P) under the assumption that P is the risk-neutral pricing measure, the Heston stochastic volatility model is given by the system of stochastic differential equations
{ dS(t)=rS(t)dt+Y(t)S(t)dW(t),dY(t)=α(θY(t))dt+σY(t)dZ(t),
(1)
where W(t) and Z(t) are standard Brown motions under P with dW(t),Z(t)=ρdt, α is the mean-reversion speed, θ is the long-term average variance, and σ is the volatility of Y(t). The first SDE in (1) is responsible for the stock price process S(t) under the risk neutral measure and the second one for the instantaneous variance Y(t). The variance process is known as square root diffusion in the literature. The SDE of Y(t) was first proposed for the short rate modeling in [15] and was called the CIR model in the term structure literature. Y(t) never reaches zero if 2αθ>σ2 holds. Moreover, Y(t) is represented by a non-central chi-squared distribution times a constant C. A European option for the Heston model has formal closed-formula solution by Hull and White[11].

2.2 Three-Factor Model

In the long-term investment, some factors like the national policy, economic development and stock market volatility, will cause the fluctuation of interest rate, then the assumption of the interest rate as a constant will conflict with reality.
In order to overcome the limitations of Black-Scholes model, assuming that the interest rate and volatility are not constant, we consider a mixed model with general stochastic interest rate and stochastic volatility. Under the risk neutral world, the asset price process S(t) is assumed to meet
{dS(t)=r(t)S(t)dt+f(Y(t))S(t)dW(t),dr(t)=μr(r(t))dt+σr(r(t))dZ1(t),dY(t)=μY(Y(t))dt+σY(Y(t))dZ2(t),
(2)
where μr,σr,μY,σY are deterministic functions. W(t),Z1(t),Z2(t) are the standard Brown motions, and their correlation matrix is as follows:
Σ=(dW(t)dW(t)dW(t)dZ1(t)dW(t)dZ2(t)dZ1(t)dW(t)dZ1(t)dZ1(t)dZ1(t)dZ2(t)dZ2(t)dW(t)dZ2(t)dZ1(t)dZ2(t)dZ2(t))=(1ρ1ρ2ρ11ρ12ρ2ρ121)dt.
In this paper, without loss of generality, we assume that ρ12=0. Otherwise, if ρ120, we can make an orthogonal decomposition of dZ2(t):
dZ2(t)=ρ12dZ1(t)+1ρ122dZ~2(t),cov(dZ1(t),dZ~2(t))=0.
We can obtain a new three random process dW(t),\rmdZ1(t),dZ~2(t), their correlation matrix
Σ~=(dW(t)dW(t)dW(t)dZ1(t)dW(t)dZ~2(t)dZ1(t)dW(t)dZ1(t)dZ1(t)dZ1(t)dZ~2(t)dZ~2(t)dW(t)dZ~2(t)dZ1(t)dZ~2(t)dZ~2(t))=(1ρ1ρ2~ρ110ρ2~01)dt.
So the assumption is reasonable. In order to meet the positive definite property of matrix Σ, the condition ρ12+ρ2~2<1 must be satisfied.

3 Conditional Monte Carlo Method

3.1 Conditional Monte Carlo Method Under the Heston Model

Asmussen and Glynn[16] introduced the conditional Monte Carlo (CMC) method as a variance reduction method. Supposing that X and V are two arbitrarily random variables, as the fact that E[V]=E[E[V|X]], and because of the conditional variance formula:
 Var[V]=Var[E[V|X]]+E[Var[V|X]]Var[E[V|X]],
(3)
it's apparently that CMC method can reduce variance of the estimation. One of the difficulties of CMC method is how to choose a suitable condition random variable X, such that the conditional expectation E[V|X] is easy to compute.
For the case of (1), we do orthogonal decomposition of dW(t)
dW(t)=ρdZ(t)+1ρ2dZ~(t),
where cov(dZ(t),dZ~(t))=0. Applying the Ito Lemma to dlnS(t), we can get
 dlnS(t)=(r12Y(t))dt+ρY(t)dZ(t)+1ρ2Y(t)dZ~(t),
(4)
integrating above equality in interval [t,T] implies that
 S(T)=S(t)ξ(t,T)exp[r(Tt)1ρ22tTY(s)ds+1ρ2tTY(s)dZ~(t)],
(5)
where
ξ(t,T)=exp[ρ22tTY(s)ds+ρtTY(s)dZ(s)].
Let the mean variance of the fluctuation rate in the interval [t,T] be
σ¯2(t,T)=1TttTY(s)ds.
Given the stochastic process Z(s), tsT, and Y(t), ξ(t,T), σ¯(t,T) are known, the expectation of the random process Z(s) can be calculated, then we can get the conditional expectation formula of the European call option's price as
V(t,S(t),Y(t))=EZ[VBS(t,S(t)ξ(t,T);r,1ρ2σ¯(t,T))],
(6)
where VBS(t,S;r,σ) is the value of a Vanilla European call option at time t, for stock price S with interest-rate r, market volatility σ, strike price K and expire time T.
VBS(t,S;r,σ)=SN(d1)Ker(Tt)N(d2),
where
d1=lnSK+(r+σ22)(Tt)σTt,d2=d1σTt.
In particular, when t=0, the price of European call option is available as
 V(0,S0,Y0)=EZ[VBS(0,S0ξ(0,T);r,1ρ2σ¯(0,T))].
(7)
By this conditional expectation formula, the value of European option can be calculated, only need to simulate the movement of Brown motion Z(t), so the original two-dimensional problem reduced to a one-dimensional problem, thereby reducing the error and the amount of calculation simultaneously. Similar formula as (7) was obtained by Hull and White[11] for ρ=0.

3.2 An Efficient Control Variate Method Based on the Conditional Expectation by Regression Under the Heston Model

The control variate method is a variance reduction technique used widely in the Monte Carlo simulation. It exploits information about the error of the estimator known quantities to reduce the error of estimator unknown quantity.
Supposing that we want to calculate the expectation of random variable V, X=(X1,,Xd)T is a random vector, and E[X]=(E[X1],E[X2],,E[Xd])T is a known vector. We can generate samples (X(j),V(j))(j=1,,M), which are independent and identical distribution samples from (X,V). Then for fixed bRd, the control variate estimator V(b) is
 V(b)=VbT(XE[X]).
(8)
Glasserman[17] proved that V(b) is an unbiased and consistent estimator of the expectation E[V]. When we choose b=X1XV, the variance of V(b) achieves the minimal value with maximum ratio of variance
 Var[V(b)]=Var[V](1ρXV2)<Var[V],
(9)
where
ρXV2=XVTX1XVσV2,
d×d matrix X is the covariance matrix of X, and XV is the covariance matrix of (X,V).
This suggests that a rather high degree of correlation is needed for a control variate to yield substantial benefits. The efficiency of a control variate depends on the correlation between the option payoff function and the control variable, and the larger correlation coefficient accompanied with the better variance reduction effect. In practice, the optimal vector of coefficient b may be estimated by samples of simulation.
The control variate previously proposed in [18] needs the probability distribution function (pdf) of the volatility process or the variance process, which is not easily obtained by most of the models. Our proposed method is suitable for more general volatility process and can be generalized for the multi-factor model, or even more complex models. Glasserman[17] pointed out that the control variable can be regarded as a regression function. Inspired by this, we propose a control variate method by regression based on conditional expectation formula derived as above.
For simulating the conditional expectation form of the European option
VBS(0,S0ξ(0,T);r,1ρ2σ¯(0,T)),
we introduce a binary function as the control variate by regression idea. Let
X=(ξ(0,T),σ¯2(0,T))T,
where
ξ(0,T)=exp[ρ220TY(s)ds+ρ0TY(s)dZ(s)],
σ¯2(0,T)=1T0TY(s)ds.
The first function is an exponential martingale, and so we have
E[ξ(0,T)]=ξ(0,0)=1.
In the Heston model, by taking the expectation on both side of the second equity of (1), we deduce that
E[Y(t)]=θ+(Y0θ)eat,
thus we have
E[σ¯2(0,T)]=1TEZ[0TY(s)ds]=1T0TE[Y(s)]ds=θ+(Y0θ)1eaTaT.
When the optimal control coefficient b1=b1,b2=b2 are determined, we can use the new estimator with control variate to reduce variance of simulations by:
V(b)=VBS(0,S0ξ(0,T);r,1ρ2σ¯(0,T))b1(ξ(0,T)E[ξ(0,T)])b2(σ¯2(0,T)E[σ2(0,T)]).
(10)
Remark 3.1 When ρ=0,ξ(0,T)=ξ(0,0)=1, so ξ(0,T) can not served as a control variable, we can use σ¯2(0,T) as a single control variable instead to reduce the simulation variance.

3.3 An Efficient Control Variance Method Based on the Conditional Expectation by Regression Under the Three-Factor Model

In this subsection we assume that the interest rate follows the CIR model, and the volatility satisfies the Heston model. Asset price S(t) is governed by the following stochastic differential equations:
{dS(t)=rS(t)dt+Y(t)S(t)dW(t),dr(t)=a1(θ1r(t))dt+σ^1r(t)dZ1(t),dY(t)=a2(θ2Y(t))dt+σ^2Y(t)dZ2(t),
(11)
where a1,θ1,σ^1,a2,θ2,σ^2 are all constants. If the payoff function of a option expired at T is h(S(T)), then the price of the European option at time t, V(t,s,r,y) is defined by
V(t,s,r,y)=E[etTr(s)dsh(S(T))S(t)=s,r(t)=r,Y(t)=y],
(12)
where E represents the expectation under the risk neutral world. Given the initial values of S0,r0 and Y0, we can get the initial price of the European call optiozn with strike price K,
V0=V(0,S0,r0,Y0)=E[e0Tr(s)ds(S(T)K)+].
For convenience, we denote the price of a vanilla European call option under geometry Brown motion model with strike K and expiring time T at time t by VBS(t,S;r,σ). In order to derive a conditional Monte Carlo expressions of (12), we make an orthogonal decomposition of W(t) as
dW(t)=ρ1dZ1(t)+ρ2dZ2(t)+1ρ12ρ22dZ~(t),
where cov(dZ1(t),dZ~(t))=0,\rmcov(dZ2(t),dZ~(t))=0,cov(\rmdZ1(t),dZ2(t))=0, as discussed in the previous section. Applying the Ito formula to dlnS(t),
dlnS(t)=(r(t)12Y(t))dt+ρ1Y(t)dZ1(t)+ρ2Y(t)dZ2(t)+1ρ12ρ22Y(t)dZ~(t).
Then after integrating on both sides of the above equality from t to T, we obtain
S(T)=S(t)ξ1(t,T)ξ2(t,T)exp[tT(r(s)1ρ12ρ222Y(s))ds+1ρ12ρ22tTY(s)dZ~(s)],
(13)
where
ξ1(t,T)=exp[ρ122tTY(s)ds+ρ1tTY(s)dZ1(s)],
ξ2(t,T)=exp[ρ222tTY(s)ds+ρ2tTY(s)dZ2(s)].
Given the information of random process Z(s)=(Z1(s),Z2(s))(tsT) and Y(t), random variables Y(s)(tsT), so ξ1(t,T),ξ2(t,T) are known, then we can get the conditional expectation formula of the European option's price by taking expecting on random process Z~(s) and Z(s) for (S(T)K)+,
V(t,S(t),r(t),Y(t))=EZ[VBS(t,S(t)ξ1(t,T)ξ2(t,T);r¯(t,T),1ρ12ρ22σ¯(t,T))],
(14)
where
r¯(t,T)=1TttTr(s)ds,σ¯2(t,T)=tTY(s)ds
are the average interest rate and the average volatility in interval [t,T] respectively. In particular, when t=0, the European call option's price at the initial time is
V(0,S0,r0,Y0)=EZ[VBS(0,S0ξ1(0,T)ξ2(0,T);r¯(0,T),1ρ12ρ22σ¯(0,T))].
(15)
We now simulate the right hand of above equality using simulation with control variable method by taking the following regression function as the multiple control variates. Let
X=(ξ1(0,T),ξ2(0,T),r¯(0,T),σ¯2(0,T))T,
(16)
where ξ1(0,T),ξ2(0,T) are exponential martingales,
E[ξ1(0,T)]=ξ1(0,0)=1,E[ξ2(0,T)]=ξ2(0,0)=1.
The distribution of r(t) and Y(t) are all proportional to a non center χ2 distribution. Therefore, their expectation are
E[r(t)]=θ1+(r0θ1)ea1t,E[Y(t)]=θ2+(Y0θ2)ea2t.
Thus
E[r¯(0,T)]=θ1+(r0θ1)1ea1ta1T,E[σ¯2(0,T)]=θ2+(Y0θ2)1ea2ta2T.
Using the sample values to determine the optimal coefficients b1,b2,b3,b4, we can take the new estimator of simulation:
V(b)=VBS(0,S0ξ1(0,T)ξ2(0,T);r¯(0,T),1ρ12ρ22σ¯(0,T))b1(ξ1(0,T)E[ξ1(0,T)])b2(ξ2(0,T)E[ξ2(0,T)])b3(r¯(0,T)E[r¯(0,T)])b4(σ¯2(0,T)E[σ¯2(0,T)]).
(17)
It is an unbiased estimation of VBS(0,S0ξ1(0,T)ξ2(0,T),r¯(0,T),1ρ12ρ22σ¯(0,T)). Similarly, when ρ1=0 or ρ2=0, the variables ξ1(t,T) or ξ2(0,T) are constants, which can not be taken as control variables, but we can take σ¯(0,T) and r¯(0,T) as control variates instead to reduce simulation variance.

4 Calculation of Greeks for Option Price

One of the important financial applications of derivatives pricing theory is to be able to manage risk via a liquid options market. Such a market provides the capability for firms and individuals to tailor their risk exposure depending upon their hedging or speculation requirements. In order to effectively assess such risks, it is necessary to calculate the sensitivity of an options price to the factors that affect it, such as the underlying asset price, volatility and time to option expiry, etc. Since all of the sensitivities are commonly denoted by letters of the Greek, called as Greeks.
The first method to calculate the Greeks is to use the finite difference method, but it may meet with the difficulty of non-smooth payoff function[19]. The second method is to use the likelihood ratio method, which is suitable for these density function exists[17]. Fournie, et al.[20] firstly proposed the Greeks calculation formula by using the method of Malliavin segment integral in Wiener space. The main idea of Malliavin integral algorithm is to turn Greeks into a product of profit function and weight function. This method shows great superiority in the calculation of exotic options Greeks[19, 21]. There are also some classical methods to calculate it, such as track simulation, Monte Carlo method, Quasi Monte Carlo method, and so on.
This section will give a simple method for calculation the Greeks of European options based on the conditional expectation formulas.

4.1 Calculation of Greeks for European Call Option Under the Heston Model

In this subsection, we will give a direct simulation method for calculation of Greeks based on the conditional expectation formula under the Heston model. By (6), we can get
V=V(t,s,y)=EZ[VBS(t,sξ(t,T);r,1ρ2σ¯(t,T))]=E[Sξ(t,T)N(d1^)Ker(Tt)N(d2^)],
(18)
where N() is the standard normal distribution function,
d1^=lnSξ(t,T)K+(r+1ρ22σ¯2)(Tt)σ¯(1ρ2)(Tt),d2^=d1^σ¯(1ρ2)(Tt),
ξ(t,T)=exp[ρ22tTY(s)ds+ρtTY(s)dZ(s)],σ¯(t,T)=1TttTY(s)ds.
We will obtain the sensitivity formula by a direct calculation for (18).
Theorem 4.1 The Greeks for a European call option under the Heston model are given by
ΔHE=VS=E[ξ(t,T)N(d1^)],
ΓHE=2VS2=E[ξ(t,T)N(d1^)σ¯S(1ρ2)(Tt)],
ΘHE=Vt=E[rKer(Tt)N(d2^)1ρ2σ¯Sξ(t,T)N(d1^)2Tt].

4.2 Calculation of Greeks for European Call Option Under the Three-Factor Model

Similarly, by the conditional expectation formula of option price (14), we have that
V=V(0,S0,r0,Y0)=EZ[VBS(0,S0ξ1(0,T)ξ2(0,T);r¯(0,T),1ρ12ρ22σ¯(0,T))],
where
Y¯(t)=0TY(s)ds,Y~1(t)=0tY(s)dZ1(s),Y~2(t)=0tY(s)dZ2(s),
ξ1(0,T)=exp[12ρ12Y¯(T)+ρ1Y~1(T)],ξ2(0,T)=exp[12ρ22Y¯(T)+ρ2Y~2(T)],
r¯(0,T)=1T0tr(s)ds,σ¯(0,T)=Y¯(T)T.
We will obtain the following theorem.
Theorem 4.2 Under the three-factor model, the Greeks of European call option can be expressed as the following
ΔTF=VS=E[ξ1(0,T)ξ2(0,T)N(d1^)],ΓTF=2VS2=E[ξ1(0,T)ξ2(0,T)N(d1^)σ¯S0(1ρ12ρ22)T],ΘTF=Vt=E[r¯Ker¯TN(d2^)1ρ12ρ22σ¯S0ξ1(0,T)ξ2(0,T)N(d1^)2T],
(19)
where
d1^=lnS0ξ1(0,T)ξ2(0,T)K+(r¯+1ρ12ρ222σ¯2)Tσ¯(1ρ12ρ22)T,d2^=d1^σ¯(1ρ12ρ22)T.

5 Numerical Results

5.1 Numerical Results Under the Heston Model

5.1.1 Efficient European Call Option Pricing Under the Heston Model

Based on the conditional expectation formula (7), we give the following conditional Monte Carlo simulation algorithm for pricing European call option under the Heston model. In this paper, we use MATLAB for numerical simulations.
Algorithm1(CMC) 1) Divide [0,T] into N parts with mesh size Δt=TN,0=t0<t1<<tN=T.
2) Discreting Y(t),Y¯(t),Y~(t) by the Euler discretizing schemes,
{Y(ti)=Y(ti1)+a(θY+(ti1))Δt+σY+(ti1)ΔtZ(ti),Y¯(ti)=Y¯(ti1)+Y(ti1)Δt,Y~(ti)=Y~(ti1)+Y(ti1)ΔtZ(ti),
for i=1,2,,N, where {Z(ti)} are the standard normal random variables. Y(0)=Y0,Y¯(0)=0,Y~(0)=0. Then
ξ(0,T)=e12ρ2Y¯(tN)+ρY~(tN),σ¯(0,T)=Y¯(tN)/T.
3) Sampling Z(j)(ti) (j=1,2,,M) from Z(ti), we can simulate the European option price by
V(0,S0,Y0)=1Mj=1MVBS[0,S0ξ(j)(0,T);r,1ρ2σ¯(j)(0,T)].
We give Algorithm 2 for pricing European call option by conditional Monte Carlo method with control variate.
Algorithm2(CMCC) 1) Simulating the random variables ξ(0,T),σ¯2(0,T) by {ξ(j)(0,T)}j=1M and {(σ¯(j)(0,T))2}j=1M as in Algorithm 1.
2) Determining the value of b=b^ by simulations: X=(ξ(0,T),σ¯2(0,T))T,V=VBS, we set
b=(b1,b2)=(X)1XV.
3) Simulating option price by
V0=1Mj=1M{VBS(0,S0ξ(j)(0,T);r,1ρ2σ¯(j)(0,T))b1(ξ(j)(0,T)E[ξ(0,T)])b2((σ¯(j))2(0,T)E[σ¯2(0,T)])}.
Experiment1 This experiment will demonstrate the efficiency of Algorithm 1 and Algorithm 2 for pricing a European call option under the Heston model. The parameters of the Heston model are set as α=2,θ=0.01,σ=0.05,Y0=0.015, interest rate r=0.05, initial price of asset S0=30, strike price K=30, maturity T=1 year, simulation path times M=2500, time steps N=50. VMC and VCMCC are the option prices estimated by the classical Monte Carlo method and the conditional Monte Carlo method with control variate respectively. Numerical results are listed in Table 1, where std and t mean standard derivations and total times for simulating the price of options by different methods, respectively. R and stac represent the ratio of standard derivation reduction and comprehensive speedup factor for different values of correlation coefficient ρ in the Heston model (1), defined by
Table 1 Comparisons of simulation error between the MC method and the CMCC method for different values of ρ
ρ VMC stdMC tMC VCMCC stdCMCC tCMCC R stac
0.75 2.1673 0.0226 2.1690 2.1660 0.0049 45.9705 52.92 18.31
0.50 2.1608 0.0236 2.1630 2.1590 0.0023 46.9173 52.25 90.38
0.25 2.1536 0.0237 2.1527 2.1558 0.0006 45.0514 52.52 1795.35
0.01 2.1458 0.0261 2.1647 2.1471 0.0001 45.8104 52.22 460976.35
0.25 2.1383 0.0243 2.1393 2.1373 0.0004 45.6547 52.86 3633.04
0.50 2.1305 0.0277 2.1310 2.1267 0.0018 45.8702 52.38 203.42
0.75 2.1225 0.0271 2.1213 2.1237 0.0043 45.9665 52.22 34.11
R=stdMCstdCMCC,stac=stdMC2tMCstdCMCC2tCMCC=R2tMCtCMCC.
We can see from Table 1 that CMCC can reduce variance obviously when ρ[0.75,0.75], the R is about 2.5700, and the stac can reach about 10106. The more smaller the value of ρ, the higher efficient of the acceleration effect of the Algorithm 1. When ρ∣=1, the comprehensive acceleration effect of CMCC is relatively small, no more than 10 times. The main reason is that we can only remove the random variable 1ρ2Z~ from (5). So the CMC method may no longer work in the case, although the method of control variate still reduces some variance of simulation.
For the purpose of giving a comparison of simulation efficiency between the Algorithm 1 and Algorithm 2, we do the following experiment.
Experiment2 For fixed ρ=0.2 and various values of S0, the other parameters are the same as in Experiment 1. We compare the errors of simulation between Algorithm 1 and Algorithm 2, and the results are listed in Table 2,
R1=stdMCstdCMC,R2=stdMCstdCMCC.
Table 2 Comparisons of simulation error for the MC method, CMC method and the CMCC method for different values of S0
S0 VMC stdMC VCMC stdCMC VCMCC stdCMCC R1 R2
28 0.9917 0.015886 0.9923 0.002251 0.9918 0.000075 7.06 29.97
29 1.5056 0.020967 1.5063 0.002926 1.5061 0.000047 7.17 62.86
30 2.1417 0.024438 2.1436 0.003282 2.1430 0.000030 7.45 107.77
31 2.8911 0.031882 2.8870 0.003576 2.8868 0.000025 8.91 141.65
32 3.7179 0.031834 3.7156 0.003637 3.7158 0.000023 8.75 154.86
Obviously, the ratios of standard deviations between the classical Monte Carlo method and the conditional Monte Carlo method are relatively stable, with a slight increase, however, the ratios of standard deviations between the classical Monte Carlo method and the conditional Monte Carlo method with control variable increase obviously with the increase of the initial value of stock price S0. So the CMCC method peforms better than the CMC method, although it costs a little more extra computing time than the classical Monte Carlo simulations.

5.1.2 Greeks Under the Heston Model

Experiment3 This subsection will demonstrate the difference with three different methods to calculate the value of Δ-hedge and the value of Γ for a European call option under the Heston model. We set parameters M=10000, N=100, θ=0.01, Y0=0.0015, a=2, r=0.05, T=1, ρ=0.75, σ^=0.05, K=30. The results are listed in Table 3 for various values of the stock price, FD means the value of Δ by the method of the forward finite difference with step size h=0.08, CMC and CMCC mean the value of Δ by the method of conditional Monte Carlo and the conditional Monte Carlo method with control variate, Exact means the exact value of Δ.
Table 3 The values of Δ simulated by different methods under the Heston model
S0 28.5 29 29.5 30 30.5 31 31.5
FD 0.5322 0.5975 0.6541 0.7088 0.7045 0.7956 0.8279
CMC 0.5368 0.5995 0.6570 0.7096 0.7569 0.7977 0.8338
CMCC 0.5371 0.5991 0.6569 0.7095 0.7566 0.7978 0.8340
Exact 0.5371 0.5991 0.6569 0.7095 0.7566 0.7980 0.8339
From the simulation results in Tables 35, we can see that, the method of finite difference is the most time-consuming. For simulating the values of Δ and Γ, the methods of CMC and CMCC are performed better than the method of finite difference. Compared to the FD method, it can also be seen that the CMC or CMCC method can not only save computing time but also be effective for reducing the error of simulations.
Table 4 The value of Γ simulated by different methods under the Heston model
S0 28.5 29 29.5 30 30.5 31 31.5
FD 0.1280 0.1266 0.1081 0.0989 0.0800 0.0642 0.0671
CMC 0.1275 0.1200 0.1106 0.0999 0.0885 0.0773 0.0664
CMCC 0.1276 0.1202 0.1107 0.0998 0.0886 0.0772 0.0664
Exact 0.1276 0.1201 0.1106 0.0998 0.0885 0.0772 0.0664
Table 5 Total simulation time of Δ,Γ under the Heston model
method FD CMC CMCC
total time 1.2094s 0.4182s 0.4187s

5.2 Numerical Results Under the Three-Factor Model

5.2.1 Efficient European Call Option Pricing Under Three-Factor Model

This subsection will demonstrate the efficiency for the simulation of option price by the method of CMC and CMCC under the three-factor model. The notations are the same as those in Subsection 4.2.
Algorithm3 (CMC) 1) Dividing [0,T] into N intervals with mesh size Δt=TN,0=t0<t1<<tN=T.
2) Simulating the paths of r(t),Y(t) by the Euler discretizing scheme with the initial values r(0)=r0, Y(0)=y0,
{r(j)(ti)=r(j)(ti1)+a1(θ1r(j)(ti1)+)Δt+σ^1r(j)(ti1)+ΔtZ(j)(ti),Y(j)(ti)=Y(j)(ti1)+a2(θ2Y(j)(ti1)+)Δt+σ^Y(j)(ti1)+ΔtZ(j)(ti),Y¯(j)(ti)=Y¯(j)(ti1)+Y(j)(ti1)Δt,Y~1(j)(ti)=Y~1(j)(ti1)+Y(j)(ti1)ΔtZ1(j)(ti),Y~2(j)(ti)=Y~2(j)(ti1)+Y(j)(ti1)ΔtZ2(j)(ti),r~(j)(ti)=r~(j)(ti1)+r(ti)Δt,
for i=1,2,,N, with Y¯(0)=0,Y~1=Y~2(0)=0,r~(0)=0, and {Z(j)(ti)} are the samples of standard normal random variable. Then
ξ1(j)(0,T)=exp[12ρ12Y¯(j)(tN)+ρ1Y~1(j)(tN)],ξ2(j)(0,T)=exp[12ρ22Y¯(j)(tN)+ρ2Y~2(j)(tN)],
r¯(j)(0,T)=i=1Nr~(j)(ti)T,(σ¯(j))2(0,T)=i=1NY¯(j)(ti)T.
So the j-th sampling value of the European call option under three-factor is
V(j)=VBS(0,S0ξ1(j)(0,T)ξ2(j)(0,T);r¯(j)(0,T),1ρ12ρ22σ¯(j)(0,T)).
3) Summarizing the above simulation values, the option is approximated by
V(0,S0,Y0)=1Mj=1MV(j).
We now give the following conditional Monte Carlo method with the control variable.
Algorithm4 (CMCC) 1) Simulate the values of option price {V(j)}j=1M as we do in Algorithm 3. Meanwhile, we can obtain a set of path samples of processes
ξ1(0,T),ξ2(0,T),r¯(0,T),σ¯2(0,T).
2) Determine the value of b=b^ by simulation: Giving
X=(ξ1(0,T),ξ2(0,T),r¯(0,T),σ¯2(0,T)),V=VBS,
seting
b=(b1,b2,b3,b4)=(X)1XY.
3) Simulate option price by
V=1Mi=1MV(i).
We now take serval experiments to demonstrate the efficiency of the methods of CMC and CMCC.
Experiment4 We set parameters θ1=0.05, θ2=0.01, Y0=0.015, a1=2, a2=2, r0=0.05, K=S0=30, σ^1=0.1, σ^2=0.05, T=1, the remaining parameters are M=10000, N=50.
From Tables 67, we can see that the smaller ρ12+ρ22 is, the better the efficiency of variance reduction of simulation is. When ρ12+ρ22 is closed to 1, CMC method almost has no acceleration effect, and the CMCC method at this time has more than 2 times the standard deviation reduction ratios. This is also an obvious advantage of the CMCC method with control variable over the CMC method. The above simulation results can also be explained by the theoretical deduction in Subsection 3.3, where the conditional Monte Carlo method only eliminates the randomness which is uncorrelated to the stochastic interest rate and the stochastic volatility process. This quantity is proportional to 1ρ12ρ22, so the greater ρ12+ρ22, the smaller elimination of the randomness, and the smaller variance can be expected to reduce.
Table 6 The ratios of standard variance between the MC method and CMC method
ρ1|ρ2 0.7 0.5 0.3 0.1 0.1 0.3 0.5 0.7
0.7 1.08 1.21 1.27 1.24 1.89 1.39 1.35 0.93
0.5 1.30 1.49 2.21 2.16 2.32 2.04 1.45 1.27
0.3 1.21 2.23 2.78 5.14 3.76 2.82 2.08 1.38
0.1 1.75 2.17 4.03 14.11 8.93 3.21 1.99 1.76
0.1 1.46 2.48 3.86 6.06 5.55 3.77 1.86 1.36
0.3 1.31 1.85 2.42 2.49 3.35 2.40 1.52 1.32
0.5 1.20 1.49 1.63 2.17 1.96 1.52 1.57 1.00
0.7 1.08 1.32 1.33 1.76 1.40 1.31 1.09 1.18
Table 7 The ratios of standard variance between the MC method and CMCC method
ρ1|ρ2 0.7 0.5 0.3 0.1 0.1 0.3 0.5 0.7
0.7 2.51 3.85 4.56 7.25 8.58 5.89 4.52 2.93
0.5 3.83 5.10 9.66 13.34 18.92 11.69 6.78 4.30
0.3 4.39 8.49 16.29 44.39 43.57 24.00 11.43 4.74
0.1 5.13 10.60 22.82 188.78 349.31 39.39 13.67 6.55
0.1 4.90 12.09 26.45 65.67 124.82 34.07 11.49 6.92
0.3 4.00 6.45 12.30 23.31 25.97 12.33 7.51 4.59
0.5 3.59 4.29 7.64 10.14 10.83 7.28 5.73 2.93
0.7 2.52 3.72 3.71 6.31 5.78 4.81 3.83 3.37

5.2.2 Greeks Under the Three-Factor Model

Experiment5 For simulating the values of Greeks: Δ,Γ under the three-factor model, we set parameters as M=20000, N=200, θ1=0.05, θ2=0.01, Y0=0.015, a1=2, a2=2, r0=0.05, ρ1=ρ2=0.10, σ^1=0.1, σ^2=0.05, K=30, T=1, h=0.6. The exact values are obtained by 40000 times simulations by the method of CMCC. The simulation results are listed in Tables 810.
Table 8 The values of Δ simulated by different methods under the three-factor model
S0 28.5 29 29.5 30 30.5 31 31.5
FD 0.5195 0.5820 0.6413 0.6965 0.7469 0.7917 0.8300
CMC 0.5197 0.5827 0.6426 0.6980 0.7482 0.7927 0.8316
CMCC 0.5197 0.5827 0.6425 0.6979 0.7482 0.7928 0.8316
Exact 0.5197 0.5827 0.6425 0.6980 0.7482 0.7928 0.8317
Table 9 The values of Γ simulated by different methods under the three-factor model
S0 28.5 29 29.5 30 30.5 31 31.5
FD 0.1284 0.1237 0.1158 0.1074 0.0949 0.0819 0.0698
CMC 0.1283 0.1233 0.1156 0.1059 0.0949 0.0834 0.0720
CMCC 0.1283 0.1233 0.1156 0.1059 0.0949 0.0834 0.0720
Exact 0.1283 0.1233 0.1156 0.1059 0.0950 0.0834 0.0720
Table 10 Total simulation time of Δ,Γ under the three-factor model
method FD CMC CMCC
total time 2.1454s 0.6673s 0.6944s
We can see from the Tables 89 that the methods of CMC and CMCC exhibit smaller simulation errors than the method of FDM. Moreover, the calculation time of finite difference method is much larger than the other two methods. So we can see that the method of CMCC or CMC peforms much better than that of FDM.
Remark 5.1 From above numerical results as in Tables 810, it is clear that CMCC method has almost the same variance reduction effect as the method of CMC, and requires more time for the calculation of Greeks.
Overall, the methods of CMC and CMCC can reduce the variance, that we have demonstrated both theoretically and empirically. First, we used formulas (3) and (9) to explain that the variance obtained by CMC and CMCC are significantly smaller than before. Second, in numerical simulation, based on the empirical research results on China's securities market, the parameters selected are consistent with the parameter ranges of existing research results[22-24], we also proved the effectiveness of these methods in reducing variance.

6 Conclusion

In this paper, we have introduced the method of conditional Monte Carlo and the method of conditional Monte Carlo with control variate, for acceleration of pricing European options under the Heston model and the three-factor model. The CMCC method combining the control variate method with conditional Monte Carlo method, which can be easily obtained by the idea of regression. Besides of option pricing, CMC method can also be applied to sensitivity analysis. Numerical results indicate that the CMCC method and CMCC method have efficiency variance reduction effect for options pricing, besides the CMC method is also effective for the calculation of Greeks under the Heston model and the three-factor model.
The idea of conditional Monte Carlo method and method of conditional Monte Carlo with control variate introduced in the paper may also be applicable for pricing options written on several assets' price, such as the Basket option, the Asian option with discrete observing points and the Barrier option with discrete observing points and other path dependent options under the Heston model and the three-factor model respectively. For example, for pricing a call Basket option with several assets, the payoff function is (AK)+, where A=1ni=1nSiT is the arithmetic average price of asset prices SiT(i=1,2,,n) at the terminal time T. Then the price of Basket option can be expressed as
V=erTE[(AK)+]=erT(E[(AK)1{G>K}]+E[(AK)1{G<KA}])ξ1+ξ2,
as suggested by Curran[25], where G=Πi=1nSiTn is the geometric average price of SiT. Since the expression ξ1 can be rewritten as a conditional expectation formula, then the conditional Monte Carlo simulation can be used for the purpose of variance reduction. ξ2 is the expectation of a small probability problem, which can be simulated by importance technique naturally, we will report the result in the future.

References

1
Chen N, Huang Z. Localization and exact simulation of Brownian motion-driven stochastic differential equations. Mathematics of Operations Research, 2013, 38(3): 591- 616.
2
Duffie D, Glynn P W. Efficient Monte Carlo simulation of security prices. The Annals of Applied Probability, 1995, 5(4): 897- 905.
3
Kemna A G Z, Vorst A C F. A pricing method for options based on average asset values. Journal of Banking & Finance, 1990, 14(1): 113- 129.
4
Longstaff F A, Schwartz E S. Valuing American options by simulation: A simple least-squares approach. The Review of Financial Studies, 2001, 14(1): 113- 147.
5
Fouque J P, Tullie T A. Variance reduction for Monte Carlo simulation in a stochastic volatility environment. Quantitative Finance, 2002, 2(1): 24- 30.
6
Bolia N, Juneja S, Glasserman P. Function-approximation-based importance sampling for pricing American options. Winter Simulation Conference, 2004, 1(4): 604- 611.
7
Broadie M, Kaya O. Exact simulation of stochastic volatility and other affine jump diffusion processes. Operations Research, 2006, 54(2): 217- 231.
8
Heston S L. A closed-form solution for options with stochastic volatility with applications to bond and currency options. The Review of Financial Studies, 1993, 6(2): 327- 343.
9
Amin K I, Jarrow R A. Pricing options on risky assets in a stochastic interest rate economy. Mathematical Finance, 1992, 2(4): 217- 237.
10
Hull J. Options futures and other derivatives. 6th ed. Prentice Hall, New Jersey, 2006.
11
Hull J, White A. The pricing of options on assets with stochastic volatilities. The Journal of Finance, 1987, 42(2): 281- 300.
12
Scott L O. Option pricing when the variance changes randomly: Theory, estimation, and an application. The Journal of Financial and Quantitative Analysis, 1987, 22(4): 419- 438.
13
Stein E M, Stein J C. Stock price distributions with stochastic volatility: An analytic approach. The Review of Financial Studies, 1991, 4(4): 727- 752.
14
Ball C, Roma A. Stochastic volatility option pricing. The Journal of Financial and Quantitative Analysis, 1994, 29(4): 589- 607.
15
Cox J C, Ingesoll J E, Ross S A. A theory of the term structure of interest rates. Econometrica: Journal of the Econometric Society, 1985, 53(2): 385- 407.
16
Asmussen S, Glynn P W. Stochastic simulation: Algorithms and analysis. New York: Springer, 2007.
17
Glasserman P. Monte Carlo methods in financial engineering. Springer, 2003.
18
Fouque J, Han C. A martingale control variate method for option pricing with stochastic volatility. ESAIM Probability and Statistics, 2007, 11, 40- 54.
19
Xu Y J, Lai Y Z, Yao H. Efficient simulation of Greeks of multiasset European and Asian style options by Malliavin calculus and quasi-Monte Carlo methods. Applied Mathematics and Computation, 2014, 236(1): 493- 511.
20
Fournie E, Lasry J M, Lebuchoux J, et al. Applications of Malliavin calculus to Monte Carlo methods in finance Ⅰ. Finance Stoch, 1999, 3(4): 391- 412.
21
Fournie E, Lasry J M, Lebuchoux J, et al. Applications of Malliavin calculus to Monte Carlo methods in finance Ⅱ. Finance Stoch, 2001, 5(2): 201- 236.
22
Chao W, Qian X T. Simulation of catastrophe insurance bankruptcy probability based on Monte Carlo algorithm under random interest rates. Journal of Central University for Nationalities (Natural Science Edition), 2021, 30(2): 47- 51.
23
Jiang L, Lin H X. Parameter estimation of short-term interest rate models based on P-spline method. Journal of Engineering Mathematics, 2015, 32(4): 485- 496.
24
Zhang X J, Chen H Z, Jiang L. Parameter estimation of stochastic volatility model based on generalized moment method. Journal of Engineering Mathematics, 2019, 36(6): 611- 626.
25
Curran M. Valuing Asian and portfolio options by conditioning on the geometric mean price. Management Science, 1994, 40(12): 1705- 1711.
PDF(184 KB)

572

Accesses

0

Citation

Detail

Sections
Recommended

/