Stochastic expansion for the pricing of call options with

discrete dividends
P. tor and E. Gobet
Grenoble University and CNRS
Laboratoire Jean Kuntzmann
51 rue des Mathmatiques
F - 38041 Grenoble cedex 9
pierre.etore@imag.fr, emmanuel.gobet@imag.fr
March 1, 2012

Abstract. In the context of an asset paying affine-type discrete dividends, we present

closed analytical approximations for the pricing of European vanilla options in the Black-
Scholes model with time-dependent parameters. They are obtained using a stochastic
Taylor expansion around a shifted lognormal proxy model. The final formulae are respec-
tively first, second and third order approximations w.r.t. the fixed part of the dividends.
Using Cameron-Martin transformations, we provide explicit representations of the correc-
tion terms as Greeks in the Black-Scholes model. The use of Malliavin calculus enables
us to provide tight error estimates for our approximations. Numerical experiments show
that the current approach yields very accurate results, in particular compared with known
approximations of [BGS03, VW09], and quicker than the iterated integration procedure of
[HHL03] or than the binomial tree method of [VN06].

Keywords: equity option, discrete dividend, stochastic approximation, analytic formula.

This research was partially supported by University Joseph Fourier (MSTIC grant REFINE) and by
Fondation du Risque - Chaire Risques Financiers Socit Gnrale - Ecole Polytechnique - ENPC (second
author).
The authors thank Hans Buehler, Adil Reghai and Mike Staunton for their valuable comments, and Haibiao
Huang for pointing us errors regarding the numerical results in the first version of this work.

1
29 pages 2

Classification. MSC (2000): 41XX - 60HXX. JEL: G13.

Premia 14
Introduction
Usually, stocks pay dividends, which modeling is a non-trivial issue. This has also some
implications regarding the computational point of view, that is to efficiently price vanilla
options written on the stock and to quickly calibrate the stock model. If we use a deter-
ministic continuously paid dividend yield and assume that the asset dynamics is lognormal,
then we can extend the classical pricing Black-Scholes formula by minor modifications (see
equation (1.6)). Assuming continuous dividends is an approximation that can be justified
if one considers a large portfolio of stocks paying individually discrete dividends. However,
for a single stock, considering discrete dividends is more realistic and this is our frame-
work. Actually, our aim is to provide efficient approximation formulae for Call options
written on a single asset paying discrete dividends. For this, we follow an approach based
on stochastic expansions, using stochastic analysis tools, approach that has been similarly
developed in a series of papers [BGM09, BGM11, BGM10b, BGM10a].
In the literature, several works handle the issues of numerical computation of the call
price when dividends are discrete. Of course, a Monte Carlo approach is still possible,
whatever the asset model and the dividend model are, but usually it is not competitive
compared with analytical approximations or one-dimensional tree methods. Several works
[HHL03, VN06, VW09] rely on the dynamic programming equation between two successive
dividend dates, say ti and ti+1 . Namely, denote by C(t, S) the option price function at
time t for an asset equal to S, write di(S) for the (known) dividend policy modeling the
dividend paid at time ti (it depends on the asset): then, for a Markovian price process
(St )0tT and deterministic interest rates (rt )0tT we have
R ti+1
rs ds
C(ti , Sti ) = E(e ti
C(ti+1 , St di+1 (St ))|Sti ), (0.1)
i+1 i+1

the expectation being computed under the risk-neutral pricing measure. In [HHL03], the
authors discuss in details the proper choice of dividend policy. In addition, they compute
the price function C(., .) using integration methods to compute the expectation in (0.1),
for tractable dynamics of S (lognormal for instance). This numerical approach is exact
(up to integration error) but it is computationally intensive. In [VN06], for a piecewise
lognormal asset, the authors design a binomial tree method to solve (0.1). The main
difficulty in using a tree method is the a priori non-recombination of the nodes at the
dividend dates. The authors overcome this problem by using interpolation techniques
29 pages 3

between nodes. They also prove the convergence of their approximation, as the number of
steps in the tree method goes to infinity. Finally, in [VW09], still for a piecewise lognormal
model and for a fixed dividend policy di+1 (S) = i+1 , the authors expand the equality (0.1)
w.r.t. (i )i and provide an approximation formula involving the Black-Scholes formula
and its Greeks w.r.t. the spot. For n dividend dates, the number of BS price/greeks to
compute grows exponentially like 3n , as the number of dividend dates increase (in their
tests, the authors take n = 7, giving 2187 terms to evaluate); it may be very costly.
Another approach is developed by Bos, Gairat and Shepeleva in [BGS03]: they give an
approximation formula for the equivalent implied Black-Scholes volatility, in order to take
into account the dividends. It is obtained by a suitable average of the instantaneous
volatility of the asset paying dividends. After the submission of this article, we have been
aware of another work [SG10] based on a matching method in the case of fixed dividend
policy and lognormal model.
In this work, we derive an alternative expansion of the price function w.r.t. the div-
idends. The resulting approximation also writes as a combination of BS formulae and
Greeks w.r.t. the strike (and not the spot). Compared with [VW09], our second order
approximation formula requires the evaluation of only 45 BS price/greeks for 7 dividend
dates. Thus, at least regarding the computational cost, it improves [VW09] and it is similar
to [BGS03]. Moreover, our assumption on the dividend policy is less restrictive, see below.
In addition, the numerical results show an excellent accuracy of our formulae.
In the current work, the model for S is a piecewise lognormal model (with time-
dependent parameters) and the dividend policy is affine in S, i.e. including a fixed and a
proportional part:
di (S) = i + yiS.
One drawback of this model is that after a dividend payment, the asset price may become
negative because the relation di (S) S may be violated for small S. However, in most
of our numerical tests, the probability of such event is very small (see Tables 5 and 6);
presumably, it has a very small impact on the call price. Although this model of dividend
policy is quite simple, it is often used by practitioners. We also mention that for such
affine dividends, Buehler [Bue10] gives the necessary form of stochastic model for S, that
is coherent with the no-arbitrage assumption; this excludes our piecewise lognormal model.
Actually, our purpose is rather to expose our approximation method in a simple case, before
extending it (in further works) to more general models and piecewise affine dividends (that
are coherent with no-arbitrage, see [Reg09] for instance).
To obtain our approximations, we choose a model proxy obtained by averaging the
future dividends. Then we use stochastic expansion techniques in the spirit of the work
[BGM09, BGM11, BGM10b, BGM10a]. A significative part of effort is made to derive non
asymptotic error estimates, justifying the first order, second or third order approximation.
29 pages 4

This approach is quite flexible and we believe that this work paves the way for future
research in order to obtain analytical approximations of call price with discrete dividends
including Heston or local volatilities, or stochastic interest rate as well.
The organization of the paper is the following. In the next section, we define the model
and notations used throughout the work. In Section 2, we state our main approximation
results about first, second and third order approximation formulae for the call price. Ex-
tensions to the computation of the Delta are given as well. Section 3 is devoted to the
proof of technical results involving Malliavin calculus. Numerical tests are presented in
Section 4.

1 Model and notations

