URL: http://www.emath.fr/ps/
A MARTINGALE CONTROL VARIATE METHOD FOR OPTION PRICING
WITH STOCHASTIC VOLATILITY
∗, ∗∗Jean-Pierre Fouque
1and Chuan-Hsiang Han
2Abstract. A generic control variate method is proposed to price options under stochastic volatility models by Monte Carlo simulations. This method provides a constructive way to select control variates which are martingales in order to reduce the variance of unbiased option price estimators. We apply a singular and regular perturbation analysis to characterize the variance reduced by martingale control variates. This variance analysis is done in the regime where time scales of associated driving volatility processes are well separated. Numerical results for European, Barrier, and American options are presented to illustrate the effectiveness and robustness of this martingale control variate method in regimes where these time scales are not so well separated.
1991 Mathematics Subject Classification. 65C05, 62P05. September 26, 2005. Revised May 2, 2006.
Dedicated to Nicole Elkaroui in honor of her 60th birthday
Introduction
Monte Carlo pricing for options is a popular approach in particular since efficient algorithms have been developed for optimal stopping problems, see for example [10] . The advantage of Monte Carlo simulations is that it is less sensitive to dimensionality of the pricing problems and suitable for parallel computation; the main disadvantage is that the rate of convergence is limited by the central limit theorem so it is slow.
To increase the efficiency besides parallel computing, Quasi Monte Carlo and variance reduction techniques are two possible approaches. We refer to [9] for an extensive review. Quasi Monte Carlo, unlike pseudo-random number generators, forms a class of methods where low-discrepancy numbers are generated in deterministic ways. Its efficiency heavily relates to the regularity of the option payoffs, which in most cases are poorly posted. The pros are that such an approach can be always implemented regardless to the pricing problems and it is easy to combine with other sampling techniques such as those involving the Brownian bridge. On the other hand, variance reduction methods seek probabilistic ways to reformulate the pricing problem considered in order to gain significant variance reduction. For example control variate methods take into account correlation properties of random variables, and importance sampling methods utilize changes of probability measures. The cons are that the efficiency of these techniques is often restricted to certain pricing problems.
Keywords and phrases: ...
∗The work of the first author is supported by NSF grant DMS-0455982
∗∗ The work of the second author is supported by NSC grant 94-2119-M-007-008, Taiwan
1Department of Statistics and Applied Probability, University of California, Santa Barbara, CA 93106-3110, [email protected] 2 Department of Quantitative Finance, National Tsing Hua University, Hsinchu, Taiwan, 30013, ROC, [email protected]
c
Stochastic volatility models have been an important class of diffusions extending the Black-Scholes model, see [6] for details. Under multifactor stochastic volatility models, this paper aims at generalizing the control variate method proposed by the authors in [4], and studying its variance analysis. Since the proposed control variates are (local) martingales, we shall call this method “martingale control variate method”. The pricing problems of European, Barrier and American options are considered in order to demonstrate the effectiveness of our method for a broad range of problems.
The martingale control variate method can be well understood in finance terminology. The constructed control variate corresponds to a continuous (non-self-financing) delta hedge strategy taken by a trader who sells an option. Though perfect replication by delta hedging under stochastic volatility models is impossible, the variance of replication error is directly related to the variance induced by the martingale control variate method. This method is also potentially useful to study contracts dealing with volatility or variance risks such as variance swaps.
A variance analysis, presented in the Appendix, deduces an asymptotic result for the variance reduced by martingale control variates. It is based on the singular and regular perturbation method presented in [8]. The paper is organized as follows. In Section 1 we introduce the basic Monte Carlo pricing mechanism and review the martingale control variate method for European options. Section 2 and 3 extends the method to Barrier and American options, respectively. Numerical experiments are included and we conclude this paper in Section 4.
1.
Monte Carlo Pricing under Multiscale Stochastic Volatility Models
Under a risk-neutral pricing probability measure IP⋆parametrized by the combined market prices of volatility risk (Λ1, Λ2) , we consider the following class of multiscale stochastic volatility models:
dSt = rStdt + σtStdWt(0)∗, (1) σt = f (Yt, Zt), dYt = 1 εc1(Yt) + g1(Yt) √ε Λ1(Yt, Zt) dt + g1√(Yt) ε ρ1dWt(0)∗+ q 1− ρ2 1dWt(1)∗ , dZt = h δc2(Zt) + √ δg2(Zt)Λ2(Yt, Zt) i dt + √δg2(Zt) ρ2dWt(0)∗+ ρ12dWt(1)∗+ q 1− ρ2 2− ρ212dWt(2)∗ ,
where Stis the underlying asset price process with a constant risk-free interest rate r. Its stochastic volatility σt
is driven by two stochastic processes Ytand Ztvarying on the time scales ε and 1/δ, respectively (ε is intended
to be a short time scale while 1/δ is thought as a longer time scale). The vectorWt(0)∗, Wt(1)∗, Wt(2)∗
consists of three independent standard Brownian motions. The instant correlation coefficients ρ1, ρ2, and ρ12 satisfy
|ρ1| < 1 and |ρ22+ ρ212| < 1. The volatility function f is assumed to be smooth bounded and bounded away from
0. The coefficient functions of Yt, namely c1and g1, are assumed to be such that under the physical probability
measure, Yt is ergodic. The Ornstein-Uhlenbeck (OU) process is a typical example by defining c1(y) = m1− y
and g1(y) = ν1√2 such that 1/ε is the rate of mean reversion, m1 is the long run mean, and ν1is the long run
standard deviation. Its invariant distribution isN m1, ν12 .
The coefficient functions of Zt, namely c2and g2are assumed to be smooth enough in order to satisfy existence
and uniqueness conditions for diffusions [11]. The combined risk premia Λ1 and Λ2 are assumed to be smooth,
bounded, and depending on the variables y and z only. Within this setup, the joint process (St, Yt, Zt) is
Markovian. We refer to [8] for a detailed discussion on this class of models.
payoff function H and expiry T is given by Pε,δ(t, x, y, z) = IE⋆t,x,y,z n e−r(T −t)H(ST) o , (2)
where IE⋆t,x,y,z denotes the expectation with respect to IP⋆ conditioned on the current states St = x, Yt =
y, Zt= z. A basic Monte Carlo simulation estimates the option price Pε,δ(0, S0, Y0, Z0) at time 0 by
1 N N X i=1 e−rTH(ST(i)), (3)
where N is the total number of independent sample paths and ST(i) denotes the i-th simulated stock price at time T .
Assuming that the European option price Pε,δ(t, S
t, Yt, Zt) is smooth enough, we apply Ito’s lemma to its
discounted price e−rtPε,δ, and then integrate from time 0 to the maturity T . The following martingale repre-sentation is obtained Pε,δ(0, S0, Y0, Z0) = e−rTH(ST)− M0(Pε,δ)− 1 √εM1(Pε,δ)− √ δM2(Pε,δ), (4)
where centered martingales are defined by M0(Pε,δ) = Z T 0 e−rs∂P ε,δ ∂x (s, Ss, Ys, Zs)f (Ys, Zs)SsdW (0)∗ s , (5) M1(Pε,δ) = Z T 0 e−rs∂P ε,δ ∂y (s, Ss, Ys, Zs)g1(Ys)d ˜W (1)∗ s , (6) M2(Pε,δ) = Z T 0 e−rs∂P ε,δ ∂z (s, Ss, Ys, Zs)g2(Zs)d ˜W (2)∗ s , (7)
with the Brownian motions ˜ Ws(1)∗ = ρ1Ws(0)∗+ q 1− ρ2 1Ws(1)∗, ˜ Ws(2)∗ = ρ2Ws(0)∗+ ρ12Ws(1)∗+ q 1− ρ2 1− ρ212Ws(2)∗.
These martingales play the role of “perfect” control variates for Monte Carlo simulations and their integrands would be the perfect Delta hedges if Pε,δwere known and volatility factors traded. Unfortunately neither the
option price Pε,δ(s, S
s, Ys, Zs) nor its gradient at any time 0≤ s ≤ T are in any analytic form even though all
the parameters of the model have been calibrated as we suppose here.
One can choose an approximate option price to substitute Pε,δ used in the martingales (5, 6, 7) and still
retain martingale properties. When time scales ε and 1/δ are well separated, namely 0 < ε ≪ 1 ≪ 1/δ, an approximation of the Black-Scholes type is derived in [8]:
Pε,δ(t, x, y, z)≈ PBS(t, x; ¯σ(z)) (8)
with an accuracy of order O√ε,√δ for continuous payoffs. We denote by PBS(t, x; ¯σ(z)) the solution of
the Black-Scholes partial differential equation with the terminal condition PBS(T, x) = H(x). The z-dependent
effective volatility ¯σ(z) is defined by ¯
where Φ(y) is the invariant distribution of the fast varying process Yt. In the OU case, the density Φ is simply
the Gaussian density with mean m1 and variance ν12. Note that the approximate option price PBS(t, x; ¯σ(z)) is
independent of the variable y. A martingale control variate estimator is formulated as 1 N N X i=1 h e−rTH(ST(i))− M(i)0 (PBS)− √ δM(i)2 (PBS) i . (10)
This is the approach taken by Fouque and Han [4], in which the proposed martingale control variate method is empirically superior to an importance sampling [5] for pricing European options. As control variates M0
and M2are martingales, we shall call them martingale control variates afterwards. Note that there is noM1
martingale term since the approximation PBS does not depend on y and the y-derivative cancels in (6) with
Pε,δ replaced by P BS.
1.1. Variance Analysis of Martingale Control Variates
Since M2(PBS) is small of order√δ, in a first approximation we can neglectM2(PBS) in (10). Hence we
reduce the number of stochastic integrals or martingale control variates from 2 to 1 and formulate the following unbiased estimator: 1 N N X i=1 h e−rTH(S(i) T )− M (i) 0 (PBS) i , (11) where M0(PBS) = Z T 0 e−rs∂PBS ∂x (s, Ss; ¯σ(Zs))f (Ys, Zs)SsdW (0)∗ s .
For the sake of simplicity, we first assume that the instant correlation coefficients, ρ1, ρ2 and ρ12 in (1), are
zero. From (4), the variance of the controlled payoff
e−rTH(ST)− M0(PBS) (12)
is simply the sum of quadratic variations of martingales: V ar e−rTH(ST)− M0(PBS) (13) = IE⋆0,t,x,y,z ( Z T 0 e−2rs ∂Pε,δ ∂x − ∂PBS ∂x 2 (s, Ss, Ys, Zs)f2(Ys, Zs)Ss2ds + 1 ε Z T 0 e−2rs ∂P ε,δ ∂y 2 (s, Ss, Ys, Zs)g21(Ys)ds + δ Z T 0 e−2rs ∂P ε,δ ∂z 2 (s, Ss, Ys, Zs)g22(Zs)ds ) . (14)
As in the numerical experiments implemented in [4] and in next Sections, we assume that the driving volatility processes Ytand Ztare of OU type; namely c1(y) = (m1−y), c1(z) = (m2−z), g1(y) = ν1√2, and g2(z) = ν2√2.
The volatility premia Λ1 and Λ2 are assumed to be smooth and bounded.
Theorem 1.1. Under the assumptions made above and the payoff function H being continuous piecewise smooth as a call (or a put), for any fixed initial state (0, x, y, z), there exists a constant C > 0 such that for ε≤ 1, δ ≤ 1,
V ar e−rTH(S
The proof of Theorem 1.1 is given in the Appendix. We comment this theorem:
(1) The assumption of zero instant correlations is not necessary. One can still obtain the same accuracy result with additional cross-variation terms appearing in equation (13).
(2) Adding the next order corrections in√ε and√δ to (8), as suggested in [8], and using two martingale control variates as in (10), we would obtain that the variance associated with the estimator is still of the same order as in the Theorem. One can obtain next order accurate result for Lemma A.1. However there is no accuracy gain for Lemma A.2 because the next order price approximation is still independent of the fast varying y-variable [8].
Several variance reduction results for pricing European call options can be found in [4], where the martingale control variate method does demonstrate significant variance reduction performance when time scales are well separated.
From the computational viewpoint, since calculating each stochastic integral along a sample path is time con-suming, it is useful to reduce the number of stochastic integrals from (10) to (11) and still retain considerable accuracy for the reduced variance. From the finance point of view, the martingale control variate M0(PBS)
represents that a trader, who sells an option, uses the delta hedge strategy continuously. By doing so, the induced error of replicated discounted-payoff e−rTH(ST)− M0(PBS) and its statistical property can be studied
through the Monte Carlo simulations (11). Since the martingale control variate method is associated with hedging strategies, it should, in principle, work for all other derivatives pricing problems provided the delta is easy enough to be computed or effectively approximated.
In the next sections, we generalize this method to Barrier and American option pricing problems under stochastic volatility models.
2.
Barrier Options
The payoff of a barrier option depends on whether the trajectory of the underlying stock hits a pre-specified level or not before the maturity T . For instance a down and out call option with the barrier B and the strike K has a payoff
ST − K+ I{τ >T },
where we denote by I the indicator function and by τ the first hitting time τ = inf{0 ≤ t ≤ T, St≤ B}.
Other popular barrier options such as down and in, up and out and up and in can be defined similarly. Under the risk-neutral probability IP⋆, a down and out barrier call option price at time t conditioning on no knock-out before time t < T is given by
Pε,δ(t, x, y, z) = IE⋆t,x,y,z
n
e−r(T −t) ST − K+ I{τ >T }
o
. (15)
The price Pε,δ(t, x, y, z) solves a boundary value problem [7]. When parameters ε and δ are small enough, the
leading order approximation to Pε,δ in (15) is given by
Pε,δ(t, x, y, z)≈ PB
BS(t, x; ¯σ(z)), (16)
where PB
BS(t, x; ¯σ(z)) solves a Black-Scholes partial differential equation for a barrier option problem with the
effective volatility ¯σ(z), and the boundary conditions PB
for x > B. It is known (see for instance [12]) that PB
BS(t, x; ¯σ(z)) admits the closed form solution
PBSB (t, x; ¯σ(z)) = CBS(t, x; ¯σ(z))− x B 1−k CBS(t, B2/x; ¯σ(z)), (17) where k = 2r/(¯σ2(z)) and C
BS(t, x; ¯σ(z)) denotes the Black-SCholes price of a European call option with strike
K, maturity T , and volatility ¯σ(z).
2.1. Martingale Control Variate Estimator for Barrier Options
Let S0> B, one can apply Ito’s lemma to the discounted barrier option price, then integrate from time 0 up
to the bounded stopping time τ∧ T so that
Pε,δ(0, S0, Y0, Z0) = e−rT(ST− B)+I{τ >T } (18) −M0(Pε,δ)− 1 √εM1(Pε,δ)− √ δM2(Pε,δ)
is deduced. The local martingales are defined as in (5, 6, 7) except that the upper bounds are replaced by τ∧T . As in Section 1.1, we use the barrier price approximation (17) to construct the following local martingale control variate M0(PBSB ) = Z τ ∧T 0 e−rs∂PBSB ∂x (s, Ss; ¯σ(Zs))f (Ys, Zs)SsdW (0)∗ s .
The unbiased martingale control variate estimator by Monte Carlo simulations for the barrier option is 1 N N X i=1 h e−rT(ST(i)− K)+I {τ(i)>T }− M (i) 0 (PBSB ) i .
The variance analysis for the estimator e−rT(S
T− K)+I{τ >T }− M0(PBSB ) can be done similarly as in Theorem
1.1. In fact one can obtain the same accuracy, namely O(ε, δ), because, as shown in [7], the accuracy of the leading order barrier option approximation in (16) is the same as for European options. All other derivations remain the same.
2.2. Numerical Results
Several numerical experiments are presented to demonstrate that the martingale control variate method is efficient and robust for barrier option problems even in the regimes where the time scales ε and 1/δ are not so well separated. Relevant parameters and volatility functions for a two-factor stochastic volatility model are chosen as in Table 1. Other values including initial conditions and option parameters are given in Table 2. Option price computations are done with various time scale parameters given in Table 3. The sample size is N = 10, 000. Simulated paths are generated based on the Euler discretization scheme [9] with time step size ∆t = 10−3. Figure 1 presents sampled barrier option prices with respect to the number of realizations. The
dash line corresponds to basic Monte Carlo simulations, while the solid line corresponds to same Monte Carlo simulations using the martingale control variateM0(PBSB ).
3.
American Options
The most important feature of an American option is that the option holder has the right to exercise the contract early. Under the stochastic volatility models considered, the price of an American option with the
Table 1. Parameters used in the two-factor stochastic volatility model (1).
r m1 m2 ν1 ν2 ρ1 ρ2 ρ12 Λ1 Λ2 f (y, z)
8% -1 -1 1 1 -0.2 -0.2 0 0 0 exp(y + z) Table 2. Initial conditions and down and out barrier call option parameters.
$S0 Y0 Z0 $K $B T years
100 -1 -1 110 80 1
Table 3. Comparison of standard errors with various ε and δ. The notation StdBMC stands
for the standard error estimated from basic Monte Carlo simulations, and StdMCV the standard
error from the same Monte Carlo simulations but using the martingale control variate. Numbers within parenthesis in the third and fourth columns are sample means estimated from the two Monte Carlo methods, respectively. The fifth column records the variance reduction ratio, which is calculated by StdBMC/StdMCV2
.
1/ε δ StdBMC StdMCV Variance Reduction Ratio
100 0.01 0.2822 (10.8153) 0.0304 (10.8497) 86 75 0.1 0.2047 (10.7652) 0.0306 (10.7594) 45 50 1 0.2455 (11.2082) 0.0474 (11.0962) 27 25 10 0.2604 (12.6212) 0.0417 (12.4372) 39 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 9 9.5 10 10.5 11 11.5 12 Number of Realizations
Down−and−Out Barrier Option Prices
Basic Monte Carlo Martingale Control Varite
Figure 1. Monte Carlo simulations for a down-and-out barrier call option price when 1/ε = 75 and δ = 0.1. Sampled prices are obtained along the number of realizations.
payoff function H is given by: Pε,δ(t, x, y, z) = (ess) sup t≤τ ≤T IE⋆t,x,y,z n e−r(τ −t)H(Sτ) o , (19)
where τ denotes any stopping time greater than t, bounded by T , adapted to the completion of the natural filtration generated by Brownian motionsWt(0)∗, Wt(1)∗, Wt(2)∗
. We consider a typical American put option pricing problem, namely H(x) = (K− x)+, and maturity T . By the connection of optimal stopping problems
and variational inequalities [11], Pε,δ(t, x, y, z) can be characterized as the solution of the following variational
inequalities L(S,Y,Z)Pε,δ(t, x, y, z)≤ 0 Pε,δ(t, x, y, z) ≥ (K − x)+ L(S,Y,Z)Pε,δ(t, x, y, z)· Pε,δ(t, x, y, z)− (K − x)+ = 0,
whereL(S,Y,Z)denotes the infinitesimal generator of the Markov process (St, Yt, Zt) . The optimal stopping time
is characterized by
τ⋆(t) = inft ≤ s ≤ T, (K − S
s)+= Pε,δ(s, Ss, Ys, Zs) . (20)
When ε and δ are small enough, the leading order approximation by a formal expansion is Pε,δ(t, x, y, z)≈ PA
BS(t, x; ¯σ(z)) (21)
while PA
BS(t, x; ¯σ(z)) solves the homogenized variational inequality
LBS(¯σ(z))PBSA (t, x; ¯σ(z))≤ 0 PA BS(t, x; ¯σ(z))≥ (K − x)+ LBS(¯σ(z))PBSA (t, x; ¯σ(z))· PBSA (t, x; ¯σ(z))− (K − x)+ = 0, (22)
where LBS(¯σ(z)) denotes the Black-Scholes operator with the constant volatility ¯σ(z). In contrast to typical
European and barrier options, there is no closed-form solution for the American put option price under a constant volatility. The derivation of the accuracy of the approximation (21) is still an open problem.
As in the previous sections, we assume that the discounted American option price e−rtPε,δ(t, S
t, Yt, Zt) before
exercise is smooth enough to apply Ito’s lemma, then we integrate from time 0 to the (bounded) optimal stopping time τ⋆ such that we obtain
Pε,δ(0, S0, Y0, Z0)
= e−rT(K− Sτ⋆)+− M0(Pε,δ)−√1
εM1(P
ε,δ)
−√δM2(Pε,δ). (23)
The local martingales are defined as in (5, 6, 7) except that the upper bounds are replaced by the optimal stopping time τ⋆.
3.1. Martingale Control Variates for American Options
As revealed in (20), the characterization of the optimal stopping time τ⋆(t) does depend on the American
option price, which itself is unknown in advance. This causes an immediate difficulty to implement Monte Carlo simulations because one does not know the time to stop in order to collect the payoff along each realized sample path.
Longstaff and Schwartz [10] took a dynamic programming approach and proposed a least-square regression to estimate the continuation value at each in-the-money stock price state.
By comparing the continuation value and the instant exercise payoff, their method exploits a decision rule, denoted by τ , for early exercise along each sample path generated. It is shown in [10] that Longstaff-Schwartz’ method provides a low-biased American option price estimate for practical Monte Carlo simulations. As the number of least-square basis functions increases to infinity for discrete exercise dates, Clement et al. in [3] show that the normalized error in Longstaff-Schwartz’ method is asymptotically Gaussian.
Like in previous sections, a local martingale control variate can be in principle constructed as
M0(PBSA ; τ⋆) = Z τ⋆ 0 e−rs∂P A BS ∂x (s, Ss; ¯σ(Zs))f (Ys, Zs)SsdW (0)∗ s .
Indeed the optimal stopping time τ⋆ is not known. We approximate τ⋆ by the exercise rule τ . Note that
M0(PBSA ; τ ) may incur a bias but from the sample means in Table 6 this effect seems negligible. In fact, one
could build an approximate stopping time from Longstaff-Schwartz’ method. This can be done but will be computational expensive.
There is no closed-form solution for the homogenized American option PA
BS(t, x; ¯σ(z)) either. We introduce
an approximation proposed by Barone-Adesi and Whaley [1], denoted by PBAW
BS , which is derived from an
elliptic-type variational inequalities as an approximation to the parabolic-type variational inequalities (22). The approximation PBAW
BS admits the closed-form solution:
PBAW BS (t, x; σ) = λxα+ PE BS(t, x; σ), x > x∗ K− x, x≤ x∗, where PE
BS(t, x; σ) denotes the corresponding European put option price and the free boundary x∗ solves the
following nonlinear algebraic equation
x∗=|α| K− P E BS(t, x∗; σ) ∂PE BS(t,x∗;σ) ∂x + 1 +|α| , with α = 1− 2r σ2 − q (1−2r σ2)2+ 8(κr+1) κσ2 2 and λ = K− x∗− P E BS(t, x; σ) (x∗)α .
To summarize, we construct the following stopped martingale as a control variate
M0(PBSBAW; τ ) = Z τ 0 e−rs∂P BAW BS ∂x (s, Ss; ¯σ(Zs))f (Ys, Zs)SsdW (0)∗ s .
The Monte Carlo estimator with the martingale control variate is 1 N N X i=1 h e−rτ(K− Sτ(i))+− M (i) 0 (P BAW BS ; τ ) i .
Table 4. Parameters used in the two-factor stochastic volatility model (1).
r m1 m2 ν1 ν2 ρ1 ρ2 ρ12 Λ1 Λ2 f (y, z)
10% -1 -1 1 1 -0.3 -0.3 0 0 0 exp(y + z) Table 5. Initial conditions and American put option parameters.
$S0 Y0 Z0 $K T years
90 -1 -1 100 1
Table 6. Comparison of standard errors with various ε and δ. The notation StdBMC stands
for the standard error estimated from basic Monte Carlo simulations and StdMCV the standard
error from same Monte Carlo simulations but using the martingale control variate. Numbers within the parenthesis in the third and fourth columns are sample means estimated from the two Monte Carlo methods, respectively. The fifth column records the variance reduction ratio, which is calculated by StdBMC/StdMCV2
.
1/ε δ StdBMC StdMCV Variance Reduction Ratio
100 0.01 0.2354 (21.4340) 0.0240 (21.5942) 96 75 0.1 0.2564 (21.4791) 0.0286 (21.8001) 81 50 1 0.2571 (21.5217) 0.0350 (21.6319) 54 25 10 0.2606 (21.9621) 0.0453 (21.3243) 32
3.2. Numerical Results
We consider American put options under two factor stochastic volatility models, specified in Table 4 and Table 5. Results of variance reduction by the martingale control variate to price American put options are illustrated in Table 6 with various time scale parameters ε and δ. The discrete time step size is ∆ = 10−3
and the total sample size is N = 5, 000. Figure 2 presents sampled American put option prices with respect to the number of realizations. The dash line corresponds to basic Monte Carlo simulations, while the solid line corresponds to same Monte Carlo simulations using martingale control variateM0(PBSBAW).
4.
Conclusion
In the context of multifactor stochastic volatility models, we propose a martingale control variate method to price options by Monte Carlo simulations. A theoretical variance analysis is provided to characterize the size of the variance reduced by martingale control variate in the case that driving volatility time scales are well separated. Comparing to plain Monte Carlo simulations, significant variance reduction ratios for European, Barrier and American options are obtained even in regimes where volatility time scales are not so well separated. These results confirm the robustness of our method based on martingale control variates constructed as in delta hedging strategies. The effectiveness of our method depends on option price approximations to the pricing problem considered.
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 19.5 20 20.5 21 21.5 22 22.5 23 23.5 Number of Realizations
American Put Option Prices
Basic Monte Carlo Martingale Control Varite
Figure 2. Monte Carlo simulations for an American put option price when 1/ε = 50 and δ = 1. Sampled prices are obtained along the number of realizations.
Appendix A. Derivation of the accuracy of the variance analysis
In order to prove Theorem 1.1, we need the following three lemmas.Lemma A.1. Under the assumptions of Theorem 1.1, for any fixed initial state (0, x, y, z), there exists a positive constant C1> 0 such that for ε≤ 1 and δ ≤ 1, one has
IE⋆0,t,x,y,z ( Z T 0 e−2rs ∂P ε,δ ∂x − ∂PBS ∂x 2 (s, Ss, Ys, Zs)f2(Ys, Zs)Ss2ds ) ≤ C1max{ε, δ}
Proof: By Cauchy-Schwartz inequality we have
IE⋆0,t,x,y,z ( Z T 0 e−2rs ∂P ε,δ ∂x − ∂PBS ∂x 2 (s, Ss, Ys, Zs)f2(Ys, Zs)Ss2ds ) (24) ≤ s IE⋆ Z T 0 ∂Pε,δ ∂x − ∂PBS ∂x 4 (s, Ss, Ys, Zs)ds × s Z T 0 IE⋆nf4(Y s, Zs) (e−rsSs)4 o ds,
where we omitted the sub-scripts under the expectation IE⋆. The second factor on the right hand side of this inequality is bounded by s Z T 0 IE⋆nf4(Y s, Zs) (e−rsSs)4 o ds≤ Cf s Z T 0 IE⋆n(e−rsSs)4ods (25)
for some constant Cf, as the volatility function f is bounded. Using the notation σt= f (Yt, Zt) as in (1), and
W(0)∗ = W for simplicity, one has
e−rsS s= xe Rs 0σudWu−12 Rs 0σ 2 udu, and therefore IE⋆n e−rsSs 4o = x4IE⋆ne6R0sσ 2 udue Rs 04σudWu−12 Rs 016σ 2 udu o ≤ C′ fx4IE⋆ n eR0s4σudWu−12 Rs 016σ 2 udu o = Cf′x4,
where we have used again the boundedness of f , and the martingale property. Combined with (25) we obtain s Z T 0 IE⋆nf4(Y s, Zs) (e−rsSs)4 o ds≤ C2, (26)
for some positive constant C2.
In order to study the first factor on the right hand side of the inequality (24), we have to control the “delta” approximation, ∂Pε,δ
∂x →
∂PBS
∂x , as opposed to the option price approximation, P ε,δ
→ PBS, studied in [8] for
European options, or in [7] for digital-type options.
By pathwise differentiation (see [9] for instance), the chain rule can be applied and we obtain ∂Pε,δ ∂St (t, St, Yt, Zt) = IE⋆ e−r(T −t)I {ST>K} ∂ST ∂St | St , Yt, Zt . At time t = 0, e−rT∂ST ∂S0 = eR0TσtdWt(0)∗− 1 2 RT 0 σ 2 tdt (27)
gives an exponential martingale, and therefore one can construct a IP⋆-equivalent probability measure ˜P by Girsanov Theorem. As a result, the delta ∂P∂Sε,δt (t, St, Yt, Zt) has a probabilistic representation under the new
measure ˜P corresponding to the digital-type option ∂Pε,δ
∂St
(t, St, Yt, Zt) = ˜EI{ST>K}| St, Yt, Zt ,
where the dynamics of St becomes
dSt = r + f2(Yt, Zt) Stdt + σtStd ˜Wt,
with ˜W being a standard Brownian motion under ˜P . The dynamics of Yt and Zt remain the same as in (1)
because we have assumed here zero correlations between Brownian motions. Then one can apply the accuracy result in [7] for digital options to claim that
E˜I{ST>K}| St, Yt, Zt − ¯EI{ ¯ST>K}| ¯St= St, Zt ≤ C3(Yt) max{ √ ε,√δ}, where the constant C3 may depend on Yt, and the “homogenized” stock price ¯Stsatisfies
with ¯Wt being a standard Brownian motion [6]. In fact, the homogenized approximation ¯EI{ ¯ST>K}| ¯St, Zt
is a probabilistic representation of the homogenized “delta”, ∂PBS
∂x . Consequently, we obtain the accuracy result
for delta approximation: ∂Pε,δ ∂x − ∂PBS ∂x (t, St, Yt, Zt) ≤ C 3(Yt) max{√ε, √ δ}.
The existence of moments of Yt ensures the existence of the fourth moment of C3(Yt), and therefore the first
factor on the right hand side of (24) is bounded by v u u tIE⋆ ( Z T 0 ∂Pε,δ ∂x − ∂PBS ∂x 4 (s, Ss, Ys, Zs)ds ) ≤ C4max{ε, δ}. (28)
for some positive constant C4. From (24), (28) and (26), we conclude that
IE⋆ ( Z T 0 e−2rs ∂P ε,δ ∂x − ∂PBS ∂x 2 (s, Ss, Ys, Zs)f2(Ys, Zs)Ss2ds ) ≤ C1max{ε, δ}
for some constant C1.
Lemma A.2. Under the assumptions of Theorem 1.1, for any fixed initial state (0, x, y, z), there exists ε a positive constant C such that for ε≤ 1 and δ ≤ 1, one has
Z T 0 e−2rs ∂P ε,δ ∂y 2 (s, Ss, Ys, Zs)g12(Ys)ds≤ C ε2
Proof: Conditioning on the path of volatility process and by iterative expectations, the price of a European option can be expressed as
Pε,δ(t, x, y, z) = IE⋆t,x,y,z n IE⋆ne−r(T −t)(ST − K)+| σs, t≤ s ≤ T oo = IE⋆t,x,y,z n PBS t, x; K, T ;pσ2o, (29)
where the realized variance is denoted by σ2:
σ2= 1
T − t Z T
t
f (Ys, Zs)2ds. (30)
Taking a pathwise derivative for Pε,δ[9] with respect to the fast varying variable y, we deduce by the chain rule
∂Pε,δ ∂y (t, x, y, z) = IE ⋆ t,x,y,z ( ∂PBS ∂σ t, x; K, T ; q σ2(y, z) ∂ p σ2 ∂y ) . (31)
Inside of the expectation the first derivative, known as Vega, ∂PBS ∂σ = xe−d2 1/2√T − t √ 2π ,
with d1=
log(x/K)+(r+1 2σ
2 )(T −t)
σ√T −t , is uniformly bounded in σ. Using the chain rule one obtains
∂pσ2 ∂y = 1 (T − t)pσ2 Z T t ∂f ∂y(Ys, Zs) ∂Ys ∂y + ∂f ∂z(Ys, Zs) ∂Zs ∂y f (Ys, Zs)ds. (32)
In order to control the growth rate of ∂Ys
∂y and
∂Zs
∂y we consider their dynamics:
d ds ∂Ys ∂y ∂Zs ∂y ! = − 1 ε+ ν1√√2 ε ∂Λ1 ∂y (Ys, Zs) ν1√√2 ε ∂Λ1 ∂z (Ys, Zs) √ δν2∂Λ∂y2(Ys, Zs) −δ + √ δν2√2∂Λ∂z2(Ys, Zs) ! ∂Y s ∂y ∂Zs ∂y ! (33)
with the initial condition∂Yt ∂y,
∂Zt ∂y
= (1, 0). Rescaling the system (33) by defining ˜Yε
s = Ysε and ˜Zsε= Zsε, we deduce d ds ∂ ˜Yε s ∂y ∂ ˜Zε s ∂y ! = −1 0 0 0 ∂ ˜Ysε ∂y ∂ ˜Zε s ∂y ! + √ε ν1 √ 2∂ ˜Λ1 ∂y ( ˜Ysε, ˜Zsε) ∂ ˜∂zΛ1( ˜Ysε, ˜Zsε) √ δεν2∂ ˜∂yΛ2( ˜Ysε, ˜Zsε) −δ √ε +√ δεν2√2∂ ˜∂zΛ2( ˜Ysε, ˜Zsε) ! ∂ ˜Yε s ∂y ∂ ˜Zε s ∂y ! .
The functions ˜Λ1 and ˜Λ2 are defined according to the rescaling and they are smooth and bounded as Λ’s. By
a classical stability result [2], we obtain |∂Ys
∂y| < C5e−(s−t)/ε and | ∂Zs
∂y | < C6δε for some constants C5 and C6.
Applying these estimates to (32) and by the smooth boundness of f , we obtain ∂pσ2
∂y ≤ Cε
for some C, and consequently a similar bound for ∂P∂yε,δ(t, x, y, z) in (31). Finally, as g1 = ν1√2, Lemma A.2
follows.
Lemma A.3. Under the assumptions of Theorem 1.1, for any fixed initial state (0, x, y, z), there exists a positive constant C such that for ε≤ 1 and δ ≤ 1, one has
Z T 0 e−2rs ∂P ε,δ ∂z 2 (s, Ss, Ys, Zs)g22(Zs)ds≤ C.
Proof: The proof is similar to Lemma A.2.
From the bounds in Lemma A.1, A.2 and A.3, we deduce Theorem 1.1.
Acknowledgment
The authors would like to thank Bernard Lapeyre for his helpful comments and suggestions on an earlier version of this work.
References
[1] G. Barone-Adesi and R. E. Whaley, “Efficient Analytic Approximation of American Option Values,” The Journal of Finance, Vol. XLII, No. 2, June 1987.
[3] E. Clement, D. Lamberton, P. Protter, “An Analysis of a Least Square Regression Method for American Option Pricing,” Finance and Stochastics 6:449-471, 2002.
[4] J.-P. Fouque and C.-H. Han, “A Control Variate Method to Evaluate Option Prices under Multi-Factor Stochastic Volatility Models,” submitted, 2004.
[5] J.-P. Fouque and C.-H. Han, “Variance Reduction for Monte Carlo Methods to Evaluate Option Prices under Multi-factor Stochastic Volatility Models,” Quantitative Finance 4(5), October 2004 (597-606).
[6] J.-P. Fouque, G. Papanicolaou, and R. Sircar, “Derivatives in Financial Markets with Stochastic Volatility,” Cambridge Uni-versity Press, 2000.
[7] J.-P. Fouque, R. Sircar, and K. Solna, “Stochastic Volatility Effects on Defaultable Bonds,” Applied Mathematical Finance, to appear in 2006.
[8] J.-P. Fouque, G. Papanicolaou, R. Sircar, and K. Solna, “Multiscale Stochastic Volatility Asymptotics,” SIAM Journal on Multiscale Modeling and Simulation 2(1), 2003 (22-42).
[9] P. Glasserman, “Monte Carlo Methods in Financial Engineering,” Springer Verlag, 2003.
[10] F. Longstaff and E. Schwartz, “Valuing American Options by Simulation: A Simple Least-Squares Approach,” Review of Financial Studies 14: 113-147, 2001.
[11] B. Oksendal, “Stochastic Differential Equations: An introduction with Applications,” Universitext, 5th edn, Springer, 1998. [12] P. Wilmott , S. Howison and J. Dewynne, ”Mathematics of Financial Derivatives: A Student Introduction”, Cambridge
