Exact retrospective Monte Carlo computation
of arithmetic average Asian options
Benjamin Jourdain1 and Mohamed Sbai 1
Abstract
Taking advantage of the recent literature on exact simulation algorithms (Beskos et al. [1]) and
unbiased estimation of the expectation of certain functional integrals (Wagner [23], Beskos et al. [2] and
Fearnhead et al. [6]), we apply an exact simulation based technique for pricing continuous arithmetic
average Asian options in the Black & Scholes framework. Unlike existing Monte Carlo methods, we are
no longer prone to the discretization bias resulting from the approximation of continuous time processes
through discrete sampling. Numerical results of simulation studies are presented and variance reduction
problems are considered.
Introduction
Although the Black & Scholes framework is very simple, it is still a challenging task to efficiently price
Asian options. Since we do not know explicitly the distribution of the arithmetic sum of log-normal variables,
there is no closed form solution for the price of an Asian option. By the early nineties, many researchers
attempted to address this problem and hence different approaches were studied including analytic approxima-
tions (see Turnball and Wakeman [20], Vorst [22], Levy [15] and more recently Lord [16]), PDE methods (see
Vecer [21], Rogers and Shi [18], Ingersoll [11], Dubois and Lelievre [5]), Laplace transform inversion methods
(see Geman and Yor [10], Geman and Eydeland [8]) and, of course, Monte Carlo simulation methods (see
Kemna and Vorst [13], Broadie and Glasserman [3], Fu et al. [7]).
Monte Carlo simulation can be computationally expensive because of the usual statistical error. Variance
reduction techniques are then essential to accelerate the convergence (one of the most efficient techniques is
the Kemna&Vorst control variate based on the geometric average). One must also account for the inherent
discretization bias resulting from approximating the continuous average of the stock price with a discrete
one. It is crucial to choose with care the discretization scheme in order to have an accurate solution (see
Lapeyre and Temam [14]). The main contribution of our work is to fully address this last feature by the use,
after a suitable change of variables, of an exact simulation method inspired from the recent work of Beskos
et al. [1, 2] and Fearnhead et al. [6].
In the first part of the paper, we recall the algorithm introduced by Beskos et al. [1] in order to simulate
sample-paths of processes solving one-dimensional stochastic differential equations. By a suitable change of
variables, one may suppose that the diffusion coefficient is equal to one. Then, according to the Girsanov
theorem, one may deal with the drift coefficient by introducing an exponential martingale weight. Because
of the one-dimensional setting, the stochastic integral in this exponential weight is equal to a standard
1Project team Math Fi, CERMICS, Ecole des Ponts, Paristech, supported by the ANR program ADAP’MC.
Postal address : 6-8 av. Blaise Pascal, Cité Descartes, Champs-sur-Marne, 77455 Marne-la-Vallée Cedex 2.
E-mails : jourdain@cermics.enpc.fr and sbai@cermics.enpc.fr
http://arxiv.org/abs/0704.1433v3
integral with respect to the time variable up to the addition of a function of the terminal value of the
path. Under suitable assumptions, conditionally on a Brownian path, an event with probability equal to
the normalized exponential weight can be simulated using a Poisson point process. This allows to accept or
reject this Brownian path as a path solution to the SDE with diffusion coefficient equal to one. In finance,
one is interested in computing expectations rather than exact simulation of the paths. In this perspective,
computation of the exponential importance sampling weight is enough. The entire series expansion of the
exponential function permits to replace this exponential weight by a computable weight with the same
conditional expectation given the Brownian path. This idea was first introduced by Wagner [23, 24, 25, 26]
in a statistical physics context and it was very recently revisited by Beskos et al. [2] and Fearnhead et al.
[6] for the estimation of partially observed diffusions. Some of the assumptions necessary to implement the
exact algorithm of Beskos et al. [1] can then be weakened.
The second part is devoted to the application of these methods to option pricing within the
Black & Scholes framework. Throughout the paper, St = S0 exp
σWt + (r − δ −
represents the
stock price at time t, T the maturity of the option, r the short interest rate, σ the volatility parameter, δ
the dividend rate and (W )t∈[0,T ] denotes a standard Brownian motion on the risk-neutral probability space
(Ω,F ,P). We are interested in computing the price C0 = E
e−rTf
αST + β
of a European
option with pay-off f
αST + β
assumed to be square integrable under the risk neutral measure P.
The constants α and β are two given non-negative parameters.
When α > 0, we remark that, by a change of variables inspired by Rogers and Shi [18], αST + β
has the same law as the solution at time T of a well-chosen one-dimensional stochastic differential equation.
Then it is easy to implement the exact methods previously presented. The case α = 0 of standard Asian
options is more intricate. The previous approach does not work and we propose a new change of variables
which is singular at initial time. It is not possible to implement neither the exact simulation algorithm nor
the method based on the unbiased estimator of Wagner [23] and we propose a pseudo-exact hybrid method
which appears as an extension of the exact simulation algorithm. In both cases, one first replaces the integral
with respect to the time variable in the function f by an integral with respect to time in the exponential
function. Because of the nice properties of this last function, exact computation is possible.
1 Exact Simulation techniques
1.1 The exact simulation method of Beskos et al. [1]
In a recent paper, Beskos et al. [1] proposed an algorithm which allows to simulate exactly the solution
of a 1-dimensional stochastic differential equation. Under some hypotheses, they manage to implement an
acceptance-rejection algorithm over the whole path of the solution, based on recursive simulation of a biased
Brownian motion. Let us briefly recall their methodology. We refer to [1] for the demonstrations and a
detailed presentation.
Consider the stochastic process (ξt)0≤t≤T determined as the solution of a general stochastic differential
equation of the form : {
dξt = b(ξt)dt+ σ(ξt)dWt
ξ0 = ξ ∈ R
where b and σ are scalar functions satisfying the usual Lipschitz and growth conditions with σ non vanishing.
To simplify this equation, Beskos et al. [1] suggest to use the following change of variables : Xt = η(ξt)
where η is a primitive of 1
(η(x) =
Under the additional assumption that 1
is continuously differentiable, one can apply Itô’s lemma to get
dXt = η
′(ξt)dξt +
η′′(ξt) d< ξ, ξ >t
b(ξt)
σ(ξt)
dt+ dWt −
σ′(ξt)
b(η−1(Xt))
σ(η−1(Xt))
′(η−1(Xt))
︸ ︷︷ ︸
a(Xt)
dt+ dWt
So ξt = η
−1(Xt) where (Xt)t is a solution of the stochastic differential equation
dXt = a(Xt)dt+ dWt
X0 = x.
Thus, without loss of generality, one can start from equation (2) instead of (1).
Let us denote by (W xt )t∈[0,T ] the process (Wt+x)t∈[0,T ], by QWx its law and by QX the law of the process
(Xt)t∈[0,T ]. From now on, we will denote by (Yt)t∈[0,T ] the canonical process, that is the coordinate mapping
on the set C([0, T ],R) of real continuous maps on [0, T ] (see Revuz and Yor [17] or Karatzas and Shreve
[12]).
One needs the following assumption to be true
Assumption 1 : Under QWx , the process
Lt = exp
a(Yu)dYu −
a2(Yu)du
is a martingale.
According to Rydberg [19] (see the proof of Proposition 4 where we give his argument on a specific
example), a sufficient condition for this assumption to hold is
-Existence and uniqueness in law of a solution to the SDE (2).
-∀t ∈ [0, T ],
a2(Yu)du <∞, QX and QWx almost surely on C([0, T ],R).
Thanks to this assumption, one can apply the Girsanov theorem to get that QX is absolutely continuous
with respect to QWx and its Radon-Nikodym derivative is equal to
= exp
a(Yt)dYt −
a2(Yt)dt
Consider A the primitive of the drift a, and assume that
Assumption 2 : a is continuously differentiable.
Since, by Itô’s lemma, A(W xT ) = A(x) +
a(W xt )dW
a′(W xt )dt, we have
= exp
A(YT )−A(x)−
a2(Yt) + a
′(Yt)dt
Before setting up an acceptance-rejection algorithm using this Radon-Nikodym derivative, a last step is
needed. To ensure the existence of a density h(u) proportional to exp(A(u) − (u−x)
), it is necessary and
sufficient that the following assumption holds
Assumption 3 : The function u 7→ exp(A(u)− (u−x)
) is integrable.
Finally, let us define a process Zt distributed according to the following law QZ
(W xt )t∈[0,T ]|W xT = y
h(y)dy
where the notation L(.|.) stands for the conditional law. One has
= C exp
a2(Yt) + a
′(Yt)dt
where C is a normalizing constant. At this level, Beskos et al. [1] need another assumption
Assumption 4 : The function φ : x 7→ a
2(x)+a′(x)
is bounded from below.
Therefore, one can find a lower bound k of this function and eventually the Radon-Nikodym derivative
of the change of measure between X and Z takes the form
= Ce−kT exp
φ(Yt)− k dt
The idea behind the exact algorithm is the following : suppose that one is able to simulate a continuous
path Zt(ω) distributed according to QZ and let M(ω) be an upper bound of the mapping t 7→ φ(Zt(ω))− k.
Let N be an independent random variable which follows the Poisson distribution with parameter TM(ω) and
let (Ui, Vi)i=1...N be a sequence of independent random variables uniformly distributed on [0, T ]× [0,M(ω)].
Then, the number of points (Ui, Vi) which fall below the graph {(t, φ(Zt(ω))− k); t ∈ [0, T ]} is equal to zero
with probability exp
φ(Zt(ω))− k dt
. Actually, simulating the whole path (Zt)t∈[0,T ] is not necessary.
It is sufficient to determine an upper bound for φ(Zt)− k since, as pointed out by the authors, it is possible
to simulate recursively a Brownian motion on a bounded time interval by first simulating its endpoint, then
simulating its minimum or its maximum and finally simulating the other points2. For this reason, one needs
the following assumption for the algorithm to be feasible :
Assumption 5 : Either lim sup
φ(u) < +∞ or lim sup
φ(u) < +∞.
Suppose for example that lim sup
φ(u) < +∞. The exact algorithm of Bekos et al. [1] then takes the
following form :
Algorithm 1
1. Draw the ending point ZT of the process Z with respect to the density h.
2. Simulate the minimum m of the process Z given ZT .
3. Fix an upper bound M(m) = sup{φ(u)− k;u ≥ m} for the mapping t 7→ φ(Zt)− k.
4. Draw N according to the Poisson distribution with parameter TM(m) and draw (Ui, Vi)i=1...N , a
sequence of independent variables uniformly distributed on [0, T ]× [0,M(m)].
5. Fill in the path of Z at the remaining times (Ui)i=1...N .
2In their paper, the authors explain how to do such a decomposition of the Brownian path.
6. Evaluate the number of points (Vi)i=1...N such that Vi ≤ φ(ZUi)− k.
If it is equal to zero, then return the simulated path Z.
Else, return to step 1.
This algorithm gives exact skeletons of the process X , solution of the SDE (2). Once accepted, a path can
be further recursively simulated at additional times without any other acceptance/rejection criteria. We also
point out that the same technique can be generalized by replacing the Brownian motion in the law of the
proposal Z by any process that one is able to simulate recursively by first simulating its ending point, its
minimum/maximum and then the other points. Also, the extension of the algorithm to the inhomogeneous
case, where the drift coefficient a in (2), and therefore the function φ, depend on the time variable t, is
straightforward given that the assumptions presented above are appropriately modified.
1.2 The unbiased estimator (U.E)
In finance, the pricing of contingent claims often comes down to the problem of computing an expectation
of the form
C0 = E (f(XT )) (3)
where X is a solution of the SDE (2) and f is a scalar function such that f(XT ) is square integrable. In a
simulation based approach, one is usually unable to exhibit an explicit solution of this SDE and will therefore
resort to numerical discretization schemes, such as the Euler or Milstein schemes, which introduce a bias.
Of course, the exact algorithm presented above avoids this bias. Here, we are going to present a technique
which permits to compute exactly the expectation (3) while assumptions 4 and 5 on the function a
which
appears in the Radon-Nikodym derivative are relaxed.
Using the previous results and notations, we get, under the assumptions 1 and 2, that
C0 = E
f(W xT ) exp
A(W xT )−A(x) −
a2(W xt ) + a
′(W xt )dt
. (4)
In order to implement an importance sampling method, let us introduce a positive density ρ on the real
line and a process (Zt)t∈[0,T ] distributed according to the following law QZ
(W xt )t∈[0,T ]|W xT = y
ρ(y)dy.
By (4), one has
C0 = E
ψ(ZT ) exp
φ(Zt)dt
where ψ : z 7→ f(z) e
A(z)−A(x)−
(z−x)2
2πρ(z)
and φ : z 7→ a
2(z)+a′(z)
. We do not impose ρ to be equal to the
density h of the previous section. It is a free parameter chosen in such a way that it reduces the variance of
the simulation.
In his first paper, Wagner [23] constructs an unbiased estimator of the expectation (5) when ψ is a
constant, (Zt)t∈[0,T ] is an R
d−valued Markov process with known transition function and φ is a measurable
function such that E
|φ(Zt)|dt
< +∞. His main idea is to expand the exponential term in a power
series, then, using the transition function of the underlying Markov process and symmetry arguments, he
constructs a signed measure ν on the space Y =
n=0([0, T ]× Rd)n+1 such that the expectation at hand is
equal to ν(Y). Consequently, any probability measure µ on Y that is absolutely continuous with respect to
ν gives rise to an unbiased estimator ζ defined on (Y, µ) via ζ(y) = dν
(y). In practice, a suitable way to
construct such an estimator is to use a Markov chain with an absorbing state. Wagner also discusses variance
reduction techniques, specially importance sampling and a shift procedure consisting on adding a constant c
to the integrand φ and then multiplying by the factor e−cT in order to get the right expectation. Wagner [25]
extends the class of unbiased estimators by perturbing the integrand φ by a suitably chosen function φ0 and
then using mixed integration formulas representation. Very recently, Beskos et al. [2] obtained a simplified
unbiased estimator for (5), termed Poisson estimator, using Wagner’s idea of expanding the exponential in
a power series and his shift procedure. To be specific, the Poisson estimator writes
ψ(ZT )e
cpT−cT
c− φ(ZVi)
where N is a Poisson random variable with parameter cP and (Vi)i is a sequence of independent random
variables uniformly distributed on [0, T ]. Fearnhead et al. [6] generalized this estimator allowing c and cP to
depend on Z and N to be distributed according to any positive probability distribution on N. They termed
the new estimator the generalized Poisson estimator. We introduce a new degree of freedom by allowing the
sequence (Vi)i to be distributed according to any positive density on [0, T ]. This gives rise to the following
unbiased estimator for (5) :
Lemma 1 — Let pZ and qZ denote respectively a positive probability measure on N and a positive probability
density on [0, T ]. Let N be distributed according to pZ and (Vi)i∈N∗ be a sequence of independent random
variables identically distributed according to the density qZ , both independent from each other conditionally
on the process (Zt)t∈[0,T ]. Let cZ be a real number which may depend on Z. Assume that
|ψ(ZT )|e−cZT exp
|cZ − φ(Zt)|dt
ψ(ZT )e
−cZT 1
pZ(N)N !
cZ − φ(ZVi)
qZ(Vi)
is an unbiased estimator of C0.
Proof : The result follows from
ψ(ZT )e
−cZT 1
pZ(N)N !
cZ − φ(ZVi)
qZ(Vi)
∣∣∣(Zt)t∈[0,T ]
= ψ(ZT )e
cZ − φ(Zt)dt
pZ(n)n!
pZ(n)
= ψ(ZT ) exp
φ(Zt)dt
Using (7), one is now able to compute the expectation at hand by a simple Monte Carlo simulation. The
practical choice of pZ and qZ conditionally on Z is studied in the appendix 4.1.
As pointed out in Fearnhead et al. [6], this method is an extension of the exact algorithm method since,
under assumptions 3, 4 and 5, the reinforced integrability assumption of Lemma 1 is always satisfied.
Indeed, suppose for example that lim sup
φ(u) < +∞ and let k be a lower bound of φ, mZ be the
minimum of the process Z and MZ an upper bound of {φ(u)− k, u ≥ mZ}. Then, taking cZ = MZ + k in
Lemma 1 ensures the integrability condition :
|ψ(ZT )|e−(MZ+k)T e
|MZ+k−φ(Zt)|dt
|ψ(ZT )|e−(MZ+k)T e
MZ+k−φ(Zt)dt
|ψ(ZT )|e−
φ(Zt)dt
and hence, one is allowed to write that
C0 = E
ψ(ZT )e
−(MZ+k)T 1
pZ(N)N !
MZ + k − φ(ZVi)
qZ(Vi)
Better still, the random variable ψ(ZT )e
−(MZ+k)T 1
pZ(N)N !
MZ+k−φ(ZVi )
qZ(Vi)
is square integrable when
pZ is the Poisson distribution with parameterMZT + k and qZ is the uniform distribution on [0, T ] since we
have then
ψ(ZT )e
−(MZ+k)T 1
pZ(N)N !
MZ + k − φ(ZVi)
qZ(Vi)
 = E
ψ2(ZT )
1− φ(ZVi)
MZ + k
ψ2(ZT )
The last inequality follows from the square integrability of f : whenever one is able to simulate from the
density h, introduced in the exact algorithm, by doing rejection sampling, there exists a density ρ such that
ψ, which is equal to f(ZT )
h(ZT )
ρ(ZT )
up to a constant factor, is dominated by f and so is square integrable.
The square integrability property is very important in that we use a Monte Carlo method. We see that,
whenever the exact algorithm is feasible, the unbiased estimator of lemma 1 is a simulable square integrable
random variable, at least for the previous choice of pZ and qZ .
Remark 2 — One can derive two estimators of C0 from the result of Lemma 1 :
f(ZiT )
T )−A(x)−
2πρ(ZiT )
e−cZT
pZ(N i)N i!
cZ − φ(ZiV i
f(ZiT )
T )−A(x)−
2πρ(ZiT )
pZ(N i)N i!
cZ − φ(ZiV i
T )−A(x)−
2πρ(ZiT )
pZ(N i)N i!
cZ − φ(ZiV i
2 Application : the pricing of continuous Asian options
In the Black & Scholes model, the stock price is the solution of the following SDE under the risk-neutral
measure P
= (r − δ)dt+ σdWt (8)
where all the parameters are constant : r is the short interest rate, δ is the dividend rate and σ is the
volatility.
Throughout, we denote γ = r − δ − σ
. The path-wise unique solution of (8) is
St = S0 exp(σWt + γt) .
We consider an option with pay-off of the form
αST + β
where f is a given function such that E
αST + β
< ∞, T is the maturity of the option and
α, β are two given non negative parameters3. Note that for α = 0, this is the pay-off of a standard continuous
Asian option.
The fundamental theorem of arbitrage-free pricing ensures that the price of the option under consideration
C0 = E
e−rTf
αST + β
At first sight, the problem seems to involve two variables : the stock price and the integral of the stock
price with respect to time. Dealing with the PDE associated with Asian option pricing, Rogers and Rogers
and Shi [18] used a suitable change of variables to reduce the spatial dimension of the problem to one. We
are going to use a similar idea.
αS0 + βS0
e−σWu−γudu
eσWt+γt.
We have that
ξt = αS0e
σWt+γt + βS0
eσ(Wt−Wu)+γ(t−u)du
= αS0e
σBt+γt + βS0
eσBs+γsds
where we set Bs = Wt −Wt−s, ∀s ∈ [0, t]. Clearly, (Bs)s∈[0,t] is a Brownian motion and thus the following
lemma holds
Lemma 3 — ∀t ∈ [0, T ], ξt and αSt + β
Sudu have the same law.
As a consequence
C0 = E
e−rT f(ξT )
By applying Itô’s lemma, we verify that the process (ξt)t≥0 is a positive solution of the following 1-dimensional
stochastic differential equation for which path-wise uniqueness holds
dξt = βS0dt+ ξt(σdWt + (γ +
ξ0 = αS0.
We are thus able to value C0 by Monte Carlo simulation without resorting to discretization schemes using
one of the exact simulation techniques described in the previous section. In the case α = 0, one has to deal
with the fact that ξt starts from zero which is the reason why we distinguish two cases.
3The underlying of this option is a weighted average of the stock price at maturity and the running average of the stock
price until maturity with respective weights α and βT .
2.1 The case α 6= 0
We are going to apply both the exact algorithm of Beskos et al. [1] and the method based on the unbiased
estimator of lemma 1.
We make the following change of variables to have a diffusion coefficient equal to 1 :
log(ξt)
dXt = (
+ βS0
e−σXt)dt+ dWt
X0 = x with x =
log(αS0)
C0 = E
e−rTf(eσXT )
The following proposition ensures that assumption 1 is satisfied.
Proposition 4 — The process (Lt)t∈[0,T ] defined by
Lt = exp
e−σYt) dYt −
e−σYt)2dt
is a martingale under QWx .
Proof : Under QWx , (Lt)t∈[0,T ] is clearly a non-negative local martingale and hence a super-martingale.
Then, it is a true martingale if and only if E
(LT ) = 1.
Checking the classical Novikov’s or Kamazaki’s criteria is not straightforward. Instead, we are going to
use the approach developed by Rydberg [19] (see also Wong and Heyde [27]) who takes advantage of the link
between explosions of SDEs and the martingale property of stochastic exponentials.
Let us define the following stopping times :
τn(Y ) = inf
t ∈ R+ such that
e−σYu
du ≥ n
with the convention inf{∅} = +∞.
The stopped process (Lt∧τn(Y ))t∈[0,T ] is a true martingale under QWx since Novikov’s condition is fulfilled.
According to the Girsanov theorem, one can define a new probability measure QnX , which is absolutely
continuous with respect to QWx , by its Radon-Nikodym derivative
= LT∧τn(Y ).
Hence
1{τn(Y )>T}
1{τn(Y )>T}LT∧τn(Y )
Since (τn(Y ))n∈N is a non decreasing sequence, we can pass to the limit in the right hand side We get
X (τn(Y ) > T ) = EQWx
1{τ∞(Y )>T}LT∧τ∞(Y )
where τ∞(Y ) denotes the limit of the non decreasing sequence (τn(Y ))n∈N.
Under QWx , (Yt)t∈[0,T ] has the same law as a Brownian motion starting from x so τ∞(Y ) = +∞ ,QWx
almost surely, and consequently
= lim
X (τn(Y ) > T ) .
On the other hand, the Girsanov theorem implies that, under QnX , (Yt)t∈[0,T∧τn(Y )] solves a SDE of the
form (11). To conclude the proof, it is sufficient to check that trajectorial uniqueness holds for this SDE.
Indeed, the law of (Yt)t∈[0,T∧τn(Y )] under Q
X is the same as the law of (Yt)t∈[0,T∧τn(Y )] under QX . Hence
X (τn(Y ) > T ) = QX (τn(Y ) > T ) −→
QX (τ∞(Y ) > T ) .
Clearly,
+ βS0
e−σYu
du <∞, QX almost surely, so
= QX (τ∞(Y ) > T ) = 1
as required.
In order to check trajectorial uniqueness for the SDE (11), we consider two solutions X1 and X2. We
have that
d(X1t −X2t ) =
t − e−σX
dt ⇒ d|X1t −X2t | =
sign(X1t −X2t )
t − e−σX
|X1t −X2t | =
sign(X1s −X2s )
t − e−σX
ds ≤ 0.
The last inequality follows from the fact that x 7→ e−σx is a decreasing function. Finally, almost surely,
∀t ≥ 0, X1t = X2t which leads to strong uniqueness. ✷
Consequently, thanks to the Girsanov theorem, we have
= exp
e−σYt)
︸ ︷︷ ︸
a(Yt)
dYt −
e−σYt)2dt
 . (12)
Set A(u) =
a(x)dx = γ
u+ βS0
(1− e−σu). Then
= exp
A(YT )−A(x)−
a2(Yt) + a
′(Yt)dt
The function u 7→ exp
A(u)− (u−Y0)
= exp
u+ βS0
(1− e−σu)− (u−Y0)
is clearly integrable so we
can define a new process (Zt)t∈[0,T ] distributed according to the following law QZ
(Wt)t∈[0,T ]|WT = y
h(y)dy
where the probability density h is of the form
h(u) = C exp
A(u)− (u− Y0)
with C a normalizing constant. (13)
Remark 5 — Simulating from this probability distribution is not difficult (see the appendix 4.2 for an
appropriate method of acceptance/rejection sampling).
We have
= C exp
(a2(Yt) + a
′(Yt))dt
Set φ(x) =
a2(x)+a′(x)
e−σx)2−βS0e−σx
. A direct calculation gives
φ(x) =
if 2γ ≥ σ2
log( 2βS0
σ2−2γ )
otherwise.
Set k = infx∈R φ(x). Finally, we get
= Ce−kT exp
φ(Yt)− k dt
We check that
φ(x) =
φ(x) = +∞.
Hence we can apply the algorithm 1 to simulate exactly XT and compute C0 = E
e−rTf(eσXT )
Monte Carlo. On the other hand, using (12) we get
C0 = E
e−rTf(eσW
T ) exp
A(W xT )−A(x) −
a2(W xt ) + a
′(W xt )dt
and we can also use the unbiased estimator presented in the previous section to compute this expectation.
Remark 6 — We also applied the exact algorithm based on a geometric Brownian motion instead of the
standard Brownian motion which seems more intuitive given the form of the SDE (10). The algorithm is
feasible because we can simulate recursively a drifted Brownian motion and therefore a geometric Brownian
motion by an exponential change of variables. The results we obtained were not different from the first
method.
2.1.1 Numerical computation
For numerical tests, we consider the case
f(x) = (x−K)+
which corresponds to the European call option with strikeK. Using the exact simulation algorithm presented
above, we can simulate the underlying αST + β
Stdt at maturity (see Figure 1). Then, all we have to
do is a simple Monte Carlo method to get the price of the option under consideration. Using the unbiased
estimator, we get
C0 = E
e−rT (eσZT −K)+
eA(ZT )−A(x)−
(ZT −x)
2πρ(ZT )
e−(MZ+k)T
pZ(N)N !
MZ + k − φ(ZVi)
q(Vi)
where (Zt)t∈[0,T ], ρ,MZ , k, pZ and qZ are defined as in section 1.2. In order to ensure square integrability,
we choose pZ to be a Poisson distribution with parameterMZT +k and qZ to be the uniform distribution on
[0, T ]. For the density ρ, a good choice is to consider the density that we use to simulate from the distribution
h by rejection sampling.
We test these exact methods against a standard discretization scheme with the variance reduction tech-
nique of Kemna and Vorst [13]. As pointed out by Lapeyre and Temam [14], the discretization of the integral
by a simple Riemannian sum is not efficient. Instead, we use the trapezoidal discretization. In the sequel, we
will denote this method by Trap+KV. The table 1 gives the results we obtained for the following arbitrary
set of parameters : S0 = 100, K = 100, r = 0.05, σ = 0.3, δ = 0, T = 1, α = 0.6 and β = 0.4. The
computation has been made on a computer with a 2.8 Ghz Intel Penthium 4 processor. We intentionally
choose a large number of simulations in order to show the influence of the number of time steps when using
a discretization scheme.
10 53 97 140 183 227 270 313 357 400
0.000
0.002
0.004
0.006
0.008
0.010
0.012
0.014
0.016
0.018
PSfrag replacements
αST + β
Exact Simulation of the underlying : αST + β
Figure 1: Histogram of 105 independent realizations of αST + β
Stdt for α = 0.6 and β = 0.4 compared
with the lognormal distribution of ST .
Method M N Acceptance rate Price C.I at 95% CPU
10 11.46 [11.43, 11.48] 5 s
Trap+KV 20 106 - 11.46 [11.43, 11.49] 9 s
50 11.47 [11.44, 11.5] 21 s
Exact Simulation - 106 24% 11.46 [11.43, 11.5] 81 s
U.E (cP =MZ , cZ =MZ + k) - 10
6 - 11.46 [11.43, 11.49] 17 s
U.E (cP = cZ = 1/T ) - 10
6 - 11.46 [11.43, 11.49] 6 s
Table 1: Price of the option (9) using a standard discretization technique and exact simulation methods.
Empirical evidence shows that the exact simulation method is quite slow. This is mainly due to the
fact that the rejection algorithm has a little acceptance rate (24% according to table 1). Using a geometric
Brownian motion instead of a standard Brownian motion did not improve the results. Also, simulating
recursively a Brownian path conditionally on its terminal value and its minimum is time consuming.
The unbiased estimator is more efficient, especially when we can avoid the recursive simulation of the
Brownian path. To do so, we choose for pZ a Poisson distribution with mean cPT where cP is a free
parameter. If we assume that the integrability condition in lemma 1 holds, then we can write that
C0 = E
e−rT (eσZT −K)+
eA(ZT )−A(x)−
(ZT −x)
2πρ(ZT )
ecPT−cZT
cZ − φ(ZVi)
Regarding the dependence of the exact simulation method with respect to the parameters α and β, it
is intuitive that whenever α >> β, the method performs well since the logarithm of the underlying is not
far from the logarithm of the geometric Brownian motion on which we do rejection-sampling. The table 2
confirms this intuition. We see that we cannot apply the algorithm for small values of α and then let α→ 0
to treat the case α = 0.
α + β
0.3 0.4 0.5 0.6 0.7
Acceptance Rate 0.003% 0.47% 5.66% 24.43% 53.85%
Table 2: Influence of the parameter α
on the acceptance rate of the exact algorithm.
2.2 Standard Asian options : the case α = 0 and β > 0
A standard Asian option is a European option on the average of the stock price over a determined period
until maturity. An Asian call, for example, has a pay-off of the form ( 1
Sudu−K)+. With our previous
notations, it corresponds to the case α = 0, β = 1
and f(x) = (x−K)+.
The change of variables we used above is no longer suitable because it starts from zero when α = 0.
Instead, we consider the following new definition of the process ξ
eσ(Wt−Wu)+γ(t−u)du
ξ0 = S0.
Obviously, the two variables ξT and
Sudu have the same law. Hence, the price of the Asian option
becomes
C0 = E
e−rT f
e−rTf(ξT )
Remark 7 — The pricing of floating strike Asian options is also straightforward using this method. It is
even more natural to consider these options since it unveils the appropriate change of variables as we shall
see below.
Let us consider a floating strike Asian call for example. We have to compute
C0 = E
Sudu− ST
Using S̃t = Ste
δt as a numéraire (see the seminal paper of Geman et al. [9]), we immediately obtain that
C0 = EP
−δT ( 1
du− 1
where P
is the probability measure associated to the numéraire S̃t. It is defined by its Radon-Nikodym
derivative
= eσWT−
Under P
, the process Bt =Wt − σt is a Brownian motion and we can write that
C0 = EP
−δT ( 1
eσ(Bu−BT )+(r−δ+
)(u−T )du − 1
−δT ( 1
eσ(Wu−WT )+(r−δ+
)(u−T )du− 1
ξT − S0
where ξt is the process defined by (14) but with γ = r − δ + σ
. We see therefore that the problem simplifies
to the fixed strike Asian pricing problem.
Let us write down the stochastic differential equation that rules the process (ξt)t∈[0,T ]. Using Itô’s lemma,
we get {
dξt =
ξ0−ξt
dt+ ξt
σdWt + (γ +
ξ0 = S0.
Note that we are faced with a singularity problem near 0 because of the term ξ0−ξt
. We are going to reduce
its effect using another change of variables.
Using Itô’s lemma, we show that
C0 = E
e−rTf
where Xt = log(ξt/ξ0) solves the following SDE
dXt = σdWt + γdt+
e−Xt−1
X0 = 0.
Lemma 8 — Existence and strong uniqueness hold for the stochastic differential equation (16).
Proof : Existence is obvious since we have a particular solution Xt. The diffusion coefficient being
constant and the drift coefficient being a decreasing function in the spatial variable, we have also strong
uniqueness for the SDE (see the proof of Proposition 4). ✷
Because of the singularity of the term e
−Xt−1
in the drift coefficient, the law of (Xt)t≥0 is not absolutely
continuous with respect to the law of (σWt)t≥0. That is why we now define (Zt)t≥0 by the following SDE
with an affine inhomogeneous drift coefficient :
dZt = σdWt + γdt−
Z0 = X0 = 0.
The drift coefficient exhibits the same behavior as the one in (16) in the limit t → 0 in order to ensure the
desired absolute continuity property. It is affine in the spatial variable so that (Zt)t≥0 is a Gaussian process
and as such is easy to simulate recursively.
Lemma 9 — The process
sdWs +
t (18)
is the unique solution of the stochastic differential equation (17).
Proof : Using Itô’s Lemma, we easily check that Zt given by (18) is a solution of (17). Again, constant
diffusion coefficient and decreasing drift coefficient ensures strong uniqueness. ✷
Remark 10 — For the computation of the price C0 = E
e−rT (S0e
XT −K)+
of a standard Asian call
option, the random variable e−rT (S0e
ZT − K)+ provides a natural control variate. Indeed, since ZT is a
Gaussian random variable with mean γ
T and variance σ
, one has
e−rT (S0e
ZT −K)+
= S0e
−r)TN
−Ke−rTN (d)
where N is the cumulative standard normal distribution function and d = log(S0/K)+
Notice that in Kemna and Vorst [13], the authors suggest the use of the control variate
S0 exp
σWt + γt dt
which has the same law than e−rT
ZT −K
σWt + γt dt is also a Gaussian variable with mean
T and variance σ
In order to define a new probability measure under which (Zt)t≥0 solves the SDE (16), one introduces
Lt = exp
e−Zs − 1 + Zs
dWs −
e−Zs − 1 + Zs
Because of the singularity of the coefficients in the neighborhood of s = 0, one has to check that the integrals
in Lt are well defined. This relies on the following lemma
Lemma 11 — Let ǫ > 0. In a random neighborhood of s = 0, we have
|Zs| ≤ cs
−ǫ and |Xs| ≤ cs
where c is a constant depending on σ,γ and ǫ.
Since ∀ǫ > 0,
∀z ≤ cs 12−ǫ,
e−z − 1 + z
≤ Cs−4ǫ,
we can choose ǫ < 1
to deduce that Lt is well defined.
Proof : We easily check that the Gaussian process (Bt)t∈[0,T ] defined by Bt =
∫ (3t) 13
sdWs is a standard
Brownian motion. Thanks to the law of iterated logarithm for the Brownian motion (see for example Karatzas
and Shreve [12] p. 112), there exists t1(ω) such that
∀t ≤ t1(ω), |Bt(ω)| ≤ t
4ω is an element of the underlying probability space Ω.
Therefore,
∀t ≤ (3t1(ω))
3 , |Zt(ω)| =
(ω) +
∣∣ ≤ σ
Taking c = max( σ
) yields
∀t ≤ (3t1(ω))
3 ∧ 1, |Zt(ω)| ≤ ct
On the other hand, recall that Xt = log(ξt/ξ0) = log
eσWt+γt
e−σWu−γudu
. So, using the law of
iterated logarithm for the Brownian motion, we deduce that there exists t2(ω) such that
∀t ≤ t2(ω), 0 ≤
eσWt(ω)+γt
e−σWu(ω)−γudu ≤ 1
−ǫ−γudu.
Denote g(t) = 1
−ǫ−γudu and let us investigate the order in time near zero of this
function. We have that
+γt = 1 + σt
−ǫ +O(t1−2ǫ)∫ t
−ǫ−γudu = t+
−ǫ +O(t2−2ǫ)
hence
g(t) = 1 + (σ +
−ǫ +O(t1−2ǫ),
so Xt(ω) ≤ log (g(t)) ∼
−ǫ, which ends the proof for Xt. ✷
Proposition 12 — (Lt)t∈[0,T ] is a martingale and, consequently, for all g : C([0, T ]) → R measurable, the
random variables g((Xt)0≤t≤T ) and g((Zt)0≤t≤T )LT are simultaneously integrable and then
g((Xt)0≤t≤T )
g((Zt)0≤t≤T )LT
Proof : The proof is similar to the proof of Proposition 4.
We have already shown existence and strong uniqueness for both SDE (16) and (17). Showing that the
stopping time
τn(Y ) = inf
t ∈ R+ such that
e−Ys − 1 + Ys
ds ≥ n
, with the convention inf{∅} = +∞,
have infinite limits when n tends to +∞, QX and QZ almost surely, follows from the previous lemma.
One has
LT = exp
e−Zt − 1 + Zt
dZt −
e−Zt − 1 + Zt
e−Zt − 1 + Zt
+ γ − Zt
Set A(t, z) =
1− z + z2
− e−z
. The function A : ]0, T ]×R→ R is continuously differentiable in time and
twice continuously differentiable in space. So, we can apply Itô’s Lemma on the interval [ǫ, T ] for ǫ > 0 :
A(T, ZT ) = A(ǫ, Zǫ) +
e−Zt − 1 + Zt
dZt −
1− Zt + Z
− e−Zt
1− e−Zt
Using the lemma 9, we let ǫ→ 0 to obtain
A(T, ZT ) =
e−Zt − 1 + Zt
dZt −
1− Zt + Z
− e−Zt
1− e−Zt
LT = exp
A(T, ZT )−
φ(t, Zt)dt
where φ is the mapping
φ(t, z) =
e−z − 1 + z − z
1− e−z
e−z − 1 + z
e−z − 1 + z
+ γ − z
. (19)
By (15) and Proposition 12, we get
C0 = E
e−rTf(S0e
ZT ) exp
A(T, ZT )−
φ(t, Zt)dt
. (20)
Since for each t > 0, lim
φ(t, z) = +∞ and lim
φ(t, z) = −∞, it is not possible to apply the exact
algorithm. One can use the unbiased estimator, at least theoretically, if there exists a random variable cZ
measurable with respect to Z such that
eA(T,ZT )−(r+cZ)T |f(S0eZT )|e
|cZ−φ(t,Zt)|dt
Unfortunately, this reinforced integrability condition is never satisfied :
Lemma 13 — Assume that f is a non identically zero function. Let pZ and qZ denote respectively a positive
probability measure on N and a positive probability density on [0, T ]. Let N be distributed according to pZ
and (Ui)i∈N∗ be a sequence of independent random variables identically distributed according to the density
qZ , both independent conditionally on the process (Zt)t∈[0,T ]. Then the random variable
eA(T,ZT )−rT f(S0e
pZ(N)N !
−φ(Ui, ZUi)
qZ(Ui)
is non integrable.
Proof : By conditioning on Z, one has
∆ := E
eA(T,ZT )−rT |f(S0eZT )|
pZ(N)N !
|φ(Ui,ZUi )|
qZ (Ui)
eA(T,ZT )−rT |f(S0eZT )|e
|φ(t,Zt)|dt
eA(T,ZT )−rT |f(S0eZT )|e
|φ(t,Zt)|dt
One can easily show that, ∀z < 0 and ∀t ∈ [T
, T ], φ(t, z) ≥ φ(z) where
φ(z) =
e−z − 1 + z − z2
e−z − 1 + z
e−z − 1 + z
+ γ+ − 2 z
Since φ(z) ∼
σ2T 2
, there exists c < 0 such that for all z < c, φ(z) ≥ e−2z
σ2T 2
. Hence,
∆ ≥ E
eA(T,ZT )−rT |f(S0eZT )|e
e−2Zt1{Zt<c}dt
eA(T,ZT )−rT |f(S0eZT )|e−
2σ2T e
e−2Ztdt
Using Jensen’s inequality we get
∆ ≥ E
eA(T,ZT )−rT |f(S0eZT )|e−
2σ2T exp
We have seen in the proof of lemma 11 that Zt =
t where (Bt)t≥0 is a standard Brownian motion.
So, conditionally on ZT ,
Ztdt is a gaussian random variable and hence ∆ = +∞.
We are in a situation where eA(T,ZT )−rT |f(S0eZT )|E
[∣∣∣ 1pZ(N)N !
−φ(Ui,ZUi )
qZ (Ui)
∣∣∣(Zt)t∈[0,T ]
is non inte-
grable while eA(T,ZT )−rT |f(S0eZT )|
pZ (N)N !
−φ(Ui,ZUi )
qZ(Ui)
∣∣∣(Zt)t∈[0,T ]
]∣∣∣ is integrable since
e−rT |f(S0eZT )| exp
A(T, ZT )−
φ(t, Zt)dt
<∞. Then, a natural idea would consist in considering,
for a given n ∈ N∗, the random variable
eA(T,ZT )−rT |f(S0eZT )|E
∣∣∣∣∣∣
pZ(Nj)Nj !
−φ(U ji , ZUj
∣∣∣∣∣∣
∣∣∣(Zt)t∈[0,T ]
where (Nj)1≤j≤n are independent variables having the same law as N and
i )i∈N∗
1≤j≤n
are independent
sequences having the same law as (Ui)i∈N∗ , both independent conditionally on the process (Zt)t∈[0,T ]. The
following general result tells us that this is not sufficient to circumvent integrability problems.
Lemma 14 — Let Y and Z be two real random variables and g : R → R a given measurable function.
Assume that g(Z)E (Y |Z) is integrable while g(Z)E (|Y | |Z) is non integrable. Then, when (Yi)1≤i≤n is
a sequence of independent random variables having the same law as Y , ∀n ∈ N∗, the random variable
g(Z)E
i=1 Yi| |Z
is non integrable.
Proof : Denote by e, e1 and en three functions satisfying
∀z ∈ R, e(z) = E (Y |Z = z) , e1(z) = E (|Y1| |Z = z) and en(z) = E
(∣∣∣∣∣
∣∣∣∣∣
|Z = z
On the one hand, since
|g(z)| |e(z)|PZ(dz) <∞ and
|g(z)| e1(z)PZ(dz) = +∞ , where PZ is the law of
Z, we have that
|g(z)| e1(z)1{e1(z)≥2|e(z)|}PZ(dz) = +∞.
On the other hand, ∀z ∈ R,
en(z) ≥
(∣∣∣∣∣
∣∣∣∣∣
1{∀2≤j≤n,Yj≥0}|Z = z
(∣∣∣∣∣
∣∣∣∣∣
1{∀2≤j≤n,Yj<0}|Z = z
Y +1 |Z = z
P (Y1 ≥ 0|Z = z)n−1 + E
Y −1 |Z = z
P (Y1 < 0|Z = z)n−1
e1(z)+e(z)
P (Y1 ≥ 0|Z = z)n−1 + e1(z)−e(z)2 P (Y1 < 0|Z = z)
e1(z)
1{e1(z)≥2|e(z)|}P (Y1 ≥ 0|Z = z)
e1(z)
1{e1(z)≥2|e(z)|}P (Y1 < 0|Z = z)
≥ e1(z)
1{e1(z)≥2|e(z)|}
Hence, E
g(Z)E
(∣∣ 1
i=1 Yi
∣∣ |Z
|g(z)|en(z)PZ(dz) = +∞. ✷
There is still hope yet. In the proof of Lemma 13, we saw that integrability problems appear when
Zt takes large negative values so that φ(t, Zt) tends rapidly towards +∞. Since lim
φ(t, z) = −∞, one
possible issue is to split the function φ(t, Zt) into a positive part and a negative part. The first term can be
handled by the exact simulation technique whereas the second term, which as we shall see in the following
section presents no integrability problems, can be handled by the unbiased estimator technique.
2.2.1 An hybrid pseudo-exact method
We rewrite (20) in the following form
C0 = E
eA(T,ZT )−rTf(S0e
ZT )e
φ−(t,Zt)dte−
φ+(t,Zt)dt
. (22)
Let pZ and qZ denote respectively a positive probability measure on N and a positive probability density on
[0, T ]. Let N be distributed according to pZ and (Ui)i∈N∗ be a sequence of independent random variables
identically distributed according to the density qZ , both independent conditionally on the process (Zt)t∈[0,T ].
Note that, since eA(T,ZT )−rTf(S0e
ZT )e
|φ−(t,Zt)|dte−
φ+(t,Zt)dt = eA(T,ZT )−rT f(S0e
ZT )e−
φ(t,Zt)dt is
integrable, one has
C0 = E
eA(T,ZT )−rT f(S0e
pZ(N)N !
φ−(Ui, ZUi)
qZ(Ui)
φ+(t,Zt)dt
. (23)
Remark 15 — There is no hope that this estimator is square integrable. Indeed, one can show as in Lemma
13 that E
−(t,Zt))
= +∞ since (φ−(t, z))2 is of order z4 for large positive z.
The idea then is to apply the exact simulation technique to simulate an event with probability e−
φ+(t,Zt)dt.
Since for each t > 0, lim
φ+(t, z) = +∞, one needs to bound from above φ+(t, z), uniformly with respect
to t ∈ [0, T ], for z > c where c < 0 is a given constant. Thanks to the following lemma, it is possible to do
so but only uniformly with respect to t ∈ [ǫ, T ] for all ǫ > 0 :
Lemma 16 — For all 0 < t ≤ T ,
φ+(t, z) ≤ γ
∀c < 0, sup
z∈[c,0]
φ+(t, z) ≤ e
−c − 1 + c
(1 + γ+t) +
(e−c − 1)2
2σ2t2
Proof : Let z > 0. It is useful to distinguish two cases according to the sign of γ :
1. γ ≥ 0
We rewrite φ in the following form
φ(t, z) =
e−z − 1 + z − z
1− e−z
2 − (z ∧ 1)2
2σ2t2
(e−z − 1)2 − (z ∧ 1)2
2σ2t2
First note that
e−z−1+z− z2
≤ 0, 1−e
(e−z−1)2−(z∧1)2
2σ2t2
≤ 0. Moreover,
2 − (z ∧ 1)2
2σ2t2
if γt ≤ 1
otherwise
Consequently, φ+(t, z) ≤ γ
2. γ ≤ 0
Now we rewrite φ in the following form
φ(t, z) =
e−z − 1 + z − z
e−z − 1 + z
(e−z − 1)2 − z2
2σ2t2
1− e−z
It is then easy to show that φ+(t, z) ≤ 1
Note that 1
. Hence, gathering the two cases yields the first part of the lemma.
Let now z ∈ [c, 0] for a given negative constant c. We rewrite φ in the following form
φ(t, z) =
e−z − 1 + z
(1 + γ+t) +
(e−z − 1)2
2σ2t2
σ2t2︸ ︷︷ ︸
≥0 for z<0
1− e−z
− γ− e
−z − 1 + z
σ2t︸ ︷︷ ︸
≤0 for z<0
Since ∂z
e−z − 1 + z
(1 + γ+t) +
(e−z − 1)2
2σ2t2
1− e−2z − 2z + tγ+(1− e−z)
is negative for all
z < 0, one has that
z∈[c,0]
φ+(t, z) ≤ e
−c − 1 + c
(1 + γ+t) +
(e−c − 1)2
2σ2t2
This lemma suggests to apply the exact algorithm on [ǫ, T ] for a fixed positive threshold ǫ. It remains
to handle the time interval [0, ǫ[. Thanks to the following lemma, we that φ+(t, Zt) can be approximately
bounded from above for small t, almost surely, by a function of t. The idea is then to extend the exact
simulation algorithm by simulating an inhomogeneous Poisson process. Of course, this hybrid method is no
longer exact since the positive threshold for which the upper bound holds is random.
Lemma 17 — For all η > 0, there exists a random neighborhood of t = 0 such that
φ+(t, Zt) ≤
−η (24)
where c = max( σ
Proof : We rewrite (19) this way
φ(t, z) =
1− e−z
e−z − 1 + z
1− z + z2
− e−z − 1
(e−z − 1 + z)(e−z − 1− z)
and make the following Taylor expansions
1− z + z
− e−z − 1
(e−z − 1 + z)(e−z − 1− z)
z3 +O(z4)
1− e−z
e−z − 1 + z
z +O(z2).
On the other hand, we have seen in the proof of lemma 11 that there exists a random neighborhood of zero
such that Zt ≤ ct
−η where c = max( σ
). We conclude that, in a random neighborhood of zero,
φ+(t, Zt) ≤
2.2.2 Numerical computation
For numerical computation, we are going to use the following set of parameters : S0 = 100, K = 100,
σ = 0.2, r = 0.1, δ = 0 and T = 1. To fix the ideas, let us consider a call option. The price C0 writes as
follows
C0 = E
eA(T,ZT )−rT
ZT −K
φ−(Ui, ZUi)
φ+(t,Zt)dt
where N ∼ P(cp) and (Ui)i≥1 is an independent sequence of independent random variables uniformly dis-
tributed in [0, T ]. The parameter cp > 0 is set to one in the following. We give a description of the hybrid
method we implement :
Algorithm 2
On the time interval Ij := [
1. Simulate Z T
, Z T
and a lower bound mj for the minimum of (Zt)t∈Ij (use the fact that Zt =
t where (Bt)t≥0 is a standard Brownian motion).
2. Find M j > 0 such that ∀t ∈ Ij , φ+(t, Zt) ≤M j (use Lemma 16).
3. Simulate an homogeneous spatial Poisson process on the rectangle Ij × [0,M j] and accept (respectively
reject) the trajectory simulated if the number of points falling below the graph (φ+(t, Zt))t∈Ij is equal
to (respectively different from) zero.
Carry on this acceptance rejection algorithm until reaching a time interval IJ for a chosen J ∈ N∗. On the
remaining time interval [0, T
], use the same acceptance/rejection algorithm but with an inhomogeneous
spatial Poisson process this time (use Lemma 17).
In table 3, we give the price obtained by our method for different values of the positive threshold ǫ = T
The number M of Monte Carlo simulations is equal to 105 and the true price is equal to 7.042 (computed
using a Monte Carlo method with a trapezoidal scheme and a Kemna-Vorst control variate technique).
Price CPU
ǫ = T
6.9394 7s
ǫ = T
6.9590 10s
ǫ = T
6.9703 13s
ǫ = T
6.9952 17s
ǫ = T
7.0423 21s
Table 3: Price of the Asian call using the hybrid-pseudo exact method.
Clearly, the method is not yet competitive regarding computation time. Nevertheless, unlike the usual
discretization methods, it is not prone to discretization errors.
3 Conclusion
In this article, we have applied two original Monte Carlo methods for pricing Asian like options which
have the following pay-off : (αST + β
Stdt −K)+. In the case α 6= 0, we applied both the algorithm of
Beskos et al. [1] and a method based on the unbiased estimator of Wagner [23] and more recently the Poisson
estimator of Beskos et al. [2] and the generalized Poisson estimator of Fearnhead et al. [6]. The numerical
results show that the latter performs the best. The more interesting case α = 0, which corresponds to usual
continuously monitored Asian options, can not be treated using neither the exact algorithm, nor the method
of exact computation of expectation but we investigate an hybrid pseudo-exact method which combines the
two techniques. More generally, this hybrid method is an extension of the two exact methods and can be
applied in other situations.
From a practical point of view, the main contribution of these techniques is to allow Monte Carlo pricing
without resorting to discretization schemes. Hence, we are no longer prone to the discretization bias that
we encounter in standard Monte Carlo methods for pricing Asian like options. Even though these exact
methods are time consuming, they provide a good and reliable benchmark.
References
[1] A. Beskos, O. Papaspiliopoulos, and Gareth O. Roberts. Retrospective exact simulation of diffusion
sample paths. Bernoulli, 12(6), December 2006.
[2] A. Beskos, O. Papaspiliopoulos, Gareth O. Roberts, and Paul Fearnhead. Exact and computationally
efficient likelihood-based estimation for discretely observed diffusion processes. To appear in the Journal
of the Royal Statistical Society, Series B.
[3] M. Broadie and P. Glasserman. Estimating security price derivatives using simulation. Management
Science, 42(2):269–285, 1996.
[4] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E Knuth. On the lambert W function.
Advances in Computational Mathematics, 5:329–359, 1996.
[5] F. Dubois and T. Lelievre. Efficient pricing of asian options by the pde approach. Journal of Compu-
tational Finance, 8(2), 2004.
[6] Paul Fearnhead, O. Papaspiliopoulos, and Gareth O. Roberts. Particle filters for partially observed
diffusions. Working paper. Lancaster University., 2006.
[7] M. Fu, D. Madan, and T. Wang. Pricing continuous asian options : a comparison of monte carlo and
laplace transform inversion methods. Journal of Computational Finance, 2(2), 1999.
[8] H. Geman and A. Eydeland. Domino effect. Risk, pages 65–67, April 1995.
[9] H. Geman, N. El Karoui, and J.C. Rochet. Changes of numéraires, changes of probability measure and
option pricing. J. Appl. Probab., 32(2):443–458, 1995.
[10] H. Geman and M. Yor. Bessel processes, asian option and perpetuities. Mathematical Finance, 3(4),
1993.
[11] J.E. Ingersoll. Theory of Financial Decision Making. Rowman & Littlefield, 1987.
[12] I. Karatzas and Steven E. Shreve. Brownian motion and stochastic calculus. Springer-Verlag New-York,
second edition, 1991.
[13] A. Kemna and A. Vorst. A pricing method for options based on average asset values. Journal of Banking
and Finance, 14(1):113–129, 1990.
[14] B. Lapeyre and E. Temam. Competitive Monte Carlo methods for pricing asian options. Journal of
Computational Finance, 5(1), 2001.
[15] E. Levy. Pricing european average rate currency options. Journal of International Money and Finance,
11(5):474–491, October 1992.
[16] R. Lord. Partially exact and bounded approximations for arithmetic Asian options. Journal of Com-
putational Finance, 10(2), 2006.
[17] D. Revuz and M. Yor. Continuous martingales and Brownian motion. Springer-Verlag Berlin Heidelberg,
1991.
[18] L. C. G. Rogers and Z. Shi. The value of an Asian option. J. Appl. Probab., 32(4):1077–1088, 1995.
[19] T. H. Rydberg. A note on the existence of unique equivalent martingale measures in a markovian
setting. Finance and Stochastics, 1(3):251–257, 1997.
[20] S. Turnball and L. Wakeman. A quick algorithm for pricing european average options. Journal of
Financial and Quantitative Analysis, 16:377–389, 1991.
[21] J. Vecer. A new pde approach for pricing arithmetic asian options. Journal of Computational Finance,
4(4), 2001.
[22] T. Vorst. Prices and hedge ratios of average exchange rate options. International Review of Financial
Analysis, 1(3):179–193, 1992.
[23] W. Wagner. Unbiased Monte Carlo evaluation of certain functional integrals. J. Comput. Phys.,
71(1):21–33, 1987.
[24] W. Wagner. Monte Carlo evaluation of functionals of solutions of stochastic differential equations.
Variance reduction and numerical examples. Stochastic Anal. Appl., 6(4):447–468, 1988.
[25] W. Wagner. Unbiased multi-step estimators for the Monte Carlo evaluation of certain functional inte-
grals. J. Comput. Phys., 79(2):336–352, 1988.
[26] W. Wagner. Unbiased Monte Carlo estimators for functionals of weak solutions of stochastic differential
equations. Stochastics Stochastics Rep., 28(1):1–20, 1989.
[27] B. Wong and C. C. Heyde. On the martingale property of stochastic exponentials. J. Appl. Probab.,
41(3):654–664, 2004.
4 Appendix
4.1 The practical choice of p and q in the U.E method
The best choice for the probability law p of N and the common density q of the variables (Vi)i≥1
is obviously the one for which the variance of the simulation is minimum. In a very general setting, it is
difficult to tackle this issue. In order to have a first idea, we are going to restrict ourselves to the computation
p(N)N !
g(Vi)
q(Vi)
where g : [0, T ] → R.
Lemma 18 — When g is a measurable function on [0, T ] such that 0 <
|g(t)|dt < +∞, the variance of
p(N)N !
g(Vi)
q(Vi)
is minimal for
qopt(t) =
|g(t)|
|g(t)|dt
1[0,T ](t) and popt(n) =
|g(t)|dt
|g(t)|dt
Proof : Minimizing the variance in (7) comes down to minimizing the expectation of the square of
p(N)N !
g(Vi)
q(Vi)
F (p, q) = E
(p(N)N !)2
g2(Vi)
q2(Vi)
g2(t)
p(n) (n!)2
Using Cauchy-Schwartz inequality we obtain a lower bound for F (p, q)
F (p, q) =
g2(t)
p(n)n!
p(n) ≥
g2(t)
q(t)dt
|g(t)|dt
= exp
|g(t)|dt
We easily check that this lower bound is attained for qopt and popt.
The optimal probability distribution popt is the Poisson law with parameter
|g(t)|dt. This justifies our
use of a Poisson distribution for p.
4.2 Simulation from the distribution h given by (13)
Recall that
h(u) = C exp
A(u)− (u−X0)
= C exp
(1− e−σu)− (u−X0)
where C is a normalizing constant.
The expansion of the exponential e−σu at the first order yields
h(u) ≈ C exp
u− (u−X0)
= C exp
(u− (X0 + T (γ+βS0)σ ))
This suggests to do rejection sampling using the normal distribution with mean X0 +
T (γ+βS0)
variance T as prior. Unfortunately, for a standard set of parameters, this method gives bad results. Even a
second order expansion of e−σu which also modifies the variance does not work.
In order to get round this problem, we evaluate the mode u∗ of h. We have
h′(u∗) = C
∗ −X0
(1− e−σu
)− (u
∗ −X0)2
So, h′(u∗) = 0 if and only if
∗ −X0
which writes
σ(u∗ −X0 −
T )eσ(u
∗−X0− γσT ) = TβS0e
−σX0−γT .
The function x 7→ xex is continuous and increasing on [0,+∞[ and so is its inverse which we denote by
W . Since TβS0e
−σX0−γT ≥ 0, we deduce that h is unimodal and that its mode satisfies
γT +W
βS0Te
−γT−σX0
+ σX0
The functionW is the well-known Lambert function, also called the Omega function. It is uniquely valued
on [0,+∞[ and there are robust and fast numerical methods based on series expansion for approximating
this function (see for example Corless et al. [4]).
Numerical tests showed that performing rejection sampling using a Gaussian distribution with variance
T and mean u∗ instead of X0 +
T (γ+βS0)
gives plain satisfaction. In table 4, we see that for arbitrary choice
of the parameter α
, the acceptance rate of the algorithm is always high (of order 70%) and that the
computation time is low.
Nb of simulations Acceptance rate Computation time
0.2 61% 3s
0.5 106 68% 3s
0.8 80% 2s
Table 4: Acceptance rate of the rejection algorithm of simulating from the distribution h in (13) with
S0 = 100, σ = 0.3, T = 2 and r = 0.1.
	1 Exact Simulation techniques
	1.1 The exact simulation method of Beskos et al. BPR
	1.2 The unbiased estimator (U.E)
	2 Application : the pricing of continuous Asian options
	2.1 The case =0
	2.1.1 Numerical computation
	2.2 Standard Asian options : the case = 0 and > 0
	2.2.1 An hybrid pseudo-exact method
	2.2.2 Numerical computation
	3 Conclusion
	Bibliography
	4 Appendix
	4.1 The practical choice of p and q in the U.E method 
	4.2 Simulation from the distribution h given by (13)

  Taking advantage of the recent litterature on exact simulation algorithms
(Beskos, Papaspiliopoulos and Roberts) and unbiased estimation of the
expectation of certain fonctional integrals (Wagner, Beskos et al. and
Fearnhead et al.), we apply an exact simulation based technique for pricing
continuous arithmetic average Asian options in the Black and Scholes framework.
Unlike existing Monte Carlo methods, we are no longer prone to the
discretization bias resulting from the approximation of continuous time
processes through discrete sampling. Numerical results of simulation studies
are presented and variance reduction problems are considered.

Introduction
Although the Black & Scholes framework is very simple, it is still a challenging task to efficiently price
Asian options. Since we do not know explicitly the distribution of the arithmetic sum of log-normal variables,
there is no closed form solution for the price of an Asian option. By the early nineties, many researchers
attempted to address this problem and hence different approaches were studied including analytic approxima-
tions (see Turnball and Wakeman [20], Vorst [22], Levy [15] and more recently Lord [16]), PDE methods (see
Vecer [21], Rogers and Shi [18], Ingersoll [11], Dubois and Lelievre [5]), Laplace transform inversion methods
(see Geman and Yor [10], Geman and Eydeland [8]) and, of course, Monte Carlo simulation methods (see
Kemna and Vorst [13], Broadie and Glasserman [3], Fu et al. [7]).
Monte Carlo simulation can be computationally expensive because of the usual statistical error. Variance
reduction techniques are then essential to accelerate the convergence (one of the most efficient techniques is
the Kemna&Vorst control variate based on the geometric average). One must also account for the inherent
discretization bias resulting from approximating the continuous average of the stock price with a discrete
one. It is crucial to choose with care the discretization scheme in order to have an accurate solution (see
Lapeyre and Temam [14]). The main contribution of our work is to fully address this last feature by the use,
after a suitable change of variables, of an exact simulation method inspired from the recent work of Beskos
et al. [1, 2] and Fearnhead et al. [6].
In the first part of the paper, we recall the algorithm introduced by Beskos et al. [1] in order to simulate
sample-paths of processes solving one-dimensional stochastic differential equations. By a suitable change of
variables, one may suppose that the diffusion coefficient is equal to one. Then, according to the Girsanov
theorem, one may deal with the drift coefficient by introducing an exponential martingale weight. Because
of the one-dimensional setting, the stochastic integral in this exponential weight is equal to a standard
1Project team Math Fi, CERMICS, Ecole des Ponts, Paristech, supported by the ANR program ADAP’MC.
Postal address : 6-8 av. Blaise Pascal, Cité Descartes, Champs-sur-Marne, 77455 Marne-la-Vallée Cedex 2.
E-mails : jourdain@cermics.enpc.fr and sbai@cermics.enpc.fr
http://arxiv.org/abs/0704.1433v3
integral with respect to the time variable up to the addition of a function of the terminal value of the
path. Under suitable assumptions, conditionally on a Brownian path, an event with probability equal to
the normalized exponential weight can be simulated using a Poisson point process. This allows to accept or
reject this Brownian path as a path solution to the SDE with diffusion coefficient equal to one. In finance,
one is interested in computing expectations rather than exact simulation of the paths. In this perspective,
computation of the exponential importance sampling weight is enough. The entire series expansion of the
exponential function permits to replace this exponential weight by a computable weight with the same
conditional expectation given the Brownian path. This idea was first introduced by Wagner [23, 24, 25, 26]
in a statistical physics context and it was very recently revisited by Beskos et al. [2] and Fearnhead et al.
[6] for the estimation of partially observed diffusions. Some of the assumptions necessary to implement the
exact algorithm of Beskos et al. [1] can then be weakened.
The second part is devoted to the application of these methods to option pricing within the
Black & Scholes framework. Throughout the paper, St = S0 exp
σWt + (r − δ −
represents the
stock price at time t, T the maturity of the option, r the short interest rate, σ the volatility parameter, δ
the dividend rate and (W )t∈[0,T ] denotes a standard Brownian motion on the risk-neutral probability space
(Ω,F ,P). We are interested in computing the price C0 = E
e−rTf
αST + β
of a European
option with pay-off f
αST + β
assumed to be square integrable under the risk neutral measure P.
The constants α and β are two given non-negative parameters.
When α > 0, we remark that, by a change of variables inspired by Rogers and Shi [18], αST + β
has the same law as the solution at time T of a well-chosen one-dimensional stochastic differential equation.
Then it is easy to implement the exact methods previously presented. The case α = 0 of standard Asian
options is more intricate. The previous approach does not work and we propose a new change of variables
which is singular at initial time. It is not possible to implement neither the exact simulation algorithm nor
the method based on the unbiased estimator of Wagner [23] and we propose a pseudo-exact hybrid method
which appears as an extension of the exact simulation algorithm. In both cases, one first replaces the integral
with respect to the time variable in the function f by an integral with respect to time in the exponential
function. Because of the nice properties of this last function, exact computation is possible.
1 Exact Simulation techniques
1.1 The exact simulation method of Beskos et al. [1]
In a recent paper, Beskos et al. [1] proposed an algorithm which allows to simulate exactly the solution
of a 1-dimensional stochastic differential equation. Under some hypotheses, they manage to implement an
acceptance-rejection algorithm over the whole path of the solution, based on recursive simulation of a biased
Brownian motion. Let us briefly recall their methodology. We refer to [1] for the demonstrations and a
detailed presentation.
Consider the stochastic process (ξt)0≤t≤T determined as the solution of a general stochastic differential
equation of the form : {
dξt = b(ξt)dt+ σ(ξt)dWt
ξ0 = ξ ∈ R
where b and σ are scalar functions satisfying the usual Lipschitz and growth conditions with σ non vanishing.
To simplify this equation, Beskos et al. [1] suggest to use the following change of variables : Xt = η(ξt)
where η is a primitive of 1
(η(x) =
Under the additional assumption that 1
is continuously differentiable, one can apply Itô’s lemma to get
dXt = η
′(ξt)dξt +
η′′(ξt) d< ξ, ξ >t
b(ξt)
σ(ξt)
dt+ dWt −
σ′(ξt)
b(η−1(Xt))
σ(η−1(Xt))
′(η−1(Xt))
︸ ︷︷ ︸
a(Xt)
dt+ dWt
So ξt = η
−1(Xt) where (Xt)t is a solution of the stochastic differential equation
dXt = a(Xt)dt+ dWt
X0 = x.
Thus, without loss of generality, one can start from equation (2) instead of (1).
Let us denote by (W xt )t∈[0,T ] the process (Wt+x)t∈[0,T ], by QWx its law and by QX the law of the process
(Xt)t∈[0,T ]. From now on, we will denote by (Yt)t∈[0,T ] the canonical process, that is the coordinate mapping
on the set C([0, T ],R) of real continuous maps on [0, T ] (see Revuz and Yor [17] or Karatzas and Shreve
[12]).
One needs the following assumption to be true
Assumption 1 : Under QWx , the process
Lt = exp
a(Yu)dYu −
a2(Yu)du
is a martingale.
According to Rydberg [19] (see the proof of Proposition 4 where we give his argument on a specific
example), a sufficient condition for this assumption to hold is
-Existence and uniqueness in law of a solution to the SDE (2).
-∀t ∈ [0, T ],
a2(Yu)du <∞, QX and QWx almost surely on C([0, T ],R).
Thanks to this assumption, one can apply the Girsanov theorem to get that QX is absolutely continuous
with respect to QWx and its Radon-Nikodym derivative is equal to
= exp
a(Yt)dYt −
a2(Yt)dt
Consider A the primitive of the drift a, and assume that
Assumption 2 : a is continuously differentiable.
Since, by Itô’s lemma, A(W xT ) = A(x) +
a(W xt )dW
a′(W xt )dt, we have
= exp
A(YT )−A(x)−
a2(Yt) + a
′(Yt)dt
Before setting up an acceptance-rejection algorithm using this Radon-Nikodym derivative, a last step is
needed. To ensure the existence of a density h(u) proportional to exp(A(u) − (u−x)
), it is necessary and
sufficient that the following assumption holds
Assumption 3 : The function u 7→ exp(A(u)− (u−x)
) is integrable.
Finally, let us define a process Zt distributed according to the following law QZ
(W xt )t∈[0,T ]|W xT = y
h(y)dy
where the notation L(.|.) stands for the conditional law. One has
= C exp
a2(Yt) + a
′(Yt)dt
where C is a normalizing constant. At this level, Beskos et al. [1] need another assumption
Assumption 4 : The function φ : x 7→ a
2(x)+a′(x)
is bounded from below.
Therefore, one can find a lower bound k of this function and eventually the Radon-Nikodym derivative
of the change of measure between X and Z takes the form
= Ce−kT exp
φ(Yt)− k dt
The idea behind the exact algorithm is the following : suppose that one is able to simulate a continuous
path Zt(ω) distributed according to QZ and let M(ω) be an upper bound of the mapping t 7→ φ(Zt(ω))− k.
Let N be an independent random variable which follows the Poisson distribution with parameter TM(ω) and
let (Ui, Vi)i=1...N be a sequence of independent random variables uniformly distributed on [0, T ]× [0,M(ω)].
Then, the number of points (Ui, Vi) which fall below the graph {(t, φ(Zt(ω))− k); t ∈ [0, T ]} is equal to zero
with probability exp
φ(Zt(ω))− k dt
. Actually, simulating the whole path (Zt)t∈[0,T ] is not necessary.
It is sufficient to determine an upper bound for φ(Zt)− k since, as pointed out by the authors, it is possible
to simulate recursively a Brownian motion on a bounded time interval by first simulating its endpoint, then
simulating its minimum or its maximum and finally simulating the other points2. For this reason, one needs
the following assumption for the algorithm to be feasible :
Assumption 5 : Either lim sup
φ(u) < +∞ or lim sup
φ(u) < +∞.
Suppose for example that lim sup
φ(u) < +∞. The exact algorithm of Bekos et al. [1] then takes the
following form :
Algorithm 1
1. Draw the ending point ZT of the process Z with respect to the density h.
2. Simulate the minimum m of the process Z given ZT .
3. Fix an upper bound M(m) = sup{φ(u)− k;u ≥ m} for the mapping t 7→ φ(Zt)− k.
4. Draw N according to the Poisson distribution with parameter TM(m) and draw (Ui, Vi)i=1...N , a
sequence of independent variables uniformly distributed on [0, T ]× [0,M(m)].
5. Fill in the path of Z at the remaining times (Ui)i=1...N .
2In their paper, the authors explain how to do such a decomposition of the Brownian path.
6. Evaluate the number of points (Vi)i=1...N such that Vi ≤ φ(ZUi)− k.
If it is equal to zero, then return the simulated path Z.
Else, return to step 1.
This algorithm gives exact skeletons of the process X , solution of the SDE (2). Once accepted, a path can
be further recursively simulated at additional times without any other acceptance/rejection criteria. We also
point out that the same technique can be generalized by replacing the Brownian motion in the law of the
proposal Z by any process that one is able to simulate recursively by first simulating its ending point, its
minimum/maximum and then the other points. Also, the extension of the algorithm to the inhomogeneous
case, where the drift coefficient a in (2), and therefore the function φ, depend on the time variable t, is
straightforward given that the assumptions presented above are appropriately modified.
1.2 The unbiased estimator (U.E)
In finance, the pricing of contingent claims often comes down to the problem of computing an expectation
of the form
C0 = E (f(XT )) (3)
where X is a solution of the SDE (2) and f is a scalar function such that f(XT ) is square integrable. In a
simulation based approach, one is usually unable to exhibit an explicit solution of this SDE and will therefore
resort to numerical discretization schemes, such as the Euler or Milstein schemes, which introduce a bias.
Of course, the exact algorithm presented above avoids this bias. Here, we are going to present a technique
which permits to compute exactly the expectation (3) while assumptions 4 and 5 on the function a
which
appears in the Radon-Nikodym derivative are relaxed.
Using the previous results and notations, we get, under the assumptions 1 and 2, that
C0 = E
f(W xT ) exp
A(W xT )−A(x) −
a2(W xt ) + a
′(W xt )dt
. (4)
In order to implement an importance sampling method, let us introduce a positive density ρ on the real
line and a process (Zt)t∈[0,T ] distributed according to the following law QZ
(W xt )t∈[0,T ]|W xT = y
ρ(y)dy.
By (4), one has
C0 = E
ψ(ZT ) exp
φ(Zt)dt
where ψ : z 7→ f(z) e
A(z)−A(x)−
(z−x)2
2πρ(z)
and φ : z 7→ a
2(z)+a′(z)
. We do not impose ρ to be equal to the
density h of the previous section. It is a free parameter chosen in such a way that it reduces the variance of
the simulation.
In his first paper, Wagner [23] constructs an unbiased estimator of the expectation (5) when ψ is a
constant, (Zt)t∈[0,T ] is an R
d−valued Markov process with known transition function and φ is a measurable
function such that E
|φ(Zt)|dt
< +∞. His main idea is to expand the exponential term in a power
series, then, using the transition function of the underlying Markov process and symmetry arguments, he
constructs a signed measure ν on the space Y =
n=0([0, T ]× Rd)n+1 such that the expectation at hand is
equal to ν(Y). Consequently, any probability measure µ on Y that is absolutely continuous with respect to
ν gives rise to an unbiased estimator ζ defined on (Y, µ) via ζ(y) = dν
(y). In practice, a suitable way to
construct such an estimator is to use a Markov chain with an absorbing state. Wagner also discusses variance
reduction techniques, specially importance sampling and a shift procedure consisting on adding a constant c
to the integrand φ and then multiplying by the factor e−cT in order to get the right expectation. Wagner [25]
extends the class of unbiased estimators by perturbing the integrand φ by a suitably chosen function φ0 and
then using mixed integration formulas representation. Very recently, Beskos et al. [2] obtained a simplified
unbiased estimator for (5), termed Poisson estimator, using Wagner’s idea of expanding the exponential in
a power series and his shift procedure. To be specific, the Poisson estimator writes
ψ(ZT )e
cpT−cT
c− φ(ZVi)
where N is a Poisson random variable with parameter cP and (Vi)i is a sequence of independent random
variables uniformly distributed on [0, T ]. Fearnhead et al. [6] generalized this estimator allowing c and cP to
depend on Z and N to be distributed according to any positive probability distribution on N. They termed
the new estimator the generalized Poisson estimator. We introduce a new degree of freedom by allowing the
sequence (Vi)i to be distributed according to any positive density on [0, T ]. This gives rise to the following
unbiased estimator for (5) :
Lemma 1 — Let pZ and qZ denote respectively a positive probability measure on N and a positive probability
density on [0, T ]. Let N be distributed according to pZ and (Vi)i∈N∗ be a sequence of independent random
variables identically distributed according to the density qZ , both independent from each other conditionally
on the process (Zt)t∈[0,T ]. Let cZ be a real number which may depend on Z. Assume that
|ψ(ZT )|e−cZT exp
|cZ − φ(Zt)|dt
ψ(ZT )e
−cZT 1
pZ(N)N !
cZ − φ(ZVi)
qZ(Vi)
is an unbiased estimator of C0.
Proof : The result follows from
ψ(ZT )e
−cZT 1
pZ(N)N !
cZ − φ(ZVi)
qZ(Vi)
∣∣∣(Zt)t∈[0,T ]
= ψ(ZT )e
cZ − φ(Zt)dt
pZ(n)n!
pZ(n)
= ψ(ZT ) exp
φ(Zt)dt
Using (7), one is now able to compute the expectation at hand by a simple Monte Carlo simulation. The
practical choice of pZ and qZ conditionally on Z is studied in the appendix 4.1.
As pointed out in Fearnhead et al. [6], this method is an extension of the exact algorithm method since,
under assumptions 3, 4 and 5, the reinforced integrability assumption of Lemma 1 is always satisfied.
Indeed, suppose for example that lim sup
φ(u) < +∞ and let k be a lower bound of φ, mZ be the
minimum of the process Z and MZ an upper bound of {φ(u)− k, u ≥ mZ}. Then, taking cZ = MZ + k in
Lemma 1 ensures the integrability condition :
|ψ(ZT )|e−(MZ+k)T e
|MZ+k−φ(Zt)|dt
|ψ(ZT )|e−(MZ+k)T e
MZ+k−φ(Zt)dt
|ψ(ZT )|e−
φ(Zt)dt
and hence, one is allowed to write that
C0 = E
ψ(ZT )e
−(MZ+k)T 1
pZ(N)N !
MZ + k − φ(ZVi)
qZ(Vi)
Better still, the random variable ψ(ZT )e
−(MZ+k)T 1
pZ(N)N !
MZ+k−φ(ZVi )
qZ(Vi)
is square integrable when
pZ is the Poisson distribution with parameterMZT + k and qZ is the uniform distribution on [0, T ] since we
have then
ψ(ZT )e
−(MZ+k)T 1
pZ(N)N !
MZ + k − φ(ZVi)
qZ(Vi)
 = E
ψ2(ZT )
1− φ(ZVi)
MZ + k
ψ2(ZT )
The last inequality follows from the square integrability of f : whenever one is able to simulate from the
density h, introduced in the exact algorithm, by doing rejection sampling, there exists a density ρ such that
ψ, which is equal to f(ZT )
h(ZT )
ρ(ZT )
up to a constant factor, is dominated by f and so is square integrable.
The square integrability property is very important in that we use a Monte Carlo method. We see that,
whenever the exact algorithm is feasible, the unbiased estimator of lemma 1 is a simulable square integrable
random variable, at least for the previous choice of pZ and qZ .
Remark 2 — One can derive two estimators of C0 from the result of Lemma 1 :
f(ZiT )
T )−A(x)−
2πρ(ZiT )
e−cZT
pZ(N i)N i!
cZ − φ(ZiV i
f(ZiT )
T )−A(x)−
2πρ(ZiT )
pZ(N i)N i!
cZ − φ(ZiV i
T )−A(x)−
2πρ(ZiT )
pZ(N i)N i!
cZ − φ(ZiV i
2 Application : the pricing of continuous Asian options
In the Black & Scholes model, the stock price is the solution of the following SDE under the risk-neutral
measure P
= (r − δ)dt+ σdWt (8)
where all the parameters are constant : r is the short interest rate, δ is the dividend rate and σ is the
volatility.
Throughout, we denote γ = r − δ − σ
. The path-wise unique solution of (8) is
St = S0 exp(σWt + γt) .
We consider an option with pay-off of the form
αST + β
where f is a given function such that E
αST + β
< ∞, T is the maturity of the option and
α, β are two given non negative parameters3. Note that for α = 0, this is the pay-off of a standard continuous
Asian option.
The fundamental theorem of arbitrage-free pricing ensures that the price of the option under consideration
C0 = E
e−rTf
αST + β
At first sight, the problem seems to involve two variables : the stock price and the integral of the stock
price with respect to time. Dealing with the PDE associated with Asian option pricing, Rogers and Rogers
and Shi [18] used a suitable change of variables to reduce the spatial dimension of the problem to one. We
are going to use a similar idea.
αS0 + βS0
e−σWu−γudu
eσWt+γt.
We have that
ξt = αS0e
σWt+γt + βS0
eσ(Wt−Wu)+γ(t−u)du
= αS0e
σBt+γt + βS0
eσBs+γsds
where we set Bs = Wt −Wt−s, ∀s ∈ [0, t]. Clearly, (Bs)s∈[0,t] is a Brownian motion and thus the following
lemma holds
Lemma 3 — ∀t ∈ [0, T ], ξt and αSt + β
Sudu have the same law.
As a consequence
C0 = E
e−rT f(ξT )
By applying Itô’s lemma, we verify that the process (ξt)t≥0 is a positive solution of the following 1-dimensional
stochastic differential equation for which path-wise uniqueness holds
dξt = βS0dt+ ξt(σdWt + (γ +
ξ0 = αS0.
We are thus able to value C0 by Monte Carlo simulation without resorting to discretization schemes using
one of the exact simulation techniques described in the previous section. In the case α = 0, one has to deal
with the fact that ξt starts from zero which is the reason why we distinguish two cases.
3The underlying of this option is a weighted average of the stock price at maturity and the running average of the stock
price until maturity with respective weights α and βT .
2.1 The case α 6= 0
We are going to apply both the exact algorithm of Beskos et al. [1] and the method based on the unbiased
estimator of lemma 1.
We make the following change of variables to have a diffusion coefficient equal to 1 :
log(ξt)
dXt = (
+ βS0
e−σXt)dt+ dWt
X0 = x with x =
log(αS0)
C0 = E
e−rTf(eσXT )
The following proposition ensures that assumption 1 is satisfied.
Proposition 4 — The process (Lt)t∈[0,T ] defined by
Lt = exp
e−σYt) dYt −
e−σYt)2dt
is a martingale under QWx .
Proof : Under QWx , (Lt)t∈[0,T ] is clearly a non-negative local martingale and hence a super-martingale.
Then, it is a true martingale if and only if E
(LT ) = 1.
Checking the classical Novikov’s or Kamazaki’s criteria is not straightforward. Instead, we are going to
use the approach developed by Rydberg [19] (see also Wong and Heyde [27]) who takes advantage of the link
between explosions of SDEs and the martingale property of stochastic exponentials.
Let us define the following stopping times :
τn(Y ) = inf
t ∈ R+ such that
e−σYu
du ≥ n
with the convention inf{∅} = +∞.
The stopped process (Lt∧τn(Y ))t∈[0,T ] is a true martingale under QWx since Novikov’s condition is fulfilled.
According to the Girsanov theorem, one can define a new probability measure QnX , which is absolutely
continuous with respect to QWx , by its Radon-Nikodym derivative
= LT∧τn(Y ).
Hence
1{τn(Y )>T}
1{τn(Y )>T}LT∧τn(Y )
Since (τn(Y ))n∈N is a non decreasing sequence, we can pass to the limit in the right hand side We get
X (τn(Y ) > T ) = EQWx
1{τ∞(Y )>T}LT∧τ∞(Y )
where τ∞(Y ) denotes the limit of the non decreasing sequence (τn(Y ))n∈N.
Under QWx , (Yt)t∈[0,T ] has the same law as a Brownian motion starting from x so τ∞(Y ) = +∞ ,QWx
almost surely, and consequently
= lim
X (τn(Y ) > T ) .
On the other hand, the Girsanov theorem implies that, under QnX , (Yt)t∈[0,T∧τn(Y )] solves a SDE of the
form (11). To conclude the proof, it is sufficient to check that trajectorial uniqueness holds for this SDE.
Indeed, the law of (Yt)t∈[0,T∧τn(Y )] under Q
X is the same as the law of (Yt)t∈[0,T∧τn(Y )] under QX . Hence
X (τn(Y ) > T ) = QX (τn(Y ) > T ) −→
QX (τ∞(Y ) > T ) .
Clearly,
+ βS0
e−σYu
du <∞, QX almost surely, so
= QX (τ∞(Y ) > T ) = 1
as required.
In order to check trajectorial uniqueness for the SDE (11), we consider two solutions X1 and X2. We
have that
d(X1t −X2t ) =
t − e−σX
dt ⇒ d|X1t −X2t | =
sign(X1t −X2t )
t − e−σX
|X1t −X2t | =
sign(X1s −X2s )
t − e−σX
ds ≤ 0.
The last inequality follows from the fact that x 7→ e−σx is a decreasing function. Finally, almost surely,
∀t ≥ 0, X1t = X2t which leads to strong uniqueness. ✷
Consequently, thanks to the Girsanov theorem, we have
= exp
e−σYt)
︸ ︷︷ ︸
a(Yt)
dYt −
e−σYt)2dt
 . (12)
Set A(u) =
a(x)dx = γ
u+ βS0
(1− e−σu). Then
= exp
A(YT )−A(x)−
a2(Yt) + a
′(Yt)dt
The function u 7→ exp
A(u)− (u−Y0)
= exp
u+ βS0
(1− e−σu)− (u−Y0)
is clearly integrable so we
can define a new process (Zt)t∈[0,T ] distributed according to the following law QZ
(Wt)t∈[0,T ]|WT = y
h(y)dy
where the probability density h is of the form
h(u) = C exp
A(u)− (u− Y0)
with C a normalizing constant. (13)
Remark 5 — Simulating from this probability distribution is not difficult (see the appendix 4.2 for an
appropriate method of acceptance/rejection sampling).
We have
= C exp
(a2(Yt) + a
′(Yt))dt
Set φ(x) =
a2(x)+a′(x)
e−σx)2−βS0e−σx
. A direct calculation gives
φ(x) =
if 2γ ≥ σ2
log( 2βS0
σ2−2γ )
otherwise.
Set k = infx∈R φ(x). Finally, we get
= Ce−kT exp
φ(Yt)− k dt
We check that
φ(x) =
φ(x) = +∞.
Hence we can apply the algorithm 1 to simulate exactly XT and compute C0 = E
e−rTf(eσXT )
Monte Carlo. On the other hand, using (12) we get
C0 = E
e−rTf(eσW
T ) exp
A(W xT )−A(x) −
a2(W xt ) + a
′(W xt )dt
and we can also use the unbiased estimator presented in the previous section to compute this expectation.
Remark 6 — We also applied the exact algorithm based on a geometric Brownian motion instead of the
standard Brownian motion which seems more intuitive given the form of the SDE (10). The algorithm is
feasible because we can simulate recursively a drifted Brownian motion and therefore a geometric Brownian
motion by an exponential change of variables. The results we obtained were not different from the first
method.
2.1.1 Numerical computation
For numerical tests, we consider the case
f(x) = (x−K)+
which corresponds to the European call option with strikeK. Using the exact simulation algorithm presented
above, we can simulate the underlying αST + β
Stdt at maturity (see Figure 1). Then, all we have to
do is a simple Monte Carlo method to get the price of the option under consideration. Using the unbiased
estimator, we get
C0 = E
e−rT (eσZT −K)+
eA(ZT )−A(x)−
(ZT −x)
2πρ(ZT )
e−(MZ+k)T
pZ(N)N !
MZ + k − φ(ZVi)
q(Vi)
where (Zt)t∈[0,T ], ρ,MZ , k, pZ and qZ are defined as in section 1.2. In order to ensure square integrability,
we choose pZ to be a Poisson distribution with parameterMZT +k and qZ to be the uniform distribution on
[0, T ]. For the density ρ, a good choice is to consider the density that we use to simulate from the distribution
h by rejection sampling.
We test these exact methods against a standard discretization scheme with the variance reduction tech-
nique of Kemna and Vorst [13]. As pointed out by Lapeyre and Temam [14], the discretization of the integral
by a simple Riemannian sum is not efficient. Instead, we use the trapezoidal discretization. In the sequel, we
will denote this method by Trap+KV. The table 1 gives the results we obtained for the following arbitrary
set of parameters : S0 = 100, K = 100, r = 0.05, σ = 0.3, δ = 0, T = 1, α = 0.6 and β = 0.4. The
computation has been made on a computer with a 2.8 Ghz Intel Penthium 4 processor. We intentionally
choose a large number of simulations in order to show the influence of the number of time steps when using
a discretization scheme.
10 53 97 140 183 227 270 313 357 400
0.000
0.002
0.004
0.006
0.008
0.010
0.012
0.014
0.016
0.018
PSfrag replacements
αST + β
Exact Simulation of the underlying : αST + β
Figure 1: Histogram of 105 independent realizations of αST + β
Stdt for α = 0.6 and β = 0.4 compared
with the lognormal distribution of ST .
Method M N Acceptance rate Price C.I at 95% CPU
10 11.46 [11.43, 11.48] 5 s
Trap+KV 20 106 - 11.46 [11.43, 11.49] 9 s
50 11.47 [11.44, 11.5] 21 s
Exact Simulation - 106 24% 11.46 [11.43, 11.5] 81 s
U.E (cP =MZ , cZ =MZ + k) - 10
6 - 11.46 [11.43, 11.49] 17 s
U.E (cP = cZ = 1/T ) - 10
6 - 11.46 [11.43, 11.49] 6 s
Table 1: Price of the option (9) using a standard discretization technique and exact simulation methods.
Empirical evidence shows that the exact simulation method is quite slow. This is mainly due to the
fact that the rejection algorithm has a little acceptance rate (24% according to table 1). Using a geometric
Brownian motion instead of a standard Brownian motion did not improve the results. Also, simulating
recursively a Brownian path conditionally on its terminal value and its minimum is time consuming.
The unbiased estimator is more efficient, especially when we can avoid the recursive simulation of the
Brownian path. To do so, we choose for pZ a Poisson distribution with mean cPT where cP is a free
parameter. If we assume that the integrability condition in lemma 1 holds, then we can write that
C0 = E
e−rT (eσZT −K)+
eA(ZT )−A(x)−
(ZT −x)
2πρ(ZT )
ecPT−cZT
cZ − φ(ZVi)
Regarding the dependence of the exact simulation method with respect to the parameters α and β, it
is intuitive that whenever α >> β, the method performs well since the logarithm of the underlying is not
far from the logarithm of the geometric Brownian motion on which we do rejection-sampling. The table 2
confirms this intuition. We see that we cannot apply the algorithm for small values of α and then let α→ 0
to treat the case α = 0.
α + β
0.3 0.4 0.5 0.6 0.7
Acceptance Rate 0.003% 0.47% 5.66% 24.43% 53.85%
Table 2: Influence of the parameter α
on the acceptance rate of the exact algorithm.
2.2 Standard Asian options : the case α = 0 and β > 0
A standard Asian option is a European option on the average of the stock price over a determined period
until maturity. An Asian call, for example, has a pay-off of the form ( 1
Sudu−K)+. With our previous
notations, it corresponds to the case α = 0, β = 1
and f(x) = (x−K)+.
The change of variables we used above is no longer suitable because it starts from zero when α = 0.
Instead, we consider the following new definition of the process ξ
eσ(Wt−Wu)+γ(t−u)du
ξ0 = S0.
Obviously, the two variables ξT and
Sudu have the same law. Hence, the price of the Asian option
becomes
C0 = E
e−rT f
e−rTf(ξT )
Remark 7 — The pricing of floating strike Asian options is also straightforward using this method. It is
even more natural to consider these options since it unveils the appropriate change of variables as we shall
see below.
Let us consider a floating strike Asian call for example. We have to compute
C0 = E
Sudu− ST
Using S̃t = Ste
δt as a numéraire (see the seminal paper of Geman et al. [9]), we immediately obtain that
C0 = EP
−δT ( 1
du− 1
where P
is the probability measure associated to the numéraire S̃t. It is defined by its Radon-Nikodym
derivative
= eσWT−
Under P
, the process Bt =Wt − σt is a Brownian motion and we can write that
C0 = EP
−δT ( 1
eσ(Bu−BT )+(r−δ+
)(u−T )du − 1
−δT ( 1
eσ(Wu−WT )+(r−δ+
)(u−T )du− 1
ξT − S0
where ξt is the process defined by (14) but with γ = r − δ + σ
. We see therefore that the problem simplifies
to the fixed strike Asian pricing problem.
Let us write down the stochastic differential equation that rules the process (ξt)t∈[0,T ]. Using Itô’s lemma,
we get {
dξt =
ξ0−ξt
dt+ ξt
σdWt + (γ +
ξ0 = S0.
Note that we are faced with a singularity problem near 0 because of the term ξ0−ξt
. We are going to reduce
its effect using another change of variables.
Using Itô’s lemma, we show that
C0 = E
e−rTf
where Xt = log(ξt/ξ0) solves the following SDE
dXt = σdWt + γdt+
e−Xt−1
X0 = 0.
Lemma 8 — Existence and strong uniqueness hold for the stochastic differential equation (16).
Proof : Existence is obvious since we have a particular solution Xt. The diffusion coefficient being
constant and the drift coefficient being a decreasing function in the spatial variable, we have also strong
uniqueness for the SDE (see the proof of Proposition 4). ✷
Because of the singularity of the term e
−Xt−1
in the drift coefficient, the law of (Xt)t≥0 is not absolutely
continuous with respect to the law of (σWt)t≥0. That is why we now define (Zt)t≥0 by the following SDE
with an affine inhomogeneous drift coefficient :
dZt = σdWt + γdt−
Z0 = X0 = 0.
The drift coefficient exhibits the same behavior as the one in (16) in the limit t → 0 in order to ensure the
desired absolute continuity property. It is affine in the spatial variable so that (Zt)t≥0 is a Gaussian process
and as such is easy to simulate recursively.
Lemma 9 — The process
sdWs +
t (18)
is the unique solution of the stochastic differential equation (17).
Proof : Using Itô’s Lemma, we easily check that Zt given by (18) is a solution of (17). Again, constant
diffusion coefficient and decreasing drift coefficient ensures strong uniqueness. ✷
Remark 10 — For the computation of the price C0 = E
e−rT (S0e
XT −K)+
of a standard Asian call
option, the random variable e−rT (S0e
ZT − K)+ provides a natural control variate. Indeed, since ZT is a
Gaussian random variable with mean γ
T and variance σ
, one has
e−rT (S0e
ZT −K)+
= S0e
−r)TN
−Ke−rTN (d)
where N is the cumulative standard normal distribution function and d = log(S0/K)+
Notice that in Kemna and Vorst [13], the authors suggest the use of the control variate
S0 exp
σWt + γt dt
which has the same law than e−rT
ZT −K
σWt + γt dt is also a Gaussian variable with mean
T and variance σ
In order to define a new probability measure under which (Zt)t≥0 solves the SDE (16), one introduces
Lt = exp
e−Zs − 1 + Zs
dWs −
e−Zs − 1 + Zs
Because of the singularity of the coefficients in the neighborhood of s = 0, one has to check that the integrals
in Lt are well defined. This relies on the following lemma
Lemma 11 — Let ǫ > 0. In a random neighborhood of s = 0, we have
|Zs| ≤ cs
−ǫ and |Xs| ≤ cs
where c is a constant depending on σ,γ and ǫ.
Since ∀ǫ > 0,
∀z ≤ cs 12−ǫ,
e−z − 1 + z
≤ Cs−4ǫ,
we can choose ǫ < 1
to deduce that Lt is well defined.
Proof : We easily check that the Gaussian process (Bt)t∈[0,T ] defined by Bt =
∫ (3t) 13
sdWs is a standard
Brownian motion. Thanks to the law of iterated logarithm for the Brownian motion (see for example Karatzas
and Shreve [12] p. 112), there exists t1(ω) such that
∀t ≤ t1(ω), |Bt(ω)| ≤ t
4ω is an element of the underlying probability space Ω.
Therefore,
∀t ≤ (3t1(ω))
3 , |Zt(ω)| =
(ω) +
∣∣ ≤ σ
Taking c = max( σ
) yields
∀t ≤ (3t1(ω))
3 ∧ 1, |Zt(ω)| ≤ ct
On the other hand, recall that Xt = log(ξt/ξ0) = log
eσWt+γt
e−σWu−γudu
. So, using the law of
iterated logarithm for the Brownian motion, we deduce that there exists t2(ω) such that
∀t ≤ t2(ω), 0 ≤
eσWt(ω)+γt
e−σWu(ω)−γudu ≤ 1
−ǫ−γudu.
Denote g(t) = 1
−ǫ−γudu and let us investigate the order in time near zero of this
function. We have that
+γt = 1 + σt
−ǫ +O(t1−2ǫ)∫ t
−ǫ−γudu = t+
−ǫ +O(t2−2ǫ)
hence
g(t) = 1 + (σ +
−ǫ +O(t1−2ǫ),
so Xt(ω) ≤ log (g(t)) ∼
−ǫ, which ends the proof for Xt. ✷
Proposition 12 — (Lt)t∈[0,T ] is a martingale and, consequently, for all g : C([0, T ]) → R measurable, the
random variables g((Xt)0≤t≤T ) and g((Zt)0≤t≤T )LT are simultaneously integrable and then
g((Xt)0≤t≤T )
g((Zt)0≤t≤T )LT
Proof : The proof is similar to the proof of Proposition 4.
We have already shown existence and strong uniqueness for both SDE (16) and (17). Showing that the
stopping time
τn(Y ) = inf
t ∈ R+ such that
e−Ys − 1 + Ys
ds ≥ n
, with the convention inf{∅} = +∞,
have infinite limits when n tends to +∞, QX and QZ almost surely, follows from the previous lemma.
One has
LT = exp
e−Zt − 1 + Zt
dZt −
e−Zt − 1 + Zt
e−Zt − 1 + Zt
+ γ − Zt
Set A(t, z) =
1− z + z2
− e−z
. The function A : ]0, T ]×R→ R is continuously differentiable in time and
twice continuously differentiable in space. So, we can apply Itô’s Lemma on the interval [ǫ, T ] for ǫ > 0 :
A(T, ZT ) = A(ǫ, Zǫ) +
e−Zt − 1 + Zt
dZt −
1− Zt + Z
− e−Zt
1− e−Zt
Using the lemma 9, we let ǫ→ 0 to obtain
A(T, ZT ) =
e−Zt − 1 + Zt
dZt −
1− Zt + Z
− e−Zt
1− e−Zt
LT = exp
A(T, ZT )−
φ(t, Zt)dt
where φ is the mapping
φ(t, z) =
e−z − 1 + z − z
1− e−z
e−z − 1 + z
e−z − 1 + z
+ γ − z
. (19)
By (15) and Proposition 12, we get
C0 = E
e−rTf(S0e
ZT ) exp
A(T, ZT )−
φ(t, Zt)dt
. (20)
Since for each t > 0, lim
φ(t, z) = +∞ and lim
φ(t, z) = −∞, it is not possible to apply the exact
algorithm. One can use the unbiased estimator, at least theoretically, if there exists a random variable cZ
measurable with respect to Z such that
eA(T,ZT )−(r+cZ)T |f(S0eZT )|e
|cZ−φ(t,Zt)|dt
Unfortunately, this reinforced integrability condition is never satisfied :
Lemma 13 — Assume that f is a non identically zero function. Let pZ and qZ denote respectively a positive
probability measure on N and a positive probability density on [0, T ]. Let N be distributed according to pZ
and (Ui)i∈N∗ be a sequence of independent random variables identically distributed according to the density
qZ , both independent conditionally on the process (Zt)t∈[0,T ]. Then the random variable
eA(T,ZT )−rT f(S0e
pZ(N)N !
−φ(Ui, ZUi)
qZ(Ui)
is non integrable.
Proof : By conditioning on Z, one has
∆ := E
eA(T,ZT )−rT |f(S0eZT )|
pZ(N)N !
|φ(Ui,ZUi )|
qZ (Ui)
eA(T,ZT )−rT |f(S0eZT )|e
|φ(t,Zt)|dt
eA(T,ZT )−rT |f(S0eZT )|e
|φ(t,Zt)|dt
One can easily show that, ∀z < 0 and ∀t ∈ [T
, T ], φ(t, z) ≥ φ(z) where
φ(z) =
e−z − 1 + z − z2
e−z − 1 + z
e−z − 1 + z
+ γ+ − 2 z
Since φ(z) ∼
σ2T 2
, there exists c < 0 such that for all z < c, φ(z) ≥ e−2z
σ2T 2
. Hence,
∆ ≥ E
eA(T,ZT )−rT |f(S0eZT )|e
e−2Zt1{Zt<c}dt
eA(T,ZT )−rT |f(S0eZT )|e−
2σ2T e
e−2Ztdt
Using Jensen’s inequality we get
∆ ≥ E
eA(T,ZT )−rT |f(S0eZT )|e−
2σ2T exp
We have seen in the proof of lemma 11 that Zt =
t where (Bt)t≥0 is a standard Brownian motion.
So, conditionally on ZT ,
Ztdt is a gaussian random variable and hence ∆ = +∞.
We are in a situation where eA(T,ZT )−rT |f(S0eZT )|E
[∣∣∣ 1pZ(N)N !
−φ(Ui,ZUi )
qZ (Ui)
∣∣∣(Zt)t∈[0,T ]
is non inte-
grable while eA(T,ZT )−rT |f(S0eZT )|
pZ (N)N !
−φ(Ui,ZUi )
qZ(Ui)
∣∣∣(Zt)t∈[0,T ]
]∣∣∣ is integrable since
e−rT |f(S0eZT )| exp
A(T, ZT )−
φ(t, Zt)dt
<∞. Then, a natural idea would consist in considering,
for a given n ∈ N∗, the random variable
eA(T,ZT )−rT |f(S0eZT )|E
∣∣∣∣∣∣
pZ(Nj)Nj !
−φ(U ji , ZUj
∣∣∣∣∣∣
∣∣∣(Zt)t∈[0,T ]
where (Nj)1≤j≤n are independent variables having the same law as N and
i )i∈N∗
1≤j≤n
are independent
sequences having the same law as (Ui)i∈N∗ , both independent conditionally on the process (Zt)t∈[0,T ]. The
following general result tells us that this is not sufficient to circumvent integrability problems.
Lemma 14 — Let Y and Z be two real random variables and g : R → R a given measurable function.
Assume that g(Z)E (Y |Z) is integrable while g(Z)E (|Y | |Z) is non integrable. Then, when (Yi)1≤i≤n is
a sequence of independent random variables having the same law as Y , ∀n ∈ N∗, the random variable
g(Z)E
i=1 Yi| |Z
is non integrable.
Proof : Denote by e, e1 and en three functions satisfying
∀z ∈ R, e(z) = E (Y |Z = z) , e1(z) = E (|Y1| |Z = z) and en(z) = E
(∣∣∣∣∣
∣∣∣∣∣
|Z = z
On the one hand, since
|g(z)| |e(z)|PZ(dz) <∞ and
|g(z)| e1(z)PZ(dz) = +∞ , where PZ is the law of
Z, we have that
|g(z)| e1(z)1{e1(z)≥2|e(z)|}PZ(dz) = +∞.
On the other hand, ∀z ∈ R,
en(z) ≥
(∣∣∣∣∣
∣∣∣∣∣
1{∀2≤j≤n,Yj≥0}|Z = z
(∣∣∣∣∣
∣∣∣∣∣
1{∀2≤j≤n,Yj<0}|Z = z
Y +1 |Z = z
P (Y1 ≥ 0|Z = z)n−1 + E
Y −1 |Z = z
P (Y1 < 0|Z = z)n−1
e1(z)+e(z)
P (Y1 ≥ 0|Z = z)n−1 + e1(z)−e(z)2 P (Y1 < 0|Z = z)
e1(z)
1{e1(z)≥2|e(z)|}P (Y1 ≥ 0|Z = z)
e1(z)
1{e1(z)≥2|e(z)|}P (Y1 < 0|Z = z)
≥ e1(z)
1{e1(z)≥2|e(z)|}
Hence, E
g(Z)E
(∣∣ 1
i=1 Yi
∣∣ |Z
|g(z)|en(z)PZ(dz) = +∞. ✷
There is still hope yet. In the proof of Lemma 13, we saw that integrability problems appear when
Zt takes large negative values so that φ(t, Zt) tends rapidly towards +∞. Since lim
φ(t, z) = −∞, one
possible issue is to split the function φ(t, Zt) into a positive part and a negative part. The first term can be
handled by the exact simulation technique whereas the second term, which as we shall see in the following
section presents no integrability problems, can be handled by the unbiased estimator technique.
2.2.1 An hybrid pseudo-exact method
We rewrite (20) in the following form
C0 = E
eA(T,ZT )−rTf(S0e
ZT )e
φ−(t,Zt)dte−
φ+(t,Zt)dt
. (22)
Let pZ and qZ denote respectively a positive probability measure on N and a positive probability density on
[0, T ]. Let N be distributed according to pZ and (Ui)i∈N∗ be a sequence of independent random variables
identically distributed according to the density qZ , both independent conditionally on the process (Zt)t∈[0,T ].
Note that, since eA(T,ZT )−rTf(S0e
ZT )e
|φ−(t,Zt)|dte−
φ+(t,Zt)dt = eA(T,ZT )−rT f(S0e
ZT )e−
φ(t,Zt)dt is
integrable, one has
C0 = E
eA(T,ZT )−rT f(S0e
pZ(N)N !
φ−(Ui, ZUi)
qZ(Ui)
φ+(t,Zt)dt
. (23)
Remark 15 — There is no hope that this estimator is square integrable. Indeed, one can show as in Lemma
13 that E
−(t,Zt))
= +∞ since (φ−(t, z))2 is of order z4 for large positive z.
The idea then is to apply the exact simulation technique to simulate an event with probability e−
φ+(t,Zt)dt.
Since for each t > 0, lim
φ+(t, z) = +∞, one needs to bound from above φ+(t, z), uniformly with respect
to t ∈ [0, T ], for z > c where c < 0 is a given constant. Thanks to the following lemma, it is possible to do
so but only uniformly with respect to t ∈ [ǫ, T ] for all ǫ > 0 :
Lemma 16 — For all 0 < t ≤ T ,
φ+(t, z) ≤ γ
∀c < 0, sup
z∈[c,0]
φ+(t, z) ≤ e
−c − 1 + c
(1 + γ+t) +
(e−c − 1)2
2σ2t2
Proof : Let z > 0. It is useful to distinguish two cases according to the sign of γ :
1. γ ≥ 0
We rewrite φ in the following form
φ(t, z) =
e−z − 1 + z − z
1− e−z
2 − (z ∧ 1)2
2σ2t2
(e−z − 1)2 − (z ∧ 1)2
2σ2t2
First note that
e−z−1+z− z2
≤ 0, 1−e
(e−z−1)2−(z∧1)2
2σ2t2
≤ 0. Moreover,
2 − (z ∧ 1)2
2σ2t2
if γt ≤ 1
otherwise
Consequently, φ+(t, z) ≤ γ
2. γ ≤ 0
Now we rewrite φ in the following form
φ(t, z) =
e−z − 1 + z − z
e−z − 1 + z
(e−z − 1)2 − z2
2σ2t2
1− e−z
It is then easy to show that φ+(t, z) ≤ 1
Note that 1
. Hence, gathering the two cases yields the first part of the lemma.
Let now z ∈ [c, 0] for a given negative constant c. We rewrite φ in the following form
φ(t, z) =
e−z − 1 + z
(1 + γ+t) +
(e−z − 1)2
2σ2t2
σ2t2︸ ︷︷ ︸
≥0 for z<0
1− e−z
− γ− e
−z − 1 + z
σ2t︸ ︷︷ ︸
≤0 for z<0
Since ∂z
e−z − 1 + z
(1 + γ+t) +
(e−z − 1)2
2σ2t2
1− e−2z − 2z + tγ+(1− e−z)
is negative for all
z < 0, one has that
z∈[c,0]
φ+(t, z) ≤ e
−c − 1 + c
(1 + γ+t) +
(e−c − 1)2
2σ2t2
This lemma suggests to apply the exact algorithm on [ǫ, T ] for a fixed positive threshold ǫ. It remains
to handle the time interval [0, ǫ[. Thanks to the following lemma, we that φ+(t, Zt) can be approximately
bounded from above for small t, almost surely, by a function of t. The idea is then to extend the exact
simulation algorithm by simulating an inhomogeneous Poisson process. Of course, this hybrid method is no
longer exact since the positive threshold for which the upper bound holds is random.
Lemma 17 — For all η > 0, there exists a random neighborhood of t = 0 such that
φ+(t, Zt) ≤
−η (24)
where c = max( σ
Proof : We rewrite (19) this way
φ(t, z) =
1− e−z
e−z − 1 + z
1− z + z2
− e−z − 1
(e−z − 1 + z)(e−z − 1− z)
and make the following Taylor expansions
1− z + z
− e−z − 1
(e−z − 1 + z)(e−z − 1− z)
z3 +O(z4)
1− e−z
e−z − 1 + z
z +O(z2).
On the other hand, we have seen in the proof of lemma 11 that there exists a random neighborhood of zero
such that Zt ≤ ct
−η where c = max( σ
). We conclude that, in a random neighborhood of zero,
φ+(t, Zt) ≤
2.2.2 Numerical computation
For numerical computation, we are going to use the following set of parameters : S0 = 100, K = 100,
σ = 0.2, r = 0.1, δ = 0 and T = 1. To fix the ideas, let us consider a call option. The price C0 writes as
follows
C0 = E
eA(T,ZT )−rT
ZT −K
φ−(Ui, ZUi)
φ+(t,Zt)dt
where N ∼ P(cp) and (Ui)i≥1 is an independent sequence of independent random variables uniformly dis-
tributed in [0, T ]. The parameter cp > 0 is set to one in the following. We give a description of the hybrid
method we implement :
Algorithm 2
On the time interval Ij := [
1. Simulate Z T
, Z T
and a lower bound mj for the minimum of (Zt)t∈Ij (use the fact that Zt =
t where (Bt)t≥0 is a standard Brownian motion).
2. Find M j > 0 such that ∀t ∈ Ij , φ+(t, Zt) ≤M j (use Lemma 16).
3. Simulate an homogeneous spatial Poisson process on the rectangle Ij × [0,M j] and accept (respectively
reject) the trajectory simulated if the number of points falling below the graph (φ+(t, Zt))t∈Ij is equal
to (respectively different from) zero.
Carry on this acceptance rejection algorithm until reaching a time interval IJ for a chosen J ∈ N∗. On the
remaining time interval [0, T
], use the same acceptance/rejection algorithm but with an inhomogeneous
spatial Poisson process this time (use Lemma 17).
In table 3, we give the price obtained by our method for different values of the positive threshold ǫ = T
The number M of Monte Carlo simulations is equal to 105 and the true price is equal to 7.042 (computed
using a Monte Carlo method with a trapezoidal scheme and a Kemna-Vorst control variate technique).
Price CPU
ǫ = T
6.9394 7s
ǫ = T
6.9590 10s
ǫ = T
6.9703 13s
ǫ = T
6.9952 17s
ǫ = T
7.0423 21s
Table 3: Price of the Asian call using the hybrid-pseudo exact method.
Clearly, the method is not yet competitive regarding computation time. Nevertheless, unlike the usual
discretization methods, it is not prone to discretization errors.
3 Conclusion
In this article, we have applied two original Monte Carlo methods for pricing Asian like options which
have the following pay-off : (αST + β
Stdt −K)+. In the case α 6= 0, we applied both the algorithm of
Beskos et al. [1] and a method based on the unbiased estimator of Wagner [23] and more recently the Poisson
estimator of Beskos et al. [2] and the generalized Poisson estimator of Fearnhead et al. [6]. The numerical
results show that the latter performs the best. The more interesting case α = 0, which corresponds to usual
continuously monitored Asian options, can not be treated using neither the exact algorithm, nor the method
of exact computation of expectation but we investigate an hybrid pseudo-exact method which combines the
two techniques. More generally, this hybrid method is an extension of the two exact methods and can be
applied in other situations.
From a practical point of view, the main contribution of these techniques is to allow Monte Carlo pricing
without resorting to discretization schemes. Hence, we are no longer prone to the discretization bias that
we encounter in standard Monte Carlo methods for pricing Asian like options. Even though these exact
methods are time consuming, they provide a good and reliable benchmark.
References
[1] A. Beskos, O. Papaspiliopoulos, and Gareth O. Roberts. Retrospective exact simulation of diffusion
sample paths. Bernoulli, 12(6), December 2006.
[2] A. Beskos, O. Papaspiliopoulos, Gareth O. Roberts, and Paul Fearnhead. Exact and computationally
efficient likelihood-based estimation for discretely observed diffusion processes. To appear in the Journal
of the Royal Statistical Society, Series B.
[3] M. Broadie and P. Glasserman. Estimating security price derivatives using simulation. Management
Science, 42(2):269–285, 1996.
[4] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E Knuth. On the lambert W function.
Advances in Computational Mathematics, 5:329–359, 1996.
[5] F. Dubois and T. Lelievre. Efficient pricing of asian options by the pde approach. Journal of Compu-
tational Finance, 8(2), 2004.
[6] Paul Fearnhead, O. Papaspiliopoulos, and Gareth O. Roberts. Particle filters for partially observed
diffusions. Working paper. Lancaster University., 2006.
[7] M. Fu, D. Madan, and T. Wang. Pricing continuous asian options : a comparison of monte carlo and
laplace transform inversion methods. Journal of Computational Finance, 2(2), 1999.
[8] H. Geman and A. Eydeland. Domino effect. Risk, pages 65–67, April 1995.
[9] H. Geman, N. El Karoui, and J.C. Rochet. Changes of numéraires, changes of probability measure and
option pricing. J. Appl. Probab., 32(2):443–458, 1995.
[10] H. Geman and M. Yor. Bessel processes, asian option and perpetuities. Mathematical Finance, 3(4),
1993.
[11] J.E. Ingersoll. Theory of Financial Decision Making. Rowman & Littlefield, 1987.
[12] I. Karatzas and Steven E. Shreve. Brownian motion and stochastic calculus. Springer-Verlag New-York,
second edition, 1991.
[13] A. Kemna and A. Vorst. A pricing method for options based on average asset values. Journal of Banking
and Finance, 14(1):113–129, 1990.
[14] B. Lapeyre and E. Temam. Competitive Monte Carlo methods for pricing asian options. Journal of
Computational Finance, 5(1), 2001.
[15] E. Levy. Pricing european average rate currency options. Journal of International Money and Finance,
11(5):474–491, October 1992.
[16] R. Lord. Partially exact and bounded approximations for arithmetic Asian options. Journal of Com-
putational Finance, 10(2), 2006.
[17] D. Revuz and M. Yor. Continuous martingales and Brownian motion. Springer-Verlag Berlin Heidelberg,
1991.
[18] L. C. G. Rogers and Z. Shi. The value of an Asian option. J. Appl. Probab., 32(4):1077–1088, 1995.
[19] T. H. Rydberg. A note on the existence of unique equivalent martingale measures in a markovian
setting. Finance and Stochastics, 1(3):251–257, 1997.
[20] S. Turnball and L. Wakeman. A quick algorithm for pricing european average options. Journal of
Financial and Quantitative Analysis, 16:377–389, 1991.
[21] J. Vecer. A new pde approach for pricing arithmetic asian options. Journal of Computational Finance,
4(4), 2001.
[22] T. Vorst. Prices and hedge ratios of average exchange rate options. International Review of Financial
Analysis, 1(3):179–193, 1992.
[23] W. Wagner. Unbiased Monte Carlo evaluation of certain functional integrals. J. Comput. Phys.,
71(1):21–33, 1987.
[24] W. Wagner. Monte Carlo evaluation of functionals of solutions of stochastic differential equations.
Variance reduction and numerical examples. Stochastic Anal. Appl., 6(4):447–468, 1988.
[25] W. Wagner. Unbiased multi-step estimators for the Monte Carlo evaluation of certain functional inte-
grals. J. Comput. Phys., 79(2):336–352, 1988.
[26] W. Wagner. Unbiased Monte Carlo estimators for functionals of weak solutions of stochastic differential
equations. Stochastics Stochastics Rep., 28(1):1–20, 1989.
[27] B. Wong and C. C. Heyde. On the martingale property of stochastic exponentials. J. Appl. Probab.,
41(3):654–664, 2004.
4 Appendix
4.1 The practical choice of p and q in the U.E method
The best choice for the probability law p of N and the common density q of the variables (Vi)i≥1
is obviously the one for which the variance of the simulation is minimum. In a very general setting, it is
difficult to tackle this issue. In order to have a first idea, we are going to restrict ourselves to the computation
p(N)N !
g(Vi)
q(Vi)
where g : [0, T ] → R.
Lemma 18 — When g is a measurable function on [0, T ] such that 0 <
|g(t)|dt < +∞, the variance of
p(N)N !
g(Vi)
q(Vi)
is minimal for
qopt(t) =
|g(t)|
|g(t)|dt
1[0,T ](t) and popt(n) =
|g(t)|dt
|g(t)|dt
Proof : Minimizing the variance in (7) comes down to minimizing the expectation of the square of
p(N)N !
g(Vi)
q(Vi)
F (p, q) = E
(p(N)N !)2
g2(Vi)
q2(Vi)
g2(t)
p(n) (n!)2
Using Cauchy-Schwartz inequality we obtain a lower bound for F (p, q)
F (p, q) =
g2(t)
p(n)n!
p(n) ≥
g2(t)
q(t)dt
|g(t)|dt
= exp
|g(t)|dt
We easily check that this lower bound is attained for qopt and popt.
The optimal probability distribution popt is the Poisson law with parameter
|g(t)|dt. This justifies our
use of a Poisson distribution for p.
4.2 Simulation from the distribution h given by (13)
Recall that
h(u) = C exp
A(u)− (u−X0)
= C exp
(1− e−σu)− (u−X0)
where C is a normalizing constant.
The expansion of the exponential e−σu at the first order yields
h(u) ≈ C exp
u− (u−X0)
= C exp
(u− (X0 + T (γ+βS0)σ ))
This suggests to do rejection sampling using the normal distribution with mean X0 +
T (γ+βS0)
variance T as prior. Unfortunately, for a standard set of parameters, this method gives bad results. Even a
second order expansion of e−σu which also modifies the variance does not work.
In order to get round this problem, we evaluate the mode u∗ of h. We have
h′(u∗) = C
∗ −X0
(1− e−σu
)− (u
∗ −X0)2
So, h′(u∗) = 0 if and only if
∗ −X0
which writes
σ(u∗ −X0 −
T )eσ(u
∗−X0− γσT ) = TβS0e
−σX0−γT .
The function x 7→ xex is continuous and increasing on [0,+∞[ and so is its inverse which we denote by
W . Since TβS0e
−σX0−γT ≥ 0, we deduce that h is unimodal and that its mode satisfies
γT +W
βS0Te
−γT−σX0
+ σX0
The functionW is the well-known Lambert function, also called the Omega function. It is uniquely valued
on [0,+∞[ and there are robust and fast numerical methods based on series expansion for approximating
this function (see for example Corless et al. [4]).
Numerical tests showed that performing rejection sampling using a Gaussian distribution with variance
T and mean u∗ instead of X0 +
T (γ+βS0)
gives plain satisfaction. In table 4, we see that for arbitrary choice
of the parameter α
, the acceptance rate of the algorithm is always high (of order 70%) and that the
computation time is low.
Nb of simulations Acceptance rate Computation time
0.2 61% 3s
0.5 106 68% 3s
0.8 80% 2s
Table 4: Acceptance rate of the rejection algorithm of simulating from the distribution h in (13) with
S0 = 100, σ = 0.3, T = 2 and r = 0.1.
	1 Exact Simulation techniques
	1.1 The exact simulation method of Beskos et al. BPR
	1.2 The unbiased estimator (U.E)
	2 Application : the pricing of continuous Asian options
	2.1 The case =0
	2.1.1 Numerical computation
	2.2 Standard Asian options : the case = 0 and > 0
	2.2.1 An hybrid pseudo-exact method
	2.2.2 Numerical computation
	3 Conclusion
	Bibliography
	4 Appendix
	4.1 The practical choice of p and q in the U.E method 
	4.2 Simulation from the distribution h given by (13)