1.1 Financial framework
We consider a standard complete financial market, with a traded risky asset on which an
European vanilla option with maturity T is written. In our study, specifically the asset
pays dividends at known dates 0 < t1 < . . . < tn T < + (n 1). We assume that the
second date t2 (whenever existing when n > 1) is larger than one year (t2 1): actually,
this assumption is only technical, it is made in order to get uniform error estimates w.r.t.
t2 (see the proof of Propositions 2.1-2.2 and Lemma 3.2). For quarterly paid dividends e.g.,
our expansion formulas are still valid but error estimates are slightly modified. We skip
details. At time ti , the amount of dividends is split into a proportional part yi [0, 1[ and
a fixed part i 0. To make clear the asset dependency w.r.t. the dividends, we denote
(y,)
by (St )t the asset price process. Then, the amount of dividend at time ti equals to
(y,)
i + yi Sti ,

which implies that the asset price jumps downwards to

(y,) (y,) (y,) (y,)
Sti = Sti [i + yi Sti ] = Sti (1 yi ) i (1.1)

just after the dividend payment.

Moreover, we assume that between two dividend dates, the asset follows an Ito dynamics
with a time-dependent volatility (t )t . Since we focus only pricing/hedging issues, we write
the dynamics of S (y,) under the (unique) risk-neutral measure Q: between two dividend
dates it writes
(y,) (y,) (y,)
dSt = t St dWt + (rt qt )St dt
where W is a standard Q-Brownian motion. In the above equation, (qt )t should be inter-
preted as a (deterministic) repo rate. The interest rate (rt )t is assumed to be deterministic.
The functions (rt )t and (qt )t are bounded.
29 pages 5

1.2 Assumptions and notations

Assumptions. It is not a practical restriction to assume that the ratio between fixed
dividends and the current asset S0 remains bounded by a constant c (likely smaller than
1 in practice): supi i /S0 c .
In addition, for some of our results we may assume that for the first dividend date,
the ratio is small enough in the sense
1
< 1. (D)
S0 (1 y1 )

For some results, we impose a non-degeneracy condition on the model (ellipticity

condition):

t [0, T ], 0 < t . (E)

Notations. For convenience, we repeatedly use the following notations.

Rt
We write Dt for the discount factor: Dt = exp( 0 (rs qs )ds).
Rt Rt
We write Mt = exp( 0 s dWs 12 0 s2 ds) for the log-normal martingale with volatility
(s )0sT .

We write S for the (fictitious) asset without dividends:

dSt = t St dWt + (rt qt )St dt,
(y,)
and its initial value is S0 = S0 . Thus,
Mt
St = S0 . (1.2)
Dt
n
We set i,n := (1 yj ) = (1 yn ) (1 yi+1 ) for 0 i n with the convention
Y

j=i+1
that j=n+1 (1 yj ) = 1 (so that n,n = 1).
Qn

For the sake of conciseness, we may use the simplified notation

Dti
i = i i,n . (1.3)
DT
For a given strike K > 0, the shifted strike K (y,) will play an important role in our
approximation formulae:
n n
Dti
K (y,) = K + =K+ (1.4)
X X
i i,n i .
i=1 DT i=1
29 pages 6

We write A c B when A cB for some constant c which depends smoothly on

the model parameters. The constant c remains bounded as the maturity T or the
parameters r, q, , y, c go to 0. The constant c may depend on the ratio / 1
and on the number n of dividend dates. The dependency w.r.t. S0 is systematically
written. When relevant, explicit dependencies w.r.t. parameters are indicated.

In the error analysis, we use repeatedly the standard estimates E(supsT Ssk ) c S0k
for any k R.

1.3 Preliminary relations

With the previous notations, we easily deduce



St if t < t1 ,
(y,) St (y,) St St


(y,)
St = Sti = (1 yi)Sti i if ti t < ti+1 for i < n (1.5)

Sti Sti Sti
or ti t T for i = n.

Then, an easy induction (detailed in Appendix) leads to the following

n
(y,) ST
= 0,n ST
X
Lemma 1.1. We have ST i i,n .
i=1 Sti
A case of special interest corresponds to proportional dividends only (i 0) for which
(y,0)
we have ST := 0,n ST . This is a lognormal random variable and an explicit formula for
the related call price is available via the Black-Scholes formula:
RT
(y,0)
E(e 0
rs ds
(ST K)+ ) = CallBS (0,n S0 , K) (1.6)
with RT RT
CallBS (x, k) = xe 0 qs ds N [d+ (x, k)] ke 0 rs ds N [d (x, k)] ,
1 x 1 Z T
1
 
d (x, k) = qR log + qR (rs qs s2 )ds,
T 2 k T 2 0 2
0 s ds 0 s ds
N being the cumulative distribution function of a standard Gaussian variable. Note that
the price depends of course of (rt )t , (qt )t and (t )t , but we choose in our notations to
highlight the dependency w.r.t. the initial value and the strike. Indeed this plays a crucial
role in our calculations.
The case i 0 is important for our study since it serves to find a proxy for the case
with fixed dividends. The proxy will not be directly given by the model with 0, but
by this model shifted by the expectation of the fixed dividends. In other words, in view of
Lemma 1.1, the proxy is defined by
n n n
(y,)
X ST  Dti
:= 0,n ST E = 0,n ST = 0,n ST (1.7)
X X
ST i i,n i i,n i ,
i=1 Sti i=1 DT i=1
29 pages 7

D
recalling the definition i = i i,n DtTi . This is a shifted lognormal random variable, thus
RT
(y,)
the computation of E(e 0 rs ds (ST K)+ ) is still explicit, by taking the shifted strike
variable K (y,) (defined in (1.4)) in the Black-Scholes formula:
RT
(y,)
 
rs ds
E(e 0 (ST K)+ ) = CallBS 0,n S0 , K (y,) . (1.8)
RT
(y,)
The above quantity stands for the main term of our expansion formula of E(e 0 rs ds (ST
K)+ ) (see Theorems 2.3 and 2.4). The asymptotics underlying the expansion is supi i /S0
0 (small fixed dividends).
Our next purpose is now twofold: first, to provide correction terms, that will enable us
to achieve a remarkable accuracy. Second, to give tight error estimates w.r.t. the model
parameters.

2 Main results
Our analysis is based on Taylor expansions and smart computations of the correction terms
(y,)
using the proxy ST . In order to study the distance to the proxy, we use Lemma 1.1 and
equality (1.7) to write
n
(y,) (y,0) Dti MT
=ST (1 + 1)
X
ST i i,n
i=1 DT Mti
n
(y,) MT
=ST i ( 1). (2.1)
X

i=1 Mti
RT
(y,)
Our ultimate purpose is to approximate E(e 0 rs ds h(ST K)) for h(x) = x+ (that is
the Call price). Actually, the derivation of the approximation and the error estimation are
simpler when the function h is smoother than for Call/Put option. We start by this case
in subsection 2.1, for the convenience of the reader. Then, handling call payoffs h(x) = x+
requires more technicalities related to Malliavin calculus and we tackle this case later in
subsection 2.2 and section 3.

2.1 Taylor expansion for smooth h

The degree k ( 1) of smoothness of h is defined as follows:

(Hk ) The function h(.) is (k1)-times continuously differentiable and the (k1)-th deriva-
tive is almost everywhere differentiable. Moreover, the derivatives are polynomially
bounded: for some positive constants C and p one has |h(x)| + kj=1 |xj h(x)|
P

C(1 + |x|p ) for any x R.

29 pages 8

RT
(y,)
First order approximation. We aim at approximating E(e 0 rs ds h(ST K)) for
functions h satisfying (H2 ) . Using a first order Taylor expansion we have
h RT i h RT i
(y,) (y,)
E e 0
rs ds
h(ST K) = E e 0
rs ds
h(ST K)
n
" #

RT
rs ds (y,) MT
(ST 1) + Error2 (h) (2.2)
X
i E e 0 h K)(
i=1 Mti

where |Error2 (h)| c (1 + S0p ) supi (i k M

MT
ti
1k3 )2 . By standard computations (see also
Lemma 3.3), we have
MT q
k 1kp cp T ti (2.3)
Mti

for any p 1. It readily follows that |Error2 (h)| c (1 + S0p ) supi (i T ti )2 .
It remains to simplify the terms in the summation of (2.2). For each 1 i n we
write
h RT
(y,) MT i h RT
(y,) MT
E e 0 rs ds h (ST K)( 1) =E e 0 rs ds h (ST K) ]
Mti Mti
h RT
(y,)
E e 0
rs ds
h (ST K)]. (2.4)

The second term on the right hand side can be rewritten using a derivative w.r.t. K (the
assumptions on h allow us to interchange derivation and expectation):
RT RT
(y,) (y,)
E[e 0
rs ds
h (ST K)] = K E[e 0
rs ds
h(ST K)]
RT 
= k E[e rs ds
h(0,n ST k)] . (2.5)
0

k=K (y,)

This representation is useful for the call/put case to interpret expansion terms as Greeks
(and thus explicit terms).
Note that we have in general, for any multiplicative constant > 0, any strike k, and
any derivative of order m N of any sufficiently smooth function h,
RT RT
rs ds (m)
E[e 0 h (ST k)] = (1)m km E[e 0
rs ds
h(ST k)]. (2.6)
MT
Regarding to the first term in the r.h.s. of (2.4), we interpret the factor Mti
as a change
Rt
of measure on FT . Under the new induced measure Q , Wt =R Wt i
0 s 1ti sT ds is a
T
s2 ds
Brownian motion. Then, ST under Qi has the same law as ST e ti
under Q. Thus,
h RT
(y,) MT i
E e 0
rs ds
h (ST K) (2.7)
Mti
RT RT n
s2 ds
=E[e rs ds
h (0,n e
X
0 ti
ST i K)]
i=1
RT RT
s2 ds

rs ds
= k E[e h(0,n e ti
ST k)]
0

k=K (y,)
29 pages 9

using (2.6) at the last line. Combining the above equality with (2.5) and (2.4), and plugging
this into (2.2), we obtain our first main result.

Theorem 2.1. For a smooth function h satisfying (H2 ) , we have

h RT i
(y,)
E e 0
rs ds
h(ST K)
h RT i
=E e 0
rs ds
h(0,n ST K (y,) )
n RT RT
s2 ds
h i
+ i k E e rs ds
X
h(0,n e ti
ST k) 
0

k=K (y,)
i=1
h RT i 
k E e rs ds
h(0,n ST k)  + Error2 (h), (2.8)
0

k=K (y,)

