Pricing High-Dimensional Bermudan Options
using Variance-Reduced Monte Carlo Methods
Peter Hepperger
∗
We present a numerical method for pricing Bermudan options depend-
ing on a large number of underlyings. The asset prices are modeled with
exponential time-inhomogeneous jump-diffusion processes. We improve
the least-squares Monte Carlo method proposed by Longstaff and Schwartz
introducing an efficient variance reduction scheme. A control variable is
obtained from a low-dimensional approximation of the multivariate Bermu-
dan option. To this end, we adapt a model reduction method called proper
orthogonal decomposition (POD), which is closely related to principal com-
ponent analysis, to the case of Bermudan options. Our goal is to make use
of the correlation structure of the assets in an optimal way. We compute
the expectation of the control variable by either solving a low-dimensional
partial integro-differential equation or by applying Fourier methods. The
POD approximation can also be used as a candidate for the minimizing
martingale in the dual pricing approach suggested by Rogers. We evaluate
both approaches in numerical experiments.
Key Words Bermudan options, dimension reduction, proper orthogonal de-
composition, regression-based Monte Carlo, Fourier methods
1 Introduction
The present article is concerned with numerical pricing of multidimensional Bermu-
dan options. We call an option multidimensional if it is written on more than one
underlying. Important examples on the stock market include index options on an av-
erage, and basket options on the minimum or maximum of a set of assets. Derivatives
on other markets depending on several prices or driving factors also fall into this cat-
egory. We will consider Bermudan options with a finite number of exercise dates. By
choosing this number sufficiently high, the method presented here can could be used
∗
Zentrum Mathematik, Technische Universität München, 85748 Garching bei München, Germany
(hepperger@ma.tum.de). Supported by the International Graduate School of Science and Engineering
(IGSSE) of Technische Universität München.
1
to approximate American options with a continuum of exercise possibilities. Note
that in fact whenever Monte Carlo (MC) simulation is used for pricing American op-
tions, this actually amounts to an approximation with a Bermudan option on the MC
time discretization grid. The approximation error may increase, however, for a larger
number of exercise points.
We assume the underlying prices to be driven by a multivariate time-inhomogeneous
jump-diffusion process. Theoretically, the same numerical pricing methods as in the
one-dimensional case can be employed here. These include partial integro-differential
equations (PIDEs) [7, 8] and Fourier transform methods [4, 5, 18]. However, they both
suffer from the curse of dimensionality which means that in practice they cannot be ap-
plied directly when the number of underlyings is large. MC methods are a feasible
alternative, whose complexity does not increase exponentially in the dimension. MC
simulations for Bermudan options are often based on dynamic programming princi-
ples. The Snell envelope is obtained through backward recursion. The conditional
expectation of future cashflows at each exercise point can be approximated using re-
gression, a method first introduced by Longstaff and Schwartz [19, 3, 15]. The option is
exercised if and only if the current intrinsic value is larger than this expectation. These
methods yield an approximation from below for the fair price. A different approach
uses a dual representation of the price, which allows for the computation of an upper
bound [20]. For an overview of Bermudan MC pricing see [16] and the references
therein.
The major drawback of all MC algorithms is their comparatively slow convergence
rate. There are several techniques, subsumed under the term variance reduction, which
can help obtaining more accurate results with fewer simulated paths. These include
antithetic variables, importance sampling, and control variables [1, 10, 11]. We will
focus on the latter and present an improved Longstaff-Schwartz algorithm using a
low-dimensional approximation of the option price as the control variable.
The dimension reduction relies on an orthogonal projection method called proper
orthogonal decomposition (POD) [17]. Similar to principal component analysis (PCA),
a small set of orthonormal vectors is found, which minimizes the projection error.
The approximation is hence optimal (in the L
2
-sense). A similar projection method
has already been successfully applied to European option pricing [13]. We extend
the concept to Bermudan options. In particular, we no longer rely on direct, accurate
computation of the option price via POD, but rather use a fast, coarse approximation
for variance reduction purposes.
An estimate for the approximation error is derived. The expectation of the low-
dimensional POD approximation can be computed efficiently with PIDE or Fourier
algorithms. Once this is done, the solution can be used for two purposes: first, it
can serve as a control variable. As the approximation is highly correlated with the
full-dimensional price process, this results in a substantially decreased variance of the
modified MC estimator. Second, the POD solution is a candidate for the minimizing
martingale in Rogers’ [20] dual MC method. We discuss both approaches. The dimen-
sion of the projection can be freely chosen. This allows for a trade-off between reduced
variance and increased effort for the computation of the low-dimensional expectation.
2
The effectiveness of the POD depends on the correlation structure of the underlying
assets. High correlation allows for a more efficient decrease of variance with fewer
POD components. We investigate the performance of the POD variance-reduction
in numerical experiments. Using different types of options, basket sizes, correlation
parameters, and numbers of exercise dates, we show that the overall computational
time can as a rule be reduced by at least 50%, often even by more than 80%.
The paper is organized as follows. In Section 2, we present the multivariate jump-
diffusion model and formulate the Bermudan option pricing problem. The dimension
reduction methods and the corresponding convergence results are derived in Section 3.
Next, we describe the algorithms for variance-reduced Bermudan MC and the dual
MC approach in detail in Section 4. Section 5 contains the numerical experiments.
We describe the test settings and analyze the computational results. Finally, Section 6
gives a short conclusion and summary of the article.
2 Bermudan Basket Options
In this section, the market model used throughout the paper is introduced. We de-
fine the driving stochastic jump-diffusion process, declare assumptions concerning the
coefficients, and state the Bermudan option pricing problem.
2.1 The asset price process
We consider a Bermudan option depending on n assets. The terminal date of maturity
is T > 0 (last exercise date). The multivariate asset price process is denoted
(2.1) S
t
:
=
S
1
( t), . . . , S
n
( t)
:
=
S
1
(0) e
X
1
(t)
, . . . , S
n
(0) e
X
n
(t)
∈ R
n
, t ∈ [0, T].
For each i = 1, . . . , n, the price S
i
of the ith asset is modeled as the product of its initial
value S
i
(0) > 0 and the ordinary exponential of a time-inhomogeneous jump-diffusion
process X
i
, given by
X
t
:
=
X
1
( t), . . . , X
n
( t)
:
=
Z
t
0
γ
s
ds +
Z
t
0
σ
s
dW( s) +
Z
t
0
Z
H
η
s
ξ
e
M(dξ, ds) ∈ R
n
, t ∈ [0, T].
(2.2)
The diffusion part is driven by an R
n
-valued Brownian motion W. The jumps are
characterized by
e
M, the compensated random measure of an R
n
-valued compound
Poisson process
J
t
=
N
t
∑
i=1
Y
i
, t ≥ 0,
which is independent of W. Here, N denotes a Poisson process with intensity λ and
Y
i
∼ P
Y
(i = 1, 2, . . .) are iid on R
n
(and independent of N). The corresponding Lévy
3
measure is denoted by ν = λP
Y
. We assume the drift γ
:
[0, T] → R
n
, the volatility
σ
:
[0, T] → R
n×n
, and the jump integrand η
:
[0, T] → R
n×n
to be deterministic
functions. We make the following assumption concerning the moments of the process.
Assumption 2.1. The second exponential moment of the jump distribution Y exists:
E[e
2
k
Y
k
R
n
] =
Z
R
n
e
2
k
ξ
k
R
n
P
Y
( dξ) < ∞.
We assume further that
Z
T
0
k
γ
t
k
2
R
n
dt < ∞,
Z
T
0
k
σ
t
k
2
R
n×n
dt < ∞, and
k
η
t
k
R
n×n
≤ 1 for a.e. t ∈ [0, T].
The interest rate r is assumed to be constant. In order to avoid a discussion of
possible measure changes, we suppose the model (2.1) to be stated under the pricing
measure. In view of the dimension reduction performed in Section 3, it turns out to
be useful to express the value of the option in terms of the centered process
Z
t
:
= X
t
− E[X
t
], t ∈ [0, T].
Since the asset price process S = S(Z) depends on Z in a deterministic way, this is
nothing more than a simple transform of variables.
2.2 Bermudan Options
A Bermudan option grants the holder the right to exercise at one of N
ex
∈ N admis-
sible dates, which we denote by 0 ≤ t
1
< t
2
< · · · < t
N
ex
= T. Let T (t, T) denote the
set of all stopping times with values in {t
i
|1 ≤ i ≤ N
ex
and t
i
≥ t}. For simplicity, we
assume a constant interest rate r > 0. The discounted value V of a Bermudan option
at time t, given that the option was not yet exercised, is the solution of the optimal
stopping problem
(2.3) V(t, z) = sup
τ∈T (t,T)
E
e
−rτ
g
S(Z
τ
)
Z
t
= z
.
The payoff g
:
R
n
→ R is determined by the asset price process S which can be
expressed in terms of the centered jump-diffusion Z. We make the following assump-
tion.
Assumption 2.2. We assume that there is a constant L
g
such that the payoff function g
satisfies the Lipschitz condition
g(S) − g(
e
S)
≤ L
g
S −
e
S
R
n
for every S,
e
S ∈ R
n
.
Remark 2.3. Assumption 2.2 is considerably weaker than the corresponding assumption in
[13, Ass. 3.4]: it refers to the payoff in terms of S instead of Z. In particular, it is satisfied for
plain vanilla call and put options on weighted averages of the asset prices. The convergence
proofs below (Lemma 3.4 and Theorem 3.5) account for this weaker assumption with additional
technical estimates.
4
This assumption is, in particular, satisfied for index options written on a weighted
average of the assets. An index put has the intrinsic value
g(S) =
K −
n
∑
i=1
w
i
S
i
!
+
,
where K is the strike and w
i
∈ R are constant weights. Other examples include
maximum or minimum options. A maximum or minimum put corresponds to the
payoff
g(S) =
K − max
i=1,...,n
S
i
+
or g(S) =
K − min
i=1,...,n
S
i
+
,
respectively.
The aim when pricing Bermudan options is to find the optimal exercise time for
(2.3). It is well known that this can be done by backward dynamic programming: at
time t = T the value of the option is
(2.4) V(T, z) = e
−rT
g
S(Z
T
)
.
For any previous exercise date t
i
, i = 0, . . . , N
ex
− 1, the value is
(2.5) V(t
i
, z) = max
e
−rt
i
g
S(z)
, E
V(t
i+1
, Z
t
i+1
)
Z
t
i
= z
.
Hence, it is optimal to exercise at time t
i
if and only if the intrinsic value g
S(z)
is
larger than or equal to the expected discounted future cash flow (given the option is
not yet exercised). Computing the conditional expectations
(2.6) E
V(t
i+1
, Z
t
i+1
)
Z
t
i
= z
for every exercise date is the basic challenge. The fair value of the option at time t = 0
is then given by V(0, 0).
3 Dimension Reduction
It is possible to derive a partial integro-differential (PIDE) equation which is satis-
fied by the conditional expectation (2.6) (see, e.g., [13]). Such a differential equation,
however, suffers from the curse of dimensionality. The same holds true for Fourier
transforms. It is hardly possible to apply these methods directly for large values of n,
say n > 10. In this section, we derive a low-dimensional approximation for Bermu-
dan options. To this end, we employ proper orthogonal decomposition (POD), which
makes use of the correlation of the individual assets. The idea is similar to principal
component analysis (PCA): we approximate the centered driving process Z with a
small set of orthonormal vectors.
5