with |Error2 (h)| c(1 + S0p ) supi (i T ti )2 .

Note that in the terms on the r.h.s. of the above equality, the function h is system-
atically evaluated at a shifted lognormal random variable. This allows for simple and
tractable one-dimensional numerical computations.
This approximation formula is a first-order expansion formula w.r.t. the fixed dividends
since the error is a O(supi i2 ).

Second order approximation. Applying the same kind of arguments, we can derive
another formula, which residual terms are of order three w.r.t. the fixed dividends.

Theorem 2.2. For a smooth function h satisfying (H3 ) , we have

h RT i
(y,)
E e 0
rs ds
h(ST K)
h RT i
=E e 0
rs ds
h(0,n ST K (y,) )
n RT RT
s2 ds
h i
rs ds
+
X
i k E e h(0,n e ti
ST k) 
0

k=K (y,)
i=1
RT !
h i
rs ds
k E e h(0,n ST k) 
0

k=K (y,)

1
RT RT RT
s2 ds+ s2 ds s2 ds
h RT i
rs ds
+ i j k2 E
X
e h(0,n e ti tj
ST k)  e ti tj
0

2 1i,jn
k=K (y,)

n n RT RT
s2 ds
X X h i
2 j i k2 E e rs ds
h(0,n e ti
ST k) 
0

k=K (y,)
j=1 i=1
n
!
X 2 h RT i
rs ds
+ j k2 E e h(0,n ST k)  + Error3 (h), (2.9)
0

k=K (y,)
j=1

with |Error3 (h)| c(1 + S0p ) supi (i T ti )3 .
29 pages 10

Proof. The proof is similar to that of Theorem 2.1, except that the equality (2.2) is replaced
by a second order Taylor expansion. It gives
h RT i h RT i
(y,) (y,)
E e 0
rs ds
h(ST K) = E e 0
rs ds
h(ST K)
n
" #

RT
rs ds (y,) MT
(ST 1)
X
i E e 0 h K)(
i=1 Mti
!2
n
1 RT
(y,)
hM
T
i
+ E e 0 rs ds h (ST K) 1
X
i
2 i=1 Mti
+ Error3 (h)
where Error3 (h) c (1+S0p ) supi (i k M
MT
ti
1k4 )3 . Then, the announced estimate on Error3 (h)
easily follows by using (2.3).
Compared with the expansion in Theorem 2.1, it remains to transform the new contri-
bution with the factor 12 . This term is equal to
!2
n n
RT
(y,) MT X
E e rs ds
h (ST K)
X
0 i i
i=1 Mti i=1
" #

RT
rs ds (y,) MT MT
= (ST
X
i j E e 0 h K)
1i,jn Mti Mtj

n n
" #

RT
rs ds (y,) MT
2 (ST
X X
j i E e 0 h K)
j=1 i=1 Mti
2
n  RT 
rs ds (y,)
+ (ST
X
j E e 0 h K)
j=1

:=T1 + T2 + T3 .
We handle separately each of the three terms above.
Term T1 . We proceed analogously to the equality (2.7) by transforming this term via
MT MT
different changes of probability measure. Indeed, note that M = exp( 0T s (1ti sT +
R
t Mt i j
1 RT RT
1tj sT )dWs 2 ti
[s (1ti sT
+ 1tj sT )] ds) exp( defines (up to the sec- 2 2
ti tj s ds)
ond exponential factor) a change of measure Q under which (Wt 0t s (1ti sT +
i,j
R
(y,)
1tj sT )ds)t0 is a Brownian motion. It means that ST under Qi,j has the same
RT RT
s2 ds+ s2 ds
law as 0,n ST e under Q. Thus, we obtain
Pn
ti tj
l=1 l
n
" RT RT # RT
s2 ds+ s2 ds s2 ds
RT
rs ds
T1 = h (0,n e
X X
i j E e 0
ti tj
ST l K) e ti tj

1i,jn l=1
RT T R RT
2 ds+ s2 ds s2 ds
h RT i
rs ds
= i j k2 E ti s
X
e h(0,n e tj
ST k)  e ti tj
,
0

k=K (y,)
1i,jn
29 pages 11

using (2.6) at the last line.

Term T2 . Similarly, we obtain

n n RT RT
s2 ds
h i
T2 = 2 i k2 E e rs ds
X X
j h(0,n e ti
ST k)  .
0

k=K (y,)
j=1 i=1

Term T3 . Clearly, we have

2
n h RT i
T3 = j k2 E e rs ds
X
h(0,n ST k)  .
0

k=K (y,)
j=1

The theorem is proved.

2.2 Expansion results for call payoff

We now extend the previous results from smooth functions h to the call option function
h(x) = x+ , using a regularization argument that is quite standard. However one has to
be careful with the error estimates since they depend on h or h in the previous case of
smooth functions. To safely pass to the limit, we impose the non-degeneracy condition
(E) on the model. The assumption (D) enables us to get error estimates uniform in t1 , as
t1 goes to 0.
We first precise the derivatives of CallBS (x, k) with respect to the strike k. We have
RT
k CallBS (x, k) = e 0
rs ds
N (d (x, k)), (2.10)
RT
e 0
rs ds
1 2
k2 CallBS (x, k) = q R e 2 d (x,k), (2.11)
T 2
k 2 0 s ds
and RT
e rs ds
d (x, k)
 
0 1 2
k3 CallBS (x, k) = q RT e 2 d (x,k) qR 1 .
T 2
k2 2 0 s2 ds 0 s ds
We now state our main results, giving a first and second order formula for the price of a
Call option written on a multi-dividend asset (a third order formula is given in Subsection
2.4).
29 pages 12

Theorem 2.3. Assume (D) and (E). We have

RT
(y,)
E(e 0
rs ds
(ST K)+ )
 
=CallBS 0,n S0 , K (y,)
n RT !
BS s2 ds (y,) BS (y,)
+ i k Call (0,n S0 e ) k Call (0,n S0 , K )
X
ti
,K
i=1

+ Error2 (Call), (2.12)

 q 2
with |Error2 (Call)| c supi i
S0
1 ti
T
S0 T .

The result below states a second order approximation result.

Theorem 2.4. Assume (D) and (E). We have

RT
(y,)
E(e 0
rs ds
(ST K)+ )
 
=CallBS 0,n S0 , K (y,)
n RT !
BS s2 ds (y,) BS (y,)
+ i k Call (0,n S0 e ) k Call (0,n S0 , K )
X
ti
,K
i=1
1
RT RT RT
s2 ds s2 ds+ s2 ds
+ k2 CallBS (0,n S0 e , K (y,) )
X
i j e ti tj ti tj

2 1i,jn
n n RT
s2 ds
X X
2 j i k2 CallBS (0,n S0 e ti
, K (y,) )
j=1 i=1
n
!
X 2
+ j k2 CallBS (0,n S0 , K (y,) ) + Error3 (Call), (2.13)
j=1

 q 3
with |Error3 (Call)| c supi i
S0
1 ti
T
S0 T .

To state the error estimates, we have taken a specific form which allows to assert that
our approximation error is of order one or two w.r.t. supi i /S
R0
. This is especially clear for
R T T
At-The-Money options, for which 0,n S0 e 0 qs ds = K (y,) e 0
rs ds
. In that case, using the

Brenner-Subrahmanyam approximation [BS88] CallBS (x, k) = 12 x( v + o(v)) as v =

k=x
RT  q 2
2
0 s ds goes to 0, we obtain that the relative ATM error is bounded by c supi i
S0
1 ti
T
 q 3
(in Theorem 2.3) or c supi Si0 1 tTi (in Theorem 2.4). This indicates that the relative
accuracy of our approximation depends mainly of the ratio supi i /S0 and not much of the
other parameters. This is confirmed by the numerical results (see Section 4).
The results for put option are simply obtained by replacing the CallBS (.) function by
the PutBS (.) function. Then we observe that these approximations verify the Call-Put
parity relation.
29 pages 13

2.3 Proof of Theorems 2.3 and 2.4

The sketch of the proof is the following: we take a sequence of smooth functions (hN )N
converging to h(x) = x+ in a suitable sense. Then, the proof is divided in two steps.
1. First, we prove that the expansion terms computed with hN converge to those com-
puted with h.

2. Second, we estimate the limsup of the error terms Error2 (hN ) and Error3 (hN ) as N
goes to infinity.
In this subsection, we only give details regarding Step 1. Step 2, involving Malliavin
calculus, is much more technical. We postpone it to the next section.
The justification of the Step 1 relies on the following Lemma.
Lemma 2.1. Assume (E), take > 0 and k > 0. Consider a sequence of measurable func-
tions (hN )N 1 and h having a polynomial growth uniformly in N, i.e. for some constants
C > 0 and p > 0 we have supxR |hN(1+|x|
(x)|+|h(x)|
p) C.
h i h i
i) Then, the functions k 7 E hN (ST k) and k 7 E h(ST k) are infinitely
continuously differentiable on ]0, [.

ii) In addition, assume that hN converges almost everywhere to h as N goes to infinity.

Then, for any m N, we have
h i h i
lim km E hN (ST k) = km E h(ST k) . (2.14)
N

Proof. Under (E), the law of ST has an explicit density w.r.t. the Lebesgue measure. It
gives
qR
S0 x T
exp(x2 /2)
RT
s2 ds 12
Z
s2 ds
h i
E h(ST k) = h( e 0 0
k) dx
R DT 2
Z
= h(z)p(z + k)dz
k
RT RT
exp([log(uDT /(S0 ))+ 12 s2 ds]2 /[2 s2 ds])
where p(u) = 1u>0 q RT
0 0
. It is easy to check that p is
u 2 s2 ds
0
infinitely continuously differentiable on ]0, [ and that its derivatives at u = 0 are equal
R m
to 0. In addition, for
h any m N, i we have 0 |u p(u)|(1 + |u| )du < . These properties
p

easily imply that E h(ST k) is smooth w.r.t. k (i.e. statement i)) and that
h i Z
km E h(ST k) = h(z)zm p(z + k)dz.
k

From this representation and by an application of the dominated convergence theorem,

statement ii) readily follows.
29 pages 14

A sequence of functions (hN )N 1 satisfying (H3 ) and converging to h(x) = x+ .

Rx Ru
For N N , take hN defined by hN (x) = N(1 N|v|)+ dvdu: each function hN
satisfies (H3 ) and it is easy to check the following properties

i) hN (x) = hN (x) = 0 for x 1/N,

ii) hN (x) = x and hN (x) = 1 for x 1/N,

iii) 0 hN (x) 1,

iv) (hN )N converges uniformly to h as N .

Owing to the above uniform convergence, we have

RT RT
(y,) (y,)
lim E(e 0
rs ds
hN (ST K)) = E(e 0
rs ds
(ST K)+ ),
N
RT RT
lim E(e 0
rs ds
hN (0,n ST K (y,) )) = E(e 0
rs ds
(0,n ST K (y,) )+ )
N
 
= CallBS 0,n S0 , K (y,)

using the Black-Scholes formula (1.8) for the last equality. Moreover, using Lemma 2.1,
we obtain for any > 0
RT RT
lim k E(e 0
rs ds
hN (ST k)) = k E(e 0
rs ds
(ST k)+ )
N
 
= k CallBS S0 , k .

Thus, we can apply Theorem 2.1 with hN and pass to the limit as N goes to infinity. It
gives the expansion of Theorem 2.3, with

lim Error2 (hN ) = Error2 (Call).

N

However, the upper bounds on Error2 (hN ) given in Theorem 2.1 involve hN and it does
not enable us to pass to the limit on the error estimates. In the next section, we prove
specific estimates using Malliavin calculus:

Proposition 2.1. Assume (D) and (E). Then, we have

2

s
i ti
|Error2 (hN )| c sup 1 S0 T ,
i S0 T

uniformly in N. Consequently, the same estimate applies to Error2 (Call).

29 pages 15

Using the above result, the proof of Theorem 2.3 is complete.

Similarly, for the second derivative, we have
RT  
lim k2 E(e 0
rs ds
hN (ST k)) = k2 CallBS S0 , k
N

Analogously to Proposition 2.1, we have

Proposition 2.2. Assume (D) and (E). Then, we have

3

s
i ti
|Error3 (hN )| c sup 1 S0 T ,
i S0 T

uniformly in N. Consequently, the same estimate applies to Error3 (Call).

Thus, we complete the proof of Theorem 2.4 as for Theorem 2.3.

2.4 Extension to the third-order approximation price formula

Using similar techniques we can state a third order formula. We leave the details of the
proof to the reader.
29 pages 16

Theorem 2.5. Assume (D) and (E). We have

RT
(y,)
E(e 0
rs ds
(ST K)+ )
 
=CallBS 0,n S0 , K (y,)
n RT !
BS s2 ds (y,) BS (y,)
+ i k Call (0,n S0 e ) k Call (0,n S0 , K )
X
ti
,K
i=1
1
RT RT RT
s2 ds s2 ds+ s2 ds
+ k2 CallBS (0,n S0 e , K (y,) )
X
i j e ti tj ti tj

2 1i,jn
n n RT
s2 ds
X X
2 j i k2 CallBS (0,n S0 e ti
, K (y,) )
j=1 i=1
n
!
X 2
+ j k2 CallBS (0,n S0 , K (y,) )
j=1

1
RT RT RT
s2 ds+ s2 ds+ s2 ds
+
X
i j l e ti tj ti tl tj tl

6 1i,j,ln
RT RT RT
s2 ds+ s2 ds+ s2 ds
k3 CallBS (0,n S0 e ti tj tl
, K (y,) )
n RT RT RT
X  s2 ds s2 ds+ s2 ds
3 k3 CallBS (0,n S0 e , K (y,) )
X
j i j e ti tj ti tj

j=1 1i,jn
n n
2 X RT
s2 ds
X
+3 j i k3 CallBS (0,n S0 e ti
, K (y,) )
j=1 i=1
n
!
X 3
j k3 CallBS (0,n S0 , K (y,) )
j=1

+ Error4 (Call),
 q 4
with |Error4 (Call)| c supi i
S0
1 ti
T
S0 T .

2.5 Extension to the approximation of the Delta

Adapting our methodology we can also derive several expansion formulas for the delta
of a Call option on a multi-dividend asset. We choose to present only the second order
approximation formula.
Let us first fix some extra notations. With the convention that t0 = 0 we set
RT
s2 ds
i,n := i,n e ti
, 0 i n,
n
Dti
i := i i,n

, 1 i n and K (y,,) := K + i .
X
DT i=1
29 pages 17

RT
(y,)
Theorem 2.6. Assume (D) and (E). Let = S0 E(e 0 rs ds (ST K)+ ) be the Delta
of the Call option of strike K on the multidividend asset. We have
(
= 0,n x CallBS (0,n S0 , K (y,,) )
n  RT
s2 ds
+ i k,x
2
CallBS (0,n S0 e , K (y,,) )
X
ti

i=1

2
k,x CallBS (0,n S0 , K (y,,))
1
RT RT RT
s2 ds s2 ds+ s2 ds

+ i j 3
CallBS (0,n S0 e , K (y,,) )
X
e ti tj
k,k,x ti tj

2 1i,jn
n n RT
s2 ds
X  X
2 j i 3
k,k,x CallBS (0,n S0 e ti
, K (y,,) )
j=1 i=1
n
X 2 )
+ j 3
k,k,x CallBS (0,n S0 , K (y,,) )
j=1

+ Error3 (Digital),
 q 3
with |Error3 (Digital)| c supi i
S0
1 ti
T
.
Remark 2.1. The third order error is denoted Error3 (Digital) because, up to multiplicative
R T
constants, the Delta is of the form, = E(e 0
rs ds
1S (y,,) >K ), that is the price of a Digital
T
Call option with S (y,,) as an asset to be described in the following sketch of proof.
RT
Remark 2.2. We recall that x CallBS (x, k) = e 0 qs ds N (d+ (x, k)). The other greeks are
easily computed from (2.10), (2.11). We skip details.
Proof. We only give the main lines. Details can be treated adapting the proofs of Theorems
2.2 and 2.4.

(y,)
Step 1. Taking into account Lemma 1.1, the pathwise derivative of ST w.r.t. S0 is
ST 0,n
0,n S0 = DT MT . Thus, interchanging derivation and expectation,
0,n RT
= E(e 0 rs ds 1S (y,) >K MT ).

DT T

Again we interpret MT as a change of measure on FT . Under the new induced measure

(y,)
Q0 , Wt = Wt 0t s ds is a brownian motion. Then ST under Q0 has the same law as
R

n n
ST ST
RT RT
(y,,) s2 ds s2 ds
:= 0,n ST e = 0,n (2.15)
X X
ST 0 i i,n e ti
ST i i,n ,
i=1 Sti i=1 Sti
under Q. Thus,
0,n RT
(y,,)
= E(e 0 rs ds h(ST

K)), (2.16)
DT
with h(x) = 1x>0 .
29 pages 18

Step 2. From (2.16), we see that the evaluation of the Delta is reduced to that of the
price of a digital Call written on an asset S (y,,) with new dividend parameters (i,n

)i
(compare Lemma 1.1 and (2.15)). Then, the derivation of an approximation formula is
similar to what we have done in Theorem 2.2 and 2.4. Briefly, we take a sequence of
smooth functions (hN := 12 (tanh(N .) + 1))N converging to h almost everywhere, and we
apply Theorem 2.2 and Lemma 2.1. It gives
"
0,n h R T rs ds
i
= Ee 0 h(0,n ST K (y,,) )
DT
n RT RT
s2 ds
h i
+ i k E e rs ds
X
h(0,n e ti
ST k) 
0

k=K (y,,)
i=1
RT !
h i
rs ds
k E e h(0,n ST k) 
0

k=K (y,,)

1
RT RT RT
s2 ds+ s2 ds s2 ds
h RT i
rs ds
+ i j k2 E
X
e 0 h(0,n e ti tj
ST k)  e ti tj

2 1i,jn
k=K (y,,)

n n RT RT
s2 ds
X X h i
2 j i k2 E e 0
rs ds
h(0,n e ti
ST k) 
k=K (y,,)
j=1 i=1
n
!#
X 2 h RT i
+ j k2 E e rs ds
h(0,n ST k)  + lim Error3 (hN ).
0

k=K (y,,) N
j=1

Similarly to Proposition 2.2, it is possible to prove that

s 3
i ti
|Error3 (hN )| c sup 1 , (2.17)
i S0 T

uniformly in N, using |hN | = 1 (see Remark 3.1). This gives the upper bound for
Error3 (Digital).

Step 3. It remains to relate the correction terms to the Black-Scholes formula. Actually,
for any multiplicative constant > 0 and any positive strike k, we have

0,n h R T rs ds i RT
0 qs ds
Ee 0 S >k = 0,n e
10,n N (d (0,n S0 , k))
DT T
RT
= 0,n e 0
qs ds
N (d+ (0,n S0 , k))
= 0,n x CallBS (0,n S0 , k).

Then, by successive differentiation w.r.t. k, we obtain the announced formula.

29 pages 19

3 Proof of Propositions 2.1 and 2.2

In the proof of Theorems 2.1 and 2.2, we have obtained that
n 2 Z !
MT 1
RT X
rs ds
Error2 (hN ) = E e 0 i ( 1) (1 )hN (FT K)d , (3.1)
i=1 Mti 0
n
(1 )2
3 Z !
MT 1
RT X
rs ds
Error3 (hN ) = E e 0 i ( 1) hN (FT K)d (3.2)
i=1 Mti 0 2
where 0 1 we define
n
(y,) Dti MT
FT := ST ( 1). (3.3)
X
i i,n
i=1 DT Mti

3.1 Technical results from Malliavin calculus

Our aim is to provide tight error estimates on supN |Error2 (hN )| and
supN |Error3 (hN )|, using |hN | = 1. For this, we use the Malliavin calculus integration by
parts to transform the above expectations. It requires the use of several lemmas stated
hereafter, that will be proved in the next subsection. The results deal with the Malliavin
MT
estimates of the random variable FT and ( M ti
1). Regarding to Malliavin calculus related
to the one-dimensional Brownian motion W , we freely adopt the notation from [Nua06].
For instance, the first Malliavin derivative of a random variable F is the H = L2 ([0, T ], dt)-
valued process denoted by DF = (Dt F )0tT . The second derivative takes values in H 2
and is denoted by D 2 F = (Ds,t 2
F )0s,tT , and so on. If the scalar product on H k is denoted
by h., .iH k , then the Malliavin covariance matrix of F is defined by F = hDF, DF iH . We
freely use the notation Dk,p (k 1, p 1) for the space of k-times Malliavin differentiable
random variables (with derivatives in Lp ) and related k.kk,p-norms (see [Nua06, Section
1.2]). We set Dk, = p1 Dk,p and D = k1 Dk, .
Lemma 3.1. For all p 1, all 0 1, FT is in D and
sup kDt (FT)kp cp S0 , 2
sup kDt,r (FT )kp cp S0 2 ,
tT t,rT
3
sup kDt,r,s (FT )kp cp S0 . 3
t,r,sT

Lemma 3.2. Assume (D) and (E). We have for all 0 1, F1 p1 Lp (). In
T
addition,
1
p 1, 0 1, kF1 k2,p cp 2 2 .
T S0 T
Lemma 3.3. Let (Nt )0tT be a Brownian martingale with a bounded bracket, and assume
that JT = NT 21 hNiT is in D2, . Then for all r 1, one has
 
keJT 1k2,r cr kJT k2,2r 1 + kJT k2,4r + kJT k21,8r e2r sup hN iT .
29 pages 20

We are now in a position to complete the proof of Propositions 2.1 and 2.2. Let us
start with Error2 (hN ): by Fubinis theorem, it is equal to
RT Z 1 h MT MT i
e rs ds
(1 )E ( 1)( 1)hN (FT K) d. (3.4)
X
0 i j
1i,jn 0 Mti Mtj

We now control (uniformly in and N) the above expectations. To remove the singularity
problem of hN , we apply an integration by parts of Malliavin calculus. Indeed, from [Nua06,
Proposition 2.1.4], one knows that for 1 i, j n and 0 1 there exists Hij1, D ,
MT MT
depending only on FT and ( M t
1)( M t
1), such that
i j

MT MT
   
E ( 1)( 1)hN (FT K) = E hN (FT K)Hij1, . (3.5)
Mti Mtj
MT MT
This is justified by the fact that ( M t
1)( M t
1) D , FT is in D and is non degenerate
i j
(Lemma 3.2) under the assumption (E). Our task then becomes to find an upper bound,
uniformly in , i, j, for kHij1, kp , for all p. Using the discussion in [Nua06, p.102] we have

MT MT
kHij1, kp cp kF1 k1,4p kD. (FT )k1,4p k( 1)( 1)k1,2p . (3.6)
T Mti Mtj

Considering Lemma 3.2 it remains to estimate the two last terms of the r.h.s. above. By
definition, we have
Z T Z T Z T q/2
kD. (FT )kq1,q = E|( |Dr (FT )|2 dr)1/2 |q + E 2
|Ds,r (FT )|2 dr ds .
0 0 0

Then by standard inequalities combined with Lemma 3.1, we get for any q

kD. (FT )k1,q cq S0 T . (3.7)

On the other hand, using Hlder type inequalities on k.kk,p -norms (see [Nua06, Proposition
1.5.6]), we have
MT MT MT MT
k( 1)( 1)k1,q c1,q k 1k1,2q k 1k1,2q .
Mti Mtj Mti Mtj

In order to apply Lemma 3.3, we define the Brownian martingale Nt = 0t s 1ti sT dWs
R

which bracket is bounded by 2 (T ti ). Using the notation of Lemma 3.3, notice that
1 MT

eJT = eNT 2 hN iT = M ti
. Clearly kJT k2,r cr T ti . Then from Lemma 3.3, it readily
follows that
MT q
k 1k2,r cr T ti . (3.8)
Mti
MT MT
q
In particular, it gives k( M t
1)( Mt
1)k 1,q c q 2
T ti T tj .
i j
29 pages 21

We combine the latter inequality with (3.6), (3.7), Lemma 3.2 and we obtain for any
i, j, p
q
T ti T tj
kHij1,kp cp ,
S0 T
uniformly in . Plugging this estimate into (3.5) (and using |hN | = 1) leads to
q
T t T tj
  
 MT M T i
E ( 1)( 1)hN (FT K)  c .
 
 Mti Mtj  S0 T

supi (i T ti )2
In view of (3.4), we have proved that |Error2 (hN )| c S0 T
. Proposition 2.1 is
proved.

The proof of Proposition 2.2 is very similar and we only give the main intermediate
estimates. Analogously to the identity (3.5), we have
RT 1 (1 )2
Z
rs ds
Error3 (h, N) =
X
e 0 i j l
1i,j,ln 0 2
MT h MT MT i
E ( 1)( 1)( 1)h
N (FT K) d
Mti Mtj Mtl
RT Z 1 (1 )2 h i
2,
= e 0 rs ds i j l E hN (FT K)Hijl (3.9)
X
.
1i,j,ln 0 2

Furthermore, applying the general estimates from [Nua06, p.102] combined with (3.8), we
obtain
2,
|Error3 (h, N)| c i j l sup kHijk
X
k1
1i,j,ln [0,1]
!
MT MT MT
i j l sup kF1 k22,8 kD. (FT )k22,8 1)( 1)( 1)k2,2
X
c k(
1i,j,ln [0,1] T Mti Mtj Mtl
1 2 3q q q
i j l 2 2 2 (S0 T ) T ti T tj T tl
X
c
1i,j,ln (S0 T )
q
c sup(i T ti )3 2 .
i S0 T
We are finished.

Remark 3.1. In the proof of Theorem 2.6 we have to control uniformly in N, |Error3 (hN )|,
with (hN )N approximating h(x) = 1x>0 , and satisfying |hN | = 1. Compared with the proof
of Proposition 2.2, we have to use
h MT MT MT i h
3,
i
E( 1)( 1)( 1)h
N (FT

K) = E hN (FT

K)H ijl .
Mti Mtj Mtl
29 pages 22

Indeed we cannot have a uniform control on hN here. Similarly to the proofs above, we
have
3,
sup[0,1] kHijk k1 c kF1 k33,16 kD. (FT )k33,16 k( M
MT
t
MT
1)( M t
MT
1)( M t
1)k3,2
T i j l

1 3 3 q
c (S 0 T ) T ti T tj T tl
(S02 2 T )3

and thus the error estimate (2.17) stated in Theorem 2.6.

3.2 Proof of technical lemmas

Proof of Lemma 3.1. Take a fixed . As FT is an affine function of the lognormal variables
MT
Mt
s (see (3.3)) it is clear that FT D . We have
i

n
MT Dt MT
Dt (FT ) = S0 0,n
X
t 1tT i i,n i t 1ti <tT ,
DT i=1 DT Mti
n
2 MT Dt MT
(FT ) = S0 0,n
X
Dt,r t r 1t,rT i i,n i t r 1ti <t,rT ,
DT i=1 DT Mti
n
3 MT Dt MT
(FT ) = S0 0,n
X
Dt,r,s t r s 1t,r,sT i i,n i t r s 1ti <t,r,sT .
DT i=1 DT Mti

Standard computations (using supi i c S0 ) lead easily to the announced estimates.

Proof of Lemma 3.2. It is enough to consider the case p 2. Take [0, 1].
Step 1. We first estimate F1 in Lp . We have
T

Z T
FT = kDFT k2H = |Ds (FT )|2 ds (3.10)
0
n
T i 2
Z 
ST2 2
X
0,n

i,n 1ti sT  ds
0 i=1 Sti
 2
 1 
ST2 2 0,n
2
t1 + 1  (T 1n=1 + t2 1n>1 t1 ) . (3.11)
 
(1 y1 )St 1


Thus, it is clear that

FT ST2 2 0,n
2
t1
inducing that FT is invertible and its inverse is in any Lp . Now, our aim is to estimate
the Lp -norm of F1 uniformly in t1 1. For this, we define the event
T

1
( )
1
A1 =
1+ (1 y1 )St1
29 pages 23

where the parameter will be set at a positive value close to 0. Then, on this event, we
have
1
1 1 > 0.
(1 y1 )St1 1+ 1+
Thus, on A1 , we obtain
FT 2
t1 + (T 1n=1 + t2 1n>1 t1 )
2
ST2 2 0,n (1 + )2
2 2
(T 1n=1 + t2 1n>1 ) (1 T )
(1 + )2 (1 + )2
using if n > 1 that t2 1.
We now estimate P(Ac1 ) by leveraging the assumption S0 (1 y1 ) < 1 . Using that St1
has a lognormal distribution, we obtain
1 (1 + )
!
P(Ac1) = P Mt1 < Dt1
S0 (1 y1 )
1 1 t1 2 1 (1 + )
" Z #!
= N R t1 2 1/2 log(Dt1 ) + ds + log( ) .
( 0 s ds) 2 0 s S0 (1 y1 )

Choose close to 0 enough to ensure that log( S01(1y

(1+)
1)
) = C < 0. Then, using N (x)
exp(x /2) for any x, we deduce
2

 #2
1 1
"
P(Ac1 ) exp 2 |r q| t1 + 2 t1 + C cp 1 ( 2 t1 )p ,
2 t1 2

for any p > 0.

Finally, bringing together our different estimates, we deduce
ST2 (1 + )2 ST2
0 F1 1A1 + 1Ac1 .
T
2
0,n 2 (1 T ) 2 2
0,n 2 t1

By Hlder inequalities, together with the fact that ST1 p1 Lp , we obtain

1 1 2 1
 
kF1 kp c 2 2 + 2 c , (3.12)
T S0 (1 T ) S02 2 T
possibly changing the value of the generic constant c at the last inequality. This proves
the first statement of the Lemma.

Step 2. We turn to estimate the Malliavin derivatives of F1 . By the chain rule, we

T
obtain
2
Ds F Ds,t F Ds FT Dt FT
Ds (F1 ) = 2 T and Ds,t
2
(F1 ) = 2 T + 2 .
T F T F F3
T T T
29 pages 24

On the one hand, by definition of the k.k2,2r -norms, we have

kF1 kp2,p = E(Fp ) + EkD(F1 )kpH + EkD 2 (F1 )kpH 2

T T T T

c E(Fp ) + E(F2p p
kDF kH )
T
+ E(F2p 2 p
kD F kH 2 )
T
+ E(F3p 2p
kDF kH ),
T
T T T T

that is

kF1 k2,p c kF1 kp + kF1 k24p kDFT kH

T T T 2p
2
+ kF1 k24p kD 2FT kH 2 + kF1 k36p kDFT kH . (3.13)
T 2p T 4p

One the other hand, using Minkowski and Hlder inequalities combined with Lemma 3.3,
from (3.10) we derive
Z T 1/2 Z T 1/2
kDF kH = |Dt FT |2 dt kDt FT k22p dt

T 2p 0 p 0
Z T Z T 1/2
k 2Ds FT Dt,s
2
FT dsk22p dt
0 0
Z T Z T 1/2
( 2
2kDs FT k4p kDt,s FT k4p ds)2 dt c S02 3 T 3/2 .
0 0

Similarly, we obtain
Z T Z T Z T 1/2
kD 2 F kH 2 ( 2
kDt,r [(Ds FT )2 ]k2p ds)2 dtdr c S02 4 T 2 .

T 2p 0 0 0

Plugging the above inequalities and the estimate (3.12) into (3.13) yields

kF1 k2,p c kF1 kp + kF1 k24p kDFT kH

T T T 2p
2
+ kF1 k24p kD 2 FT kH 2 + kF1 k36p kDFT kH .

T 2p T 4p
1 1 1
!
 
c 2 2 1 + 2 2 S02 3 T 3/2 + S02 4 T 2 + 2 2 2 (S02 4 T 2 )2
S0 T S0 T (S0 T )
1
c 2 2 .
S0 T
The proof is complete.

Proof of Lemma 3.3. The Taylor formula yields

Z 1
eJT 1 = JT euJT du.
0

Using Hlder and Minkowski inequalities, we obtain

Z 1
JT
ke 1k2,r c kJT k2,2r keuJT k2,2r du. (3.14)
0
29 pages 25

For any 0 u 1, we have

Ds euJT = ueuJT Ds JT 2 uJT

and Ds,t e = euJT (uDs,t
2
JT + u2 Ds JT Dt JT ).

Then, by definition of the k.k2,2r -norms, we obtain

keuJT k2r
2,2r = E([e
uJT 2r
] ) + EkDeuJT k2r 2 uJT 2r
H + EkD e k H 2
E(e2ruJT ) + E(e2ruJT kDJT k2r
H ) + E(e
2ruJT
(kD 2 JT kH 2 + kDJT k2H )2r )
c keuJT k2r 2r 4r
4r (1 + kJT k2,4r + kJT k1,8r ). (3.15)
1
Finally, since (epuNt 2 hpuN it )t defines an exponential martingale (for any fixed p), one has
1 1 1 2 p2
keuJT kpp = E[epu(NT 2 hN iT ) ] = E[epuNT 2 hpuN iT + 2 hN iT (pu+(pu) ) ] e 2 sup hN iT
.

Plugging this estimate into (3.15) and (3.14), we get the announced result.

4 Numerical experiments
In all our tests we use as benchmark a Monte Carlo price computed with 2.109 drawings,
and control variates (column "Monte Carlo" in the tables). The control variates consist in
European options with the same parameters except that 0 (see the discussion after
Lemma 1.1). In the tables the numbers between parentheses in the Monte Carlo columns
refer to the half width of the 95% confidence interval around the computed prices.
We wish first to compare our results with the ones of recent papers in the literature
(namely [BGS03, VN06, VW09]). In Table 1, the abbrevations EG3, VNRE, VN1000, VW
and BGS refer respectively to our method with the order three formula, the method of
Vellekoop and Nieuwenhuis with Richardson Extrapolation, their method without extrap-
olation and 1000 time steps (both in [VN06]), the method in [VW09] and the method in
[BGS03]. The example is the one treated in these last three papers: till time maturity
T = 7.0 we have 7 dividend payment dates 0 < t1 < . . . < t7 < T with ti+1 ti = 1 for
all 1 i < 7. We test the cases t1 = 0.1, 0.5 and 0.9. The successive i s are 6, 6.5, 7,
7.5, 8, 8 and 8. We have y 0. The coefficients (rt )t , (qt )t and (t )t are constant, with
r = 6%, q = 0% and = 25%. We take S0 = 100 and test the strikes K = 70, 100 and 130.

insert Table 1 about here

These tests show that the accuracy of our method is better than the one of VN1000
and BGS and similar to the one of VW. However the VNRE method seems to be the most
accurate.
29 pages 26

But note that, for n dividend payment dates, the number of terms to compute in our
order two and three formulae are respectively

(n + 2)(n + 3) (n + 2)(n2 + 7n + 12)

and .
2 6
Thus, the number of terms computed for EG3 in Table 1 is 165, which requires a small
computational time. Concerning the VNRE method the maximal number of time steps is
64000, which is fairly demanding. See also the discussion p13 in [VW09]: 2187 evaluations
of the Black-Scholes formula and any of its derivatives are computed to achieve the prices
reported in Table 1. In other words, from the computational point of view, our approach
is very competitive, compared with other existing methods.
We now test the sensitivity to the parameters of the precision of our option pricers.
As indicated by Theorems 2.3, 2.4 and 2.5, the error should increase with volatility, time
maturity and the amplitude of the i s.
In Table 2 we have r = 6%, q = 0% , S0 = 100, n = 3, and yi = 0.02 and i = 2 for all
1 i n. We have t1 = 0.5, ti+1 ti = 1 for all 1 i < n and T = 3.0. We successively
test = 15%, 25%, 45%, and compute the prices with the formulae at order one, two and
three (respectively EG1, EG2 and EG3) for various strikes. Under each price we report
the corresponding implied volatility (expressed in %).

Insert Table 2 about here

As the volatility increases we observe a loss of accuracy on the prices computed with
EG1, while for EG2 and EG3 the accuracy remains nearly the same. This suggests that
our method is quite robust to variations of the volatility.
In Table 3 we set = 25%, the other parameters as in Table 2, and test the influence
of the amplitude of the i s. We take i = for all 1 i n and test the values = 2, 6, 10.

Insert Table 3 about here

With 2 the results of EG1, EG2 and EG3 are accurate up to one basis point on
implied volatilities (even if EG1 seems to be a bit less accurate on the prices themselves).
With 6 both EG2 and EG3 match the implied volatilities, but we observe a slightly
difference of accuracy on the prices. With 10 only EG3 still performs well to match
prices and implied volatilities. Note that, as expected the solvers are always more accurate
at the money.
Finally, in Table 4 we investigate the influence of n, keeping = 25% constant, and
the other parameters as in Table 2, except 4. We choose n = 3, 5 , 10, which is related
to testing the influence of the maturity T = n.
29 pages 27

Insert Table 4 about here

With n = 3 the solvers EG2 and EG3 perfectly match the implied volatilities. The
solver EG1 is accurate up to 2 bp on implied volatilities, which is generally sufficient for
calibration purposes. As expected, with n = 10 a loss of accuracy can be observed (both
on prices and implied volatilities). Even with EG3 the implied volatilities can fail to match
the ones corresponding to Monte Carlo prices. Some computed prices are slightly outside
the Monte Carlo confidence interval (especially for in the money options).
Note that similar tests show no significative influence of the parameters yi on the results:
for = 0.25 and the other parameters as in Table 2 , EG2 and EG3 both match the implied
volatility with 0 bp error, whatever the value of the yi s.
Note also that we have used our Monte Carlo simulations to estimate the probabilities
(y,)
that ST < 0. Indeed, with the affine type dividend model there is no guarantee that
this never occurs. The numerical results show that this probability increases with and n
(see Tables 5 and 6). For n = 10 this probability is larger than 2% (in the results of Table
1 this estimated probability is also about 2%: indeed the dividends are of high amplitude
and n = 7). This suggests that the dividend model itself has to be refined as S (y,) is close
to zero.

Insert Tables 5 and 6 about here

5 Conclusion
In this work, we have derived approximation formulae for the vanilla option prices written
on an asset paying discrete dividends, under lognormality assumptions. Numerical tests
show that the second order approximation (Theorem 2.4) is accurate enough for usual val-
ues of the fixed part of dividends (that is few % of the spot value) and for maturities smaller
than five years. For larger dividends or longer maturities, the third order approximation
(Theorem 2.5) yields additional accuracy. Moreover, compared with other methods, these
expansions are quicker to evaluate (or as quick as [BGS03]). Finally, we mention several
possible extensions. Combining the stochastic expansion approaches recently developed
in [BGM09, BGM11, BGM10b, BGM10a] with the current work, we could generalize the
closed formulae to local or stochastic volatility models, including Gaussian stochastic in-
terest rates. This is left to further research.
29 pages 28

A Proof of Lemma 1.1

This is proved by induction. The result is true for n = 1, considering (1.5). Suppose it is
true for any n( 1) dates (ti )1in and consider that an extra dividend payment is made
at time tn+1 (tn , T ]. Then, we have

(y,) (y,) ST ST
ST =(1 yn+1)Stn+1 n+1
Stn+1 Stn+1
 n n  n
ST  Y Stn+1
  
=(1 yn+1) (1 yi ) Stn+1 (1 yj )
X Y
i
Stn+1 i=1 i=1 j=i+1 Sti
ST
n+1
Stn+1
 n+1 n  n S ST
T

= (1 yi ) ST i (1 yn+1 ) (1 yj )
Y X Y
n+1
i=1 i=1 j=i+1 Sti Stn+1
 n+1 n+1 n+1 S
T
 X
= (1 yi ) ST (1 yj )
Y Y
i .
i=1 i=1 j=i+1 Sti

Thus, it is proved for n + 1 dates.

References
[BGM09] E. Benhamou, E. Gobet, and M. Miri. Smart expansion and fast calibration
for jump diffusion. Finance and Stochastics, 13(4):563589, 2009. 2, 3, 27

[BGM10a] E. Benhamou, E. Gobet, and M. Miri. Expansion formulas for European options
in a local volatility model. International Journal of Theoretical and Applied
Finance, 13(4):603634, 2010. 2, 3, 27

[BGM10b] E. Benhamou, E. Gobet, and M. Miri. Time dependent Heston model. SIAM
Journal on Financial Mathematics, 1:289325, 2010. 2, 3, 27

[BGM11] E. Benhamou, E. Gobet, and M. Miri. Analytical formulas for local volatility
model with stochastic rates. To appear in Quantitative Finance, 2011. 2, 3, 27

[BGS03] R. Bos, A. Gairat, and D. Shepeleva. Dealing with discrete dividends. Risk
Magazine, pages 109112, 2003. 1, 3, 25, 27

[BS88] M. Brenner and M.G. Subrahmanyam. A simple formula to compute the im-
plied standard deviation. Financial Analysts Journal, 44:8083, Sept./Oct.
1988. 12
29 pages 29

[Bue10] H. Buehler. Volatility and dividends - Volatility modelling with

cash dividends and simple credit risk. Technical report, TU Berlin,
http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1141877, 2010. 3

[HHL03] E.G. Haug, J. Haug, and A. Lewis. Back to Basics: a New Approach to the
Discrete Dividend Problem. Wilmott Magazine, pages 3747, september 2003.
1, 2

[KR05] R. Korn and L. C. G. Rogers. Stocks paying discrete dividends: modelling and
option pricing. The Journal of Derivatives, 13(2):4448, 2005.

[Nua06] D. Nualart. Malliavin calculus and related topics. Springer Verlag, second
edition, 2006. 19, 20, 21

[Reg09] A. Reghai. Two factor stochastic volatility with em-

bedded local volatility. Technical report, Natixis,
http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1360962, 2009. 3

[SG10] F. Sahel and A. Gocsei. Matching the sensitivities to discrete dividends:

a new approach for pricing vanillas. Technical report, Socit Gnrale,
http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1663880, 2010. 3

[VN06] M. Vellekoop and J. Nieuwenhuis. Efficient Pricing of Derivatives on Assets

with Discrete Dividends. Applied Mathematical Finance, 13(3):265284, 2006.
1, 2, 25

[VW09] C. Veiga and U. Wystup. Closed formula for options with discrete dividends
and its derivatives. Applied Mathematical Finance, 16(6):517531, 2009. 1, 2,
3, 25, 26
29 pages 30

t1 K Monte Carlo EG3 VNRE VN1000 VW BGS

0.1 70 24.8962 (11.104 ) 24.8787 24.90 24.92 24.8862 24.71

100 17.4338 (12.104 ) 17.4255 17.43 17.46 17.4394 17.42
130 12.3994 (12.104 ) 12.396 12.40 12.43 12.4114 12.50

0.5 70 26.0806 (11.104 ) 26.0678 26.08 26.10 26.0752 25.87

100 18.4815 (12.104 ) 18.476 18.48 18.50 18.489 18.45
130 13.2844 (11.104 ) 13.283 13.29 13.31 13.2968 13.38

0.9 70 27.2341 (10.104 ) 27.205 27.21 27.23 27.2117 26.99

100 19.4817 (11.104 ) 19.4784 19.48 19.5 19.4905 19.43
130 14.1296 (10.104 ) 14.1293 14.13 14.16 14.1419 14.06

Table 1: European Call option prices, with = 25%, r = 6%, q = 0%, S0 = 100.
29 pages 31

K Monte Carlo EG1 EG2 EG3

0.15 40 55.3322 (104 ) 55.3319 55.3322 55.3322
15.63% 15.02% 15.63% 15.63%
60 38.7349 (104 ) 38.7315 38.735 38.735
15.55% 15.47% 15.55% 15.55%
80 23.4551 (104 ) 23.4511 23.4511 23.4552
15.49% 15.48% 15.49% 15.49%
100 12.0146 (2.104 ) 12.0127 12.0145 12.0145
15.46% 15.45% 15.46% 15.46%
120 5.3187 (104 ) 5.3164 5.3187 5.3187
15.43% 15.43% 15.43% 15.43%
180 0.2746 (5.105 ) 0.2727 0.27462 0.2746
15.38% 15.36% 15.38% 15.38%
250 0.00619 (105 ) 0.00604 0.00619 0.006187
15.34% 15.34% 15.34% 15.34%
0.25 40 55.4626 (104 ) 55.4551 55.4628 55.4627
26.08% 25.86% 26.08% 26.08%
60 40.0472 (104 ) 40.0386 40.0470 40.0472
25.92% 25.88% 25.92% 25.92%
80 27.3445 (2.104 ) 27.3398 27.3444 27.3446
25.83% 25.82% 25.83% 25.83%
100 17.9604 (2.104 ) 17.9574 17.9603 17.9604
25.77% 25.77% 25.77% 25.77%
120 11.5237 (2.104 ) 11.5203 11.5236 11.5236
25.73% 25.72% 25.73% 25.73%
180 2.90748 (104 ) 2.9028 2.90744 2.90745
25.64% 25.63% 25.64% 25.64%
250 0.60324 (5.105 ) 0.6007 0.60328 0.60325
25.58% 25.56% 25.58% 25.58%
0.45 40 57.6888 (2.104 ) 57.6699 57.6882 57.6888
47.02% 46.91% 47.02% 47.02%
60 45.9191 (2.104 ) 45.9088 45.9187 45.9191
46.73% 46.70% 46.73% 46.73%
80 36.7278 (2.104 ) 36.7213 36.7272 36.7278
46.56% 46.55% 46.56% 46.56%
100 29.599 (2.104 ) 29.5933 29.5988 29.599
46.45% 46.44% 46.45% 46.45%
120 24.0528 (2.104 ) 24.0468 24.0526 24.0528
46.36% 46.35% 46.36% 46.36%
180 13.5404 (2.104 ) 13.533 13.5403 13.5404
46.20% 46.18% 46.20% 46.20%
250 7.4835 (2.104 ) 7.4762 7.4834 7.4834
46.08% 46.07% 46.08% 46.08%

Table 2: European Call option prices, with r = 6%, q = 0%, S0 = 100, n = 3, y 0.02
and 2.
29 pages 32

K Monte Carlo EG1 EG2 EG3

2 40 55.4526 (104 ) 55.4551 55.4628 55.4627
26.08% 25.86% 26.08% 26.08%
60 40.0472 (104 ) 40.0386 40.047 40.0472
25.92% 25.88% 25.92% 25.92%
80 27.3445 (2.104 ) 27.3393 27.3444 27.3446
25.83% 25.82% 25.83% 25.83%
100 17.9604 (2.104 ) 17.9574 17.9603 17.9604
25.77% 25.77% 25.77% 25.77%
120 11.5237 (2.104 ) 11.5203 11.5236 11.5236
25.73% 25.72% 25.73% 25.73%
180 2.9074 (104 ) 2.9028 2.9074 2.9074
25.64% 25.63% 25.64% 25.64%
250 0.60324 (5.105 ) 0.6007 0.60328 0.60325
25.58% 25.56% 25.58% 25.58%
6 40 45.0211 (2.104 ) 44.9362 45.0188 45.0217
28.64% 27.50% 28.61% 28.64%
60 30.8275 (3.104 ) 30.7691 30.8228 30.8272
28.10% 27.89% 28.08% 28.10%
80 20.0422 (4.104 ) 20.0108 20.0393 20.0419
27.78% 27.72% 27.78% 27.78%
100 12.6282 (4.104 ) 12.5993 12.6258 12.628
27.57% 27.52% 27.57% 27.57%
120 7.8293 (4.104 ) 7.7927 7.8272 7.8292
27.42% 27.34% 27.41% 27.42%
180 1.8401 (2.104 ) 1.8005 1.8400 1.8401
27.12% 26.97% 27.12% 27.12%
250 0.36396 (2.104 ) 0.3449 0.36483 0.36399
26.91% 26.70% 26.92% 26.92%
10 40 35.0693 (3.104 ) 34.8603 35.046 35.0681
31.87% 30.33% 31.70% 31.86%
60 22.6743 (5.104 ) 22.5611 22.6572 22.672
30.81% 30.47% 30.76% 30.64%
80 14.0403 (7.104 ) 13.9583 14.0285 14.0388
29.12% 29.02% 29.12% 29.12%
100 8.5044 (8.104 ) 8.4067 8.4942 8.5036
29.79% 29.58% 29.77% 29.79%
120 5.1082 (6.104 ) 4.9913 5.1004 5.1082
29.49% 29.20% 29.47% 29.49%
180 1.1263 (4.104 ) 1.0274 1.1280 1.1268
28.92% 28.37% 28.93% 28.92%
250 0.21358 (2.104 ) 0.1707 0.2174 0.21368
28.53% 27.71% 28.60% 28.53%

Table 3: European Call option prices, with = 25%, r = 6%, q = 0%, S0 = 100, n = 3,
y 0.02.
29 pages 33

n K Monte Carlo EG1 EG2 EG3

3 40 50.2015 (104 ) 50.1659 50.2017 50.2017
27.28% 26.60% 27.29% 27.29%
60 35.325 (2.104 ) 35.2945 35.3236 35.325
26.95% 26.83% 26.95% 26.95%
80 23.5389 (2.104 ) 23.523 23.538 23.5389
26.76% 26.72% 26.76% 26.76%
100 15.1403 (3.104 ) 15.1283 15.1395 15.1402
26.63% 26.61% 26.63% 26.63%
120 9.5465 (3.104 ) 9.5317 9.5458 9.5464
26.53% 26.50% 26.53% 26.53%
180 2.3228 (2.104 ) 2.3045 2.3226 2.3227
26.34% 26.29% 26.34% 26.34%
250 0.4702 (104 ) 0.4609 0.4704 0.4702
26.22% 26.14% 26.22% 26.22%
5 40 45.2638 (3.104 ) 45.1382 45.2541 45.2636
29.06% 28.19% 29.00% 29.06%
60 33.6548 (4.104 ) 33.5814 33.6474 33.6542
28.45% 28.24% 28.43% 28.45%
80 24.6489 (4.104 ) 24.6068 24.6445 24.6485
28.09% 28.01% 28.08% 28.09%
100 17.9586 (5.104 ) 17.9226 17.9552 17.9583
27.85% 27.79% 27.84% 27.85%
120 13.0912 (5.104 ) 13.0496 13.0883 13.0912
27.67% 27.61% 27.66% 27.67%
180 5.2225 (4.104 ) 5.1651 5.2209 5.2225
27.33% 27.22% 27.33% 27.33%
250 1.9317 (3.104 ) 1.8848 1.9321 1.9319
27.10% 26.95% 27.10% 27.10%
10 40 36.474 (8.104 ) 36.1277 36.3948 36.456
34.64% 33.33% 34.34% 34.57%
60 29.9636 (8.104 ) 29.7413 29.9152 29.9536
33.06% 29.91% 32.94% 33.03%
80 24.7978 (9.104 ) 24.623 24.7633 24.7919
32.14% 31.81% 32.08% 32.13%
100 20.6818 (8.104 ) 20.514 20.6532 20.6782
31.53% 31.26% 31.49% 31.53%
120 17.3784 (8.104 ) 17.2 17.3528 17.3763
31.09% 30.82% 31.05% 31.09%
180 10.7324 (8.104 ) 10.5102 10.7136 10.7326
30.25% 29.92% 30.22% 30.25%
250 6.5113 (7.104 ) 6.2773 6.5018 6.5128
29.69% 29.30% 29.68% 29.70%

Table 4: European Call option, with = 25%, r = 6%, q = 0%, S0 = 100, y 0.02, 4.
29 pages 34

2 6 10
(y,)
P(ST < 0) 0 2.1019 9.106
(y,)
Table 5: P(ST < 0), with = 25%, r = 6%, q = 0%, S0 = 100, n = 3 and y 0.02.

n 3 5 10
(y,)
P(ST < 0) 0 6.106 0.03
(y,)
Table 6: P(ST < 0), with = 25%, r = 6%, q = 0%, S0 = 100, 4 and y 0.02.