VIBRATO AND AUTOMATIC DIFFERENTIATION FOR HIGH
ORDER DERIVATIVES AND SENSITIVITIES OF FINANCIAL
OPTIONS

arXiv:1606.06143v1 [q-fin.CP] 20 Jun 2016

GILLES PAGÈS

∗,

OLIVIER PIRONNEAU

† , AND

GUILLAUME SALL

‡

Abstract. This paper deals with the computation of second or higher order greeks of financial
securities. It combines two methods, Vibrato and automatic differentiation and compares with other
methods. We show that this combined technique is faster than standard finite difference, more
stable than automatic differentiation of second order derivatives and more general than Malliavin
Calculus. We present a generic framework to compute any greeks and present several applications
on different types of financial contracts: European and American options, multidimensional Basket
Call and stochastic volatility models such as Heston’s model. We give also an algorithm to compute
derivatives for the Longstaff-Schwartz Monte Carlo method for American options. We also extend
automatic differentiation for second order derivatives of options with non-twice differentiable payoff.
Key words. Financial securities, risk assessment, greeks, Monte-Carlo, automatic differentiation, vibrato.
AMS subject classifications. 37M25, 65N99

1. Introduction. Due to BASEL III regulations, banks are requested to evaluate the sensitivities of their portfolios every day (risk assessment). Some of these
portfolios are huge and sensitivities are time consuming to compute accurately. Faced
with the problem of building a software for this task and distrusting automatic differentiation for non-differentiable functions, we turned to an idea developed by Mike
Giles called Vibrato.
Vibrato at core is a differentiation of a combination of likelihood ratio method
and pathwise evaluation. In Giles [12], [13], it is shown that the computing time,
stability and precision are enhanced compared with numerical differentiation of the
full Monte Carlo path.
In many cases, double sensitivities, i.e. second derivatives with respect to parameters, are needed (e.g. gamma hedging).
Finite difference approximation of sensitivities is a very simple method but its
precision is hard to control because it relies on the appropriate choice of the increment. Automatic differentiation of computer programs bypass the difficulty and its
computing cost is similar to finite difference, if not cheaper. But in finance the payoff
is never twice differentiable and so generalized derivatives have to be used requiring
approximations of Dirac functions of which the precision is also doubtful.
The purpose of this paper is to investigate the feasibility of Vibrato for second
and higher derivatives. We will first compare Vibrato applied twice with the analytic
differentiation of Vibrato and show that it is equivalent; as the second is easier we
propose the best compromise for second derivatives: Automatic Differentiation of
Vibrato.
In [8], Capriotti has recently investigated the coupling of different mathematical
methods – namely pathwise and likelihood ratio methods – with an Automatic differ∗ Laboratoire de Probabilités et Modèles Aléatoires, UMR 7599, UPMC, Case 188, 4 pl. de Jussieu,
F-75252 Paris Cedex 5, France, gilles.pages@upmc.fr.
† Laboratoire Jacques Louis Lions, UMR 7598, Case 187, 4 pl. de Jussieu, F-75252 Paris Cedex
5, France, olivier.pironneau@upmc.fr.
‡ Laboratoire de Probabilités et Modèles Aléatoires, UMR 7599, UPMC, Case 188, 4 pl. de Jussieu,
F-75252 Paris Cedex 5, France, guillaume.sall@upmc.fr.

1

2
entiation technique for the computation of the second order greeks; here we follow the
same idea but with Vibrato and also for the computation of higher order derivatives.
Automatic Differentiation (AD) of computer program as described by Greiwank
in [19], [20], Naumann in [33] and Hascoet in [22] can be used in direct or reverse mode.
In direct mode the computing cost is similar to finite difference but with no roundoff
errors on the results: the method is exact because every line of the computer program
which implements the financial option is differentiated exactly. The computing cost
of a first derivative is similar to running the program twice.
Unfortunately, for many financial products the first or the second sensitivities do
not exist at some point, such is the case for the standard Digital option at x = K;
even the payofff of the a plain vanilla European option is not twice differentiatble at
x = K, yet the Gamma is well defined due to the regularizing effect of the Brownian
motion (or the heat kernel) which gives sense to the expectation of a Dirac as a
pointwise value of a probability density; in short the end result is well defined but the
intermediate steps of AD are not.
We tested ADOL-C [21] and tried to compute the Hessian matrix for a standard
European Call option in the Black-Scholes model but the results were wrong. So we
adapted our AD library based on operator overloading by including approximations
of Dirac functions and obtained decent results; this is the second conclusion of the
paper: AD for second sensitivities can be made to work; it is simpler than Vibrato+AD
(VAD) but it is risky and slightly more computer intensive.
More details on AD can be found in Giles et al. [11], Pironneau [35], Capriotti
[7], Homescu [26] and the references therein.
An important constraint when designing costly software for risk assessment is to
be compatible with the history of the company which contracts the software; most of
the time, this rules out the use of partial differential equations (see [1]) as most quant
companies use Monte Carlo algorithms for pricing their portfolios.
For security derivatives computed by a Monte Carlo method, the computation
of their sensitivities with respect to a parameter is most easily approximated by
finite difference (also known as the shock method) thus requiring the reevaluation
of the security with an incremented parameter. There are two problems with this
method: it is imprecise when generalized to higher order derivatives and expensive
for multidimensional problems with multiple parameters. The nth derivative of a
security with p parameters requires (n + 1)p evaluations; furthermore the choice of
the perturbation parameter is tricky.
From a semi-analytical standpoint the most natural way to compute a sensitivity
is the pathwise method described in Glasserman [15] which amounts to compute
the derivative of the payoff for each simulation path. Unfortunately, this technique
happens to be inefficient for certain types of payoffs including some often used in
quantitative finance like Digitals or Barrier options. For instance, as it is not possible
to obtain the Delta of a Digital Call that way (the derivative of the expectation of a
Digital payoff is not equal to the expectation of the derivative of the Digital payoff,
which in fact does not exist as a function), the pathwise method cannot evaluate the
Gamma of a Call option in a standard Black-Scholes model. The pathwise derivative
estimation is also called infinitesimal perturbation and there is a extensive literature
on this subject; see for example Ho et al. [24], in Suri et al. [39] and in L’Ecuyer [28].
A general framework for some applications to option pricing is given in Glasserman
[14].
There are also two well known mathematical methods to obtain sensibilities, the

3
so-called log-likelihood ratio method and the Malliavin calculus. However, like the
pathwise method, both have their own advantage and drawback. For the former,
the method consists in differentiating the probability density of the underlying and
clearly, it is not possible to compute greeks if the probability density of the underlying
is not known. Yet, the method has a great advantage in that the probability densities
are generally smooth functions of their parameters, even when payoff functions are
not. This method has been developed primarily in Glynn [17], Reiman et al. [36],
Rubinstein [37] and some financial applications in Broadie et al. [5] and Glasserman
et al. [16].
As for the Malliavin calculus, the computation of the greeks consists in writing
the expectation of the orignal payoff function times a specific factor i.e. the Malliavin
weight which is a Skorohod integral, the adjoint operator of the Malliavin derivative.
The main problem of this method is that the computation of the Malliavin weight can
be complex and/or computationally costly for a high dimensional problem. Several
articles deal with the computation of greeks via Malliavin calculus, Fournié et al.
[10], Benhamou [2] and Gobet et al. [18] to cite a few. The precision of the Malliavin
formulae also degenerates for short maturities, especially for the ∆-hedge.
Both the likelihood ratio and the Malliavin calculus are generally faster than the
pathwise or finite difference method because, once the terms in front of the payoff
function (the weight is computed analytically), the approximation of a greek in a
one-dimensional case is almost equivalent to the cost of the evaluation of the pricing
function. One systematic drawback is the implementation of these method in the
financial industry is limited by the specific analysis required by each new payoff.
The paper is organized as follows; in section 2 we begin by recalling the Vibrato
method for first order derivatives as in Giles [12] for the univariate and the multivariate
case. We then generalize the method for the second and higher order derivatives with
respect to one or several parameters and we describe the coupling to an analytical or
Automatic differentiation method to obtain an additional order of differentiation.
In section 3, we recall briefly the different methods of Automatic differentiation.
We describe the direct and the adjoint or reverse mode to differentiate a computer
program. We also explain how it can be extended to some non differentiable functions.
Section 4 deals with several applications to different derivative securities. We
show some results of second order derivatives (Gamma and Vanna) and third order
derivatives in the case of a standard European Call option: the sensitivity of the
Gamma with respect to changes in the underlying asset and a cross-derivatives with
respect to the underlying asset, the volatility and the interest rate. Also, we compare
different technique of Automatic differentiation and we give some details about our
computer implementations.
In section 5 we study some path-dependent products; we apply the combined
Vibrato plus Automatic differentiation method to the computation of the Gamma for
an American Put option computed with the Longstaff Schwartz algorithm [31]. We
also illustrate the method on a multidimensional Basket option (section 4) and on a
European Call with Heston’s model in section 6. In section 7, we study the computing
time for the evaluation of the Hessian matrix of a standard European Call Option in
the Black-Scholes model. Finally, in section 8 we compare VADs to Malliavin’s and
to the likelihood ratio method in the context of short maturities.
2. Vibrato. Vibrato was introduced by Giles in [12]; it is based on a reformulation of the payoff which is better suited to differentiation. The Monte Carlo path is

4
split into the last time step and its past. Let us explain the method on a plain vanilla
multi-dimensional option.
First, let us recall the likelihood ratio method for derivatives.
Let the parameter set Θ be a subset of Rp . Let b : Θ×Rd → Rd , σ : Θ×Rd → Rd×q be
continuous functions, locally Lipstchitz in the space variable, with linear growth, both
uniformly in θ ∈ Θ. We omit time as variable in both b and σ only for simplicity. And
let (Wt )t≥0 be a q-dimensional standard Brownian motion defined on a probability
space (Ω, A, P).
Lemma 2.1. (Log-likelihood ratio)
Let p(θ, ·) be the probability density of a random variable X(θ), which is function
of θ; consider
Z
V (y)p(θ, y)dy.
(2.1)
E[V (X(θ))] =
Rd

If θ 7→ p(θ, ·) is differentiable at θ0 ∈ Θ for all y, then, under a standard domination or
a uniform integrability assumption one can interchange differentiation and integration
: for i = 1, .., p,
i
∂ h
E[V (X(θ))] 0 =
∂θi
|θ



∂ log p
∂ log p 0
0
(θ , y)p(θ , y)dy = E V (X(θ))
(θ, X(θ))
V (y)
.
∂θi
∂θi
Rd
|θ 0
(2.2)

Z

2.1. Vibrato for a European Contract. Let X = (Xt )t∈[0,T ] be a diffusion
process, the strong solution of the following Stochastic Differential Equation (SDE)
dXt = b (θ, Xt ) dt + σ(θ, Xt )dWt , X0 = x.

(2.3)

For simplicity and without loss of generality, we assume that q = d; so σ is a square
matrix. Obviously, Xt depends on θ; for clarity, we write Xt (θ) when the context
requires it.
Given an integer n > 0, the Euler scheme with constant step h = Tn , defined below in (2.3), approximates Xt at time tnk = kh , i.e. X̄kn ≈ Xkh , and it is recursively
defined by
√
n
n
n
X̄kn = X̄k−1
+ b(θ, X̄k−1
)h + σ(θ, X̄k−1
) hZk , X̄0n = x, k = 1, . . . , n,
(2.4)
where {Zk }k=1,..,n are independent random Gaussian N (0, Id ) vectors. The relation
between W and Z is
√
Wtnk − Wtnk−1 = hZk .
(2.5)
√
Note that X̄nn = µn−1 (θ) + σn−1 (θ) hZn with
√
n
n
n
(2.6)
µn−1 (θ) = X̄n−1
(θ) + b(θ, X̄n−1
(θ))h and σn−1 (θ) = σ(θ, X̄n−1
(θ)) h.
Then, for any Borel function V : Rd → R such that E|V (X̄nn (θ))| < +∞,


 


 
n
.
E V (X̄nn (θ)) = E E V (X̄nn (θ)) | (Wtnk )k=0,...,n−1 = E E V (X̄nn (θ)) | X̄n−1
(2.7)

5
This follows from the obvious fact that the Euler scheme defines a Markov chain X̄
with respect to the filtration Fk = σ(Wtnℓ , ℓ = 0, . . . , k).
Furthermore, by homogeneity of the chain,
n
o
√
  


n
= Ex V (X̄1n (x, θ)) |x=X̄ n = E[V (µ + σ hZ)] µ = µn−1 (θ) .
E V (X̄nn (θ)) | X̄n−1
n−1

σ = σn−1 (θ)

(2.8)
Where X̄1n (x, θ) denotes the value at time tn1 of the Euler scheme with k = 1, starting
at x and where the last expectation is with respect to Z.
h
i
√
2.2. First Order Vibrato. We denote ϕ(µ, σ) = E V (µ + σ hZ) . From

(2.7) and (2.8), for any i ∈ (1, . . . , p)
#
"


o
√
∂
∂ϕ
∂ n
n
E[V (X̄n (θ))] = E
(µn−1 (θ), σn−1 (θ))
E[V (µ + σ hZ)] µ = µn−1 (θ) = E
∂θi
∂θi
∂θi
σ = σn−1 (θ)
(2.9)
and
∂µn−1 ∂ϕ
∂σn−1 ∂ϕ
∂ϕ
(µn−1 , σn−1 ) =
·
:
(µn−1 , σn−1 ) +
(µn−1 , σn−1 ) (2.10)
∂θi
∂θi
∂µ
∂θi
∂σ
where · denotes the scalar product and : denotes the trace of the product of the
matrices.
Lemma 2.2.
∂Xt
The θi -tangent process to X, Yt =
, is defined as the solution of the following
∂θi
SDE (see Kunita[27] for a proof )




∂X0
dYt = b′θi (θ, Xt ) + b′x (θ, Xt )Yt dt + σθ′ i (θ, Xt ) + σx′ (θ, Xt )Yt dWt , Y0 =
∂θi
(2.11)
where the primes denote standard derivatives. As for X̄kn in (2.3), we may discretize
(2.11) by
√



n
Ȳk+1
= Ȳkn + b′θi (θ, X̄kn ) + b′x (θ, X̄kn )Ȳkn h + σθ′ i (θ, X̄kn ) + σx′ (θ, X̄kn )Ȳkn hZ(2.12)
k+1 .
Then from (2.6),



∂µn−1
n
n
n
n
(θ)) + b′x (θ, X̄n−1
(θ))Ȳn−1
(θ)
= Ȳn−1
(θ) + h b′θi (θ, X̄n−1
∂θi
√ 

∂σn−1
n
n
n
(θ)) + σx′ (θ, X̄n−1
(θ))Ȳn−1
(θ) .
= h σθ′ i (θ, X̄n−1
∂θi

So far we have shown the following lemma.

Lemma 2.3. When Xnn (θ) is given by (2.3), then

(2.13)

∂
E[V (X̄nn (θ))] is given by
∂θi

(2.9) with (2.10), (2.13) and (2.12).
In (2.3) b and σ are constant in the time interval (kh, (k + 1)h), therefore the
n
conditional probability of X̄nn given X̄n−1
given by
T −1
1
1
p e− 2 (x−µ) Σ (x−µ)
p(x) = √
d
( 2π) |Σ|

(2.14)

6
where µ and Σ = hσσ T are evaluated at time (n − 1)h and given by (2.6). As in
Dwyer et al. [9],
∂
∂
1
1
log p(x) = Σ−1 (x − µ),
log p(x) = − Σ−1 + Σ−1 (x − µ)(x − µ)T Σ−1 ⇒
∂µ
∂Σ
2
2
Z
1
∂
∂
log p(x)|x=Xnn = σ −T √ ,
log p(x)|x=Xnn =
σ −T (ZZ T − I)σ −1 .
∂µ
2h
h ∂Σ
Finally, applying Lemma 2.3 and Lemma 2.1 yields the following proposition
Theorem 2.4. (Vibrato, multidimensional first order case)
"
#
o
√
∂ n
∂
n
E[V (µ + σ hZ)] µ = µn−1 (θ)
E[V (X̄n (θ))] = E
∂θi
∂θi
σ = σn−1 (θ)

i
h
√
1 ∂µ
· E V (µ + σ hZ)σ −T Z µ = µn−1 (θ)
=E √
σ = σn−1 (θ)
h ∂θi
+

i
h
√
1 ∂Σ
: E V (µ + σ hZ)σ −T (ZZ T − I)σ −1
2h ∂θi

µ = µn−1 (θ)
σ = σn−1 (θ)

.



(2.15)

2.3. Antithetic Vibrato. One can expect to improve the above formula – that
is, reducing its variance – by the means of antithetic transform (see section 2.6 below
for a short discussion) The following holds:
i 1 h

i
h
√
√
√
E V (µ + σ hZ)σ −T Z = E V (µ + σ hZ) − V (µ − σ hZ) σ −T Z .
2

(2.16)

similarly, using E[ZZ T − I] = 0,
h
i
√
E V (µ + σ hZ)σ −T (ZZ T − I)σ −1

i
√
√
1 h
= E V (µ + σ hZ) − 2V (µ) + V (µ − σ hZ) σ −T (ZZ T − I)σ −1 . (2.17)
2
Corollary 2.5. (One dimensional case, d=1)



 Z 
√
√
1
∂µ
∂
n
E[V (X̄n (θ))] = E
E V (µ + σ hZ) − V (µ − σ hZ) √
µ = µn−1 (θ)
∂θi
2
∂θi
σ = σn−1 (θ)
σ h




2
√
√
∂σ
Z −1
 (2.18)
√
E V (µ + σ hZ) − 2V (µ) + V (µ − σ hZ)
+
µ = µn−1 (θ)
∂θi
σ h
σ = σn−1 (θ)

Conceptual Algorithm. In figure 1 we have illustrated the Vibrato decomposition
at the path level. To implement the above one must perform the following steps:
1. Choose the number of time step n, the number of Monte-Carlo path M for
the n − 1 first time steps, the number MZ of replication variable Z for the
last time step.
2. For each Monte-Carlo path j = 1..M
• Compute {Xkn }k=1:n−1 , µn−1 , σn−1 by (2.3), (2.6).

7
1.5
1.4
1.3
1.2
1.1
1
0.9
0.8
0.7
0.6
First part of the path (n-1 th terms of the Euler scheme)

Last term

0.5
0

0.2

0.4

0.6

0.8

1

Figure 1. Scheme of simulation path of the Vibrato decomposition.

• Compute V (µn−1 )
∂σn−1
∂µn−1
and
by (2.11), (2.13) and (2.12)
• Compute
∂θi
∂θi
• Replicate MZ times the last time step, i.e.
For mZ ∈ (1, . . . , MZ )
√
√
– Compute V (µn−1 + σn−1 hZ (mZ ) ) and V (µn−1 − σn−1 hZ (mZ ) )
3. In (2.18) compute the inner expected value by averaging over all MZ results,
∂µ
∂σ
and ∂θ
and then average over the M paths.
then multiply by ∂θ
i
i
Remark 1. For simple cases such as of the sensibilities of European options, a
small MZ suffices; this is because there is another average with respect to M in the
outer loop.
Remark 2. For European options one may also use the Black-Scholes formula
for the expected value in (2.15).
2.4. Second Derivatives. Assume that X0 , b and σ depend on two parameters
(θ1 , θ2 ) ∈ Θ2 . There are two ways to compute second order derivatives. Either by
differentiating the Vibrato (2.15) while using Lemma 2.1 or by applying the Vibrato
idea to the second derivative.
2.4.1. Second Derivatives by Differentiation of Vibrato. Let us differentiate (2.15) with respect to a second parameter θj :

h
i
√
1  ∂2µ
∂2
E[V (XT )] = E √
· E V (µ + σ hZ)σ −T Z
∂θi ∂θj
h ∂θi ∂θj
i
√
∂µ ∂ h
·
E V (µ + σ hZ)σ −T Z
+
µ = µn−1 (θ)
∂θi ∂θj
σ = σn−1 (θ)
i

h
2
√
∂ Σ
1
: E V (µ + σ hZ)σ −T (ZZ T − I)σ −1
+
2h ∂θi ∂θj
i
√
∂Σ
∂ h
+
:
E V (µ + σ hZ)σ −T (ZZ T − I)σ −1
∂θi ∂θj

µ = µn−1 (θ)
σ = σn−1 (θ)




(2.19)
The derivatives can be expanded further; for instance in the one dimensional case and
after a tedious algebra one obtains:

8
Theorem 2.6. (Second Order by Differentiation of Vibrato)
"


  2 
√
√
∂2µ
∂2
∂µ
Z
Z2 − 1
√
E[V
(X
)]
=
E
E
V
(µ
+
σ
E
V
(µ
+
σ
+
hZ)
hZ)
T
∂θ2
∂θ2
∂θ
σ2h
σ h
 2 

4
2
√
∂σ
Z − 5Z + 2
+
E V (µ + σ hZ)
2
∂θ 


σ h
√
√
∂2σ
Z2 − 1
Z 3 − 3Z
∂µ ∂σ
+ 2 E V (µ + σ hZ) √
E V (µ + σ hZ)
(2.20)
+2
∂θ
∂θ ∂θ
σ2h
σ h

2.4.2. Second Derivatives by Second Order Vibrato. The same Vibrato
strategy can be applied also directly to second derivatives.
As before the derivatives are transfered to the PDF p of XT :
Z
Z
V (x) ∂ 2 p
∂ 2 ln p
∂ ln p ∂ ln p
∂2
V (x)[
E[V (XT )]=
p(x)dx =
+
]p(x)dx
∂θi ∂θj
∂θi ∂θj
Rd p(x)∂θi ∂θj
Rd  ∂θi ∂θj
∂ 2 ln p
∂ ln p ∂ ln p
= E V (x)
(2.21)
+
∂θi ∂θj
∂θi ∂θj
Then
∂2ϕ
∂2
E[V (X̄Tn (θ1 , θ2 ))] =
(µ, σ)
∂θ1 ∂θ2
∂θ1 ∂θ2
∂σ ∂σ ∂ 2 ϕ
∂ 2 µ ∂ϕ
∂µ ∂µ ∂ 2 ϕ
(µ,
σ)
+
(µ,
σ)
+
(µ, σ)
=
∂θ1 ∂θ2 ∂µ2
∂θ1 ∂θ2 ∂σ 2
∂θ1 ∂θ2 ∂µ


2
2
∂σ ∂µ
∂ σ ∂ϕ
∂ ϕ
∂µ ∂σ
+
(µ, σ) +
(µ, σ).
+
∂θ1 ∂θ2 ∂σ
∂θ1 ∂θ2
∂θ1 ∂θ2 ∂µ∂σ
∂2
∂2
µn−1 (θ1 , θ2 ) and
σn−1 (θ1 , θ2 ).
∂θ1 ∂θ2
∂θ1 ∂θ2
It requires the computation of the first derivative with respect to θi of the tangent
(2)
process Yt , that we denote Yt (θ1 , θ2 ).
Then (2.13) is differentiated and an elementary though tedious computations
yields the following proposition:
Proposition 2.7.
(i)
The θi -tangent process Yt defined above in Lemma 2.11 has a θj -tangent process
(ij)
Yt
defined by
h
(j)
(i)
(ij)
dYt
= b′′θi θj (θ1 , θ2 , Xt ) + b′′θi ,x (θ1 , θ2 , Xt )Yt + b′′θj ,x (θ1 , θ2 , Xt )Yt
i
(i) (j)
(ij)
dt
+b′′x2 (θ1 , θ2 , Xt )Yt Yt + b′x (θ1 , θ2 , Xt )Yt
h
(j)
(i)
+ σθ′′i θj (θ1 , θ2 , Xt ) + σθ′′i ,x (θ1 , θ2 , Xt )Yt + σθ′′j ,x (θ1 , θ2 , Xt )Yt
i
(ij)
(i) (j)
dWt .
+σx′′2 (θ1 , θ2 , Xt )Yt Yt + σx′ (θ1 , θ2 , Xt )Yt

We need to calculate the two new terms

Finally in the univariate case θ = θ1 = θ2 this gives
Proposition 2.8. (Second Order Vibrato)
∂2
E[V (XT )] =
∂θ2

9
"


  2 h
  2 
√
√
√
∂2µ
Z
Z2 − 1
∂σ
∂µ
√
E
E V (µ + σ hZ)
+
E V (µ + σ hZ) 2
E V (µ + σ hZ)
+
2
∂θ
∂θ
σ h
∂θ
σ h





√
√
∂2σ
Z2 − 1
Z 3 − 3Z
Z 4 − 5Z 2 + 2
∂µ ∂σ
√
+
E
V
(µ
+
σ
hZ)
E
V
(µ
+
σ
hZ)
(2.22)
+
2
σ2h
∂θ2
∂θ ∂θ
σ2 h
σ h

Remark 3. It is equivalent to Proposition 2.6 hence to the direct differentiation
of Vibrato.
2.5. Higher Order Vibrato. The Vibrato-AD method can be generalized to
higher order of differentiation of Vibrato with respect to the parameter θ with the
help of the Faà di Bruno formula and its generalization to a composite function with
a vector argument, as given in Mishkov [32].
2.6. Antithetic Transform, Regularity and Variance. In this section, we
assume d = q = 1 for simplicity.
√
Starting from Vibrato ϕ(µ, σ) = E[f (µ + σ hZ)] and assuming f Lipschitz continuous with Lipschitz coefficients [f ]Lip , we have




√
√  Z
√
∂ϕ
Z
√
√
=
E
f
(µ
+
σ
. (2.23)
hZ) − f (µ − σZ h)
(µ, σ) = E f (µ + σ hZ)
∂µ
σ h
2σ h

Therefore the variance satisfies
Var




√
f (µ + σ hZ) − f (µ − σ hZ)
√

"
#


 Z 2
√
√
Z
√
√
f (µ + σ hZ) − f (µ − σ hZ)
≤E
2σ h
2σ h
"
#
√
(2σ hZ)2 2
≤ [f ]2Lip E
Z = [f ]2Lip E[Z 4 ] = 3[f ]2Lip .
4σ 2 h

As E[Z] = 0, we also have

 Z 
√
∂ϕ
√ .
(µ, σ) = E f (µ + σ hZ) − f (µ)
∂µ
σ h

(2.24)

Then,
Var



#
"

 Z 
 Z 2
√
√
√
√
f (µ + σ hZ) − f (µ)
≤E
f (µ + σ hZ) − f (µ)
σ h
σ h
h √
i
1
2
2 2
≤ 2 [f ]Lip E (σ hZ) Z = [f ]2Lip E[Z 4 ] = 3[f ]2Lip
σ h

Remark 4. The variances of formulae (2.23) and (2.24) are equivalent but the
latter is less expensive to compute.
If f is differentiable and f ′ has polynomial
growth, we also have
√
∂ϕ
(µ, σ) = E[f ′ (µ + σ hZ)].
∂µ

(2.25)

Thus,

h
i
2 
√
√
Var f ′ (µ + σ hZ) ≤ E f ′ (µ + σ hZ)
≤kf ′ k2∞ .

Remark 5. Let f ]Lip denote the Lipschitz constant of f . If f ′ is bounded, we
have [f ]Lip = kf ′ k∞ then the expression in (2.25) has a smaller variance than (2.23)

10
and (2.24). Assume that f ′ is Lipschitz continuous with Lipschitz coefficients [f ′ ]Lip .
We can improve the efficiency of (2.25) because
h
i
h
i
√
√
Var f ′ (µ + σ hZ) = Var f ′ (µ + σ hZ) − f ′ (µ)


2
√
′
′
≤ [f ′ ]2Lip hσ 2 E[Z 2 ] ≤ [f ′ ]Lip hσ 2
≤ E f (µ + σ hZ) − f (µ)
Remark 6. Assuming that f (x) = 1{x≤K} , clearly we cannot differentiate inside
the expectation and the estimation of the variance seen previously can not be applied.
2.6.1. Indicator Function. Let us assume that f (x) = 1{x≤K} . To simplify
assume that K ≤ µ, we have
√
√
h
o − 1n
o = 1n
io ,
f (µ + σ hZ) − f (µ − σ hZ) = 1nZ≤ K−µ
µ−K
√
√
√ , √
Z∈
/ K−µ
Z≥ µ−K
σ

h

σ

h

σ

h

σ

h

hence
 Z

√
√
1
io .
f (µ + σ hZ) − f (µ − σ hZ) √ = √ |Z|1nZ ∈/ h K−µ
µ−K
√ , √
σ h σ h
σ h
σ h
For the variance, we have
√

 Z 
√
√
Var f (µ + σ hZ) − f (µ − σ hZ)
σ h #
"
2


√
√
Z
√
.
≤E
f (µ + σ hZ) − f (µ − σ hZ)
σ h


By Cauchy-Schwarz we can write
#
"
 Z 2

√
√
√
E
f (µ + σ hZ) − f (µ − σ hZ)
σ h




√
√
2
1
1
2
2 n
io
h
=
hZ)
−
f
(µ
−
σ
hZ)
=
E
Z
E
Z
1
f
(µ
+
σ
√ , µ−K
√
Z∈
/ K−µ
2σ 2 h
2σ 2 h
σ h σ h
√  
 1
 1
 

2
2
1
1
µ−K
3
K −µ µ−K
4 2
√ , √
≤
2P Z ≥ √
E[Z ]
≤
.
P Z∈
/
2h
2σ 2 h
2σ
σ h σ h
σ h

Then

√

3
2σ 2 h



√
 1

2
µ−K
6
=
2P Z ≥ √
2h
2σ
σ h

Z

+∞
µ−K
√
σ h

2

e

− u2

du
√
2π

!1

2

.

a2

e− 2
, so when a → +∞,
Now, ∀ a > 0, P(Z ≥ a) ≤ √
a 2π
Var

≤



r
(µ−K)2
−
 Z 
√
1
3 e 4σ2 h
q
√
≤ 2
f (µ + σ hZ) − f (µ − σ hZ)
σ h 2 (2π) 41 µ−K
σ h
√
√

σ

1

1
4

3
2

(2π) σ h

3
4

r

−

(µ−K)2
4σ2 h

3e
√
−→
2 µ − K σ→0



0
+∞

h

if µ 6= K
otherwise.

The fact that such estimate can be obtained with non differentiable f demonstrates
the power of the Vibrato technique.

11
3. Second Derivatives by Vibrato plus Automatic Differentiation (VAD).
The differentiation that leads to formula (2.22) can be derived automatically by AD;
then one has just to write a computer program that implements the formula of proposition 2.19 and apply automatic differentiation to the computer program. We recall
here the basis of AD.
3.1. Automatic Differentiation. Consider a function z = f (u) implemented
in C or C++ by
double f(double u){...}
To find an approximation of zu′ , one could call in C
double dxdu= (f(u + du)-f(u))/du
because
zu′ = f ′ (u) =

f (u + du) − f (u)
+ O(|du|).
du

A good precision ought to be reached by choosing du small. However arithmetic truncation limits the accuracy and shows that it is not easy to choose du appropriately
because beyond a certain threshold, the accuracy of the finite difference formula degenerates due to an almost zero over almost zero ratio. As described in Squire et al.

Figure 2. Precision (log-log plot of |dzdu −
cos(1.)| computed with the forward finite difference formula to evaluate sin′ (u) at u = 1.

Figure 3. Same as Fig. 2 but with the finite
difference which uses complex increments; both
test have been done with Maple-14

[38], one simple remedy is to use complex imaginary increments because
Re

′′′
f (u + idu) − f (u)
f (u + idu)
du2
= Re
= f ′ (u) − Ref (u + iθdu)
idu
idu
6

leads to f ′ (u) = Re[f (u + idu)/(idu)] where the numerator is no longer the result of
a difference of two terms. Indeed tests show that the error does not detoriate when
du → 0 (figure 3). Hence one can choose du = 10−8 to render the last term with a
O(10−16 ) accuracy thus obtaining an essentially exact result.
The cost for using this formula is two evaluations of f (), and the programming
requires to redefine all double as std::complex of the Standard Template Library in
C++.

12
3.2. AD in Direct Mode. A conceptually better idea is based on the fact
that each line of a computer program is differentiable except at switching points of
branching statements like if and at zeros of the sqrt functions etc.
Denoting by dx the differential of a variable x, the differential of a*b is da*b+a*db,
the differential of sin(x) is cos(x)dx, etc. . . By operator overloading, this algebra
can be built into a C++ class, called ddouble here:
class ddouble {
public: double val[2];
ddouble(double a=0, double b=0){ val[1]=b; val[0]=a; }
ddouble operator=(const ddouble& a)
{ val[1] = a.val[1]; val[0]=a.val[0]; return *this; }
ddouble operator - (const ddouble& a, const ddouble& b)
{ return ddouble(a.val[0] - b.val[0],a.val[1] - b.val[1]); }
ddouble operator * (const ddouble& a, const ddouble& b)
{ return ddouble(a.val[0] * b.val[0], a.val[1]*b.val[0]
+ a.val[0] * b.val[1]); }
... };
So all ddouble variables have a 2-array of data: val[0] contains the value of the variable and val[1] the value of its differential. Notice that the constructor of ddouble
assigns zero by default to val[1].
To understand how it works, consider the C++ example of figure 4 which calls a
function f (u, ud ) = (u−ud )2 for u = 2 and ud = 0.1. Figure 5 shows the same program
where double has been changed to ddouble and the initialization of u implies that
its differential is equal to 1. The printing statement displays now the differential of f
which is also its derivative with respect to u if all parameters have their differential
initialized to 0 except u for which has du = 1.
Writing the class double with all
double f(double u, double u_d)
{ double z = u-u_d;
return z*(u-u_d); }
int main() {
double u=2., u_d =0.1;
cout << f(u,u_d)<< endl;
return 0;
}

ddouble f(ddouble u, ddouble u_d)
{ ddouble z = u-u_d;
return z*(u-u_d); }
int main() {
ddouble u=ddouble(2.,1.), u_d = 0.1;
cout << f(u,u_d).val[1] << endl;
return 0;
}

Figure 4. A tiny C++ program to compute (u − ud )2 at u = 2, ud = 0.1.

Figure 5. The same program now comd
putes du
(u − ud )2 at u = 2, ud = 0.1.

functions and common arithmetic operators is a little tedious but not difficult. An
example can be downloaded from www.ann.jussieu.fr/pironneau.
The method can be extended to higher order derivatives easily. For second derivatives, for instance, a.val[4] will store a, its differentials with respected to the first and
second parameter, d1 a, d2 a and the second differential d12 a where the two parameters
can be the same. The second differential of a*b is a∗d12 b+d1a∗d2 b+d2 a∗d1 b+b∗d12a,
and so on.
df
Notice that du
can also be computed by the same program provided the first
d
line in the main() is replaced by ddouble u=2., u d=ddouble(0.1,1.);. However
df df
if both derivatives
are needed, then, either the program must be run twice
,
du dud
or the class ddouble must be modified to handle partial derivatives. In either case

13
the cost of computing n partial derivatives will be approximately n times that of the
original program; the reverse mode does not have this numerical complexity and must
be used when, say, n > 5 if expression templates with traits are used in the direct
mode and n > 5 otherwise [35].
3.3. AD in Reverse Mode. Consider finding Fθ′ where (u, θ) → F (u, θ) ∈ R
and u ∈ Rd and θ ∈ Rn . Assume that u is the solution of a well posed linear system
Au = Bθ + c.
The direct differentiation mode applied to the C++ program which implements
F will solve the linear system n times at the cost of d2 n operations at least.
The mathematical solution by calculus of variations starts with
Fθ′ dθ = (∂θ F )dθ + (∂u F )du with Adu = Bdθ,
then introduces p ∈ Rd solution of AT p = (∂u F )T and writes
(∂u F )du = (AT p)T du = pT Bdθ ⇒ Fθ′ dθ = (∂θ F + pT B)dθ.
The linear system for p is solved only once, i.e. performing O(d2 ) operations at least.
Thus, as the linear system is usually the costliest operation, this second method is
the most advantageous when n is large.
A C program only made of assignments can be seen as a triangular linear system
for the variables. Loops can be unrolled and seen as assignments and tests, etc. Then,
by the above method, the ith line of the program is multiplied by pi and p is computed
from the last line up; but the biggest difficulty is the book-keeping of the values of
the variables, at the time p is computed.
For instance, for the derivative of f=u+ud with respect to ud with u given by
{u=2*ud+4; u=3*u+ud;},u in the second line is not the same as u in the third line
and the program should be rewritten as u1=2*ud+4; u=3*u1+ud;. Then the system
for p is p2=1; p1=3*p2; and the derivative is 2*p1+p2+1=8.
In this study we have used the library adept 1.0 by R.J. Hogan described in
Hogan [25]. The nice part of this library is that the programming for the reverse
mode is quite similar to the direct mode presented above; all differentiable variables
have to be declared as ddouble and the variable with respect to which things are
differentiated is indicated at initialization, as above.
3.4. Non-Differentiable Functions. In finance, non-differentiability is everywhere. For instance, the second derivative Rin K of (x − K)+ does not exist at x = K
∞
as a function, yet the second derivative of 0 f (x)(x − K)+ dx is f (K). Distribution
theory extends the notion of derivative: the Heavyside function H(x) = 1{x≥0} has
the Dirac mass at zero δ(x) for derivative.
Automatic differentiation can be extended to handle this difficulty to some degree
by approximating the Dirac mass at 0 by the functions δ a (x) defined by
x2
1
δ a (x) = √ e− a .
aπ

Now, suppose f is discontinuous at x = z and smooth elsewhere; then
f (x) = f + (x)H(x − z) + f − (x)(1 − H(x − z))
hence
fz′ (x) = (f + )′z (x)H(x − z) + (f − )′z (x)(1 − H(x − z)) − (f + (z) − f − (z))δ(x − z)

14
Unless this last term is added, the computation of the second order sensitivities will
not be right.
If in the AD library the ramp function x+ is defined as xH(x) with its derivative
to be H(x), if H is defined with its derivative equal to δ a and if in the program which
computes the financial asset it is written that (x − K)+ = ramp(x − K), then the
second derivative in K computed by the AD library will be δ a (x − K). Moreover, it
will also compute
Z

0

∞

f (x)(x − K)+ dx ≈

N
1 X
f (ξi )δ a (ξi − K)
N i=1

where ξi are the N quadrature points of the integral or the Monte-Carlo points used
by the programmer to approximate the integral.
However, this trick does not solve all problems and one must be cautious; for
instance writing that (x − K)+ = (x − K)H(x − K) will not yield the right result.
Moreover, the precision is rather sensitive to the value of a.
Remark 7. Notice that finite difference (FD) is not plagued by this problem,
which means that FD with complex increment is quite a decent method for first order
sensitivities. For second order sensitivities the “very small over very small” problem
is still persistent.
4. VAD and the Black-Scholes Model. In this section, we implement and
test VAD and give a conceptual algorithm that describes the implementation of this
method (done automatically). We focus on indicators which depend on the solution
of an SDE, instead of the solution of the SDE itself. Let us take the example of a
standard European Call option in the Black-Scholes model.
4.1. Conceptual algorithm for VAD.
1. Generate M simulation paths with time step h = Tn of the underlying asset
∂X
X and its tangent process Y =
with respect to a parameter θ for k =
∂θ
0, . . . , n − 2:

√
∂X0

n
n
n
n
n
n

 X̄k+1 = X̄k + rhX̄k + X̄k σ hZk+1 , X̄0 = X0 , Ȳ0 =
∂θ


 √ 
√
∂
∂

n
n
n
n
n
n

σ h X̄k Zk+1 , .
(rh) X̄k + Ȳk σ h +
 Ȳk+1 = Ȳk + rhȲk +
∂θ
∂θ
(4.1)
2. For each simulation path
(a) Generate MZ last time steps (X̄T = X̄nn )
√
n
X̄nn = X̄n−1
(1 + rh + σ hZn ).

(4.2)

(b) Compute the first derivative with respect to θ by Vibrato using the
antithetic technique (formula (2.19) with σ(Xt ) equal Xt σ)
Zn
∂µn−1 1
∂VT
√
=
(VT+ − VT− ) n
∂θ
∂θ 2
X̄n−1 σ h
Z2 − 1
∂σn−1 1
(VT+ − 2VT• + VT− ) nn √ .
+
∂θ 2
X̄n−1 σ h

(4.3)

15
With VT±,• = (X̄T±,• − K)+ ,
(
√
n
n
n
+ rhX̄n−1
± σ X̄n−1
X̄T± = X̄n−1
hZn
n
n
X̄T• = X̄n−1 + rhX̄n−1 .

(4.4)

and
∂µn−1
∂
n
n
= Ȳn−1
(1 + rh) + X̄n−1
(rh)
∂θ
∂θ
√
√
∂
∂σn−1
n
n
= Ȳn−1
σ h + Xn−1
(σ h)
∂θ
∂θ

(4.5)

∂ −rT
(e
)VT to the result above.
∂θ
(c) Apply an Automatic Differentiation method on the computer program
that implements step 4.3 to compute the second derivative with respect
to θ at some θ∗ .
(d) Compute the mean per path i.e. over MZ .
3. Compute the mean of the resulting vector (over the M simulation paths) and
discount it.
If θ = T or θ = r, we have to add

4.2. Greeks. The Delta measures the rate of changes in the premium E[V (XT )]
with respect to changes in the spot price X0 .
The Gamma measures the rate of changes of the Delta with respect to changes
in the spot price. Gamma can be important for a Delta-hedging of a portfolio.
The Vanna is the second derivative of the premium with respect to σ and X0 . The
Vanna measures the rate of changes of the Delta with respect to changes in the volatility.
4.3. Numerical Test. For the generation of the random numbers, we chose the
standard Mersenne-Twister generator available in the version 11 of the C++ STL.
We take MZ = 1 i.e. we simulate only one last time step per path; for all the test
cases except for the European Call contract in the Black-Scholes model. However, for
the European Call in a Black-Scholes model, we used a multiple time steps with the
Euler scheme with or without a Brownian bridge.
The parameters considered in the following numerical experiments are K = 100,
σ = 20% and r = 5%, T = 1 year. The initial price of the risky asset price is varying
from 1 to 200. The Monte Carlo parameters are set to 100, 000 simulation paths, 25
time steps.
4.3.1. Preliminary Numerical Test. Here, we focus on the numerical precision of VAD on the Gamma of a standard European Call contract with constant
volatility and drift for which there is an analytical Black Scholes formula. Since
Vibrato of Vibrato is similar to Vibrato+AD (VAD) it is pointless to compare the
two.
Recall (Proposition 2.6 & 2.8) that it is equivalent to apply Vibrato to Vibrato
or to apply automatic differentiation to Vibrato. However, the computation times are
different and naturally double Vibrato is faster.
We compare the analytical solution to those obtained with VAD but now for each
new set of parameters, we reuse the same sample of the random variables.

16
On figure 6, the Gammas are compared at X0 = 120; true value of the Gamma is
Γ0 = 0.0075003. The convergence with respect to the number of paths is also displayed
for two values of MZ . The method shows a good precision and fast convergence when
the number of paths for the final time step is increased.
✥ ✥✄✁

✥ ✥✝

❱❆✞ ✟✠✡☛
❆☞✌✍✎✏✑✒✌✍ ✓♦✍✔✏✑♦☞
✥ ✥✥✆✁

✥ ✥✄

✥ ✥✥✆

✥ ✥✥☎✁

✥ ✥✂✁

✥ ✥✥☎

✥ ✥✥✄✁

✥ ✥✂

✥ ✥✥✄

✥ ✥✥✂✁

✥ ✥✥✁

✥ ✥✥✂

❱❆✡ ☛☞✌✝ ✍☛✎✮
✥ ✥✥✁✁

❱❆✡ ☛☞✌✞ ✍☛✎✮
❆✏✑✒✓✔✕✖✑✒ ✗♦✒✘✔✕♦✏

✥

✥ ✥✥✁

✥

✄✥

☎✥

✆✥

✝✥

✂✥✥

✂✄✥

✂☎✥

✂✆✥

✂✝✥

✄✥✥

✥

✁✥✥✥

✝✥✥✥✥

✝✁✥✥✥

✞✥✥✥✥

✞✁✥✥✥

✟✥✥✥✥

✟✁✥✥✥

✠✥✥✥✥

✠✁✥✥✥

✁✥✥✥✥

Figure 6. On the left the Gamma versus Price is displayed when computed by VAD; the
analytical exact Gamma is also displayed; both curves overlap. On the right, the convergence history
at one point X0 = 120 is displayed with respect to the number of Monte Carlo samples MW . This
is done for two values of MZ (the number of the final time step), MZ = 1 (low curve) and MZ = 2
(upper curve).

The L2 -error denoted by εL2 is defined by
εL 2

P
1 X i
(Γ̄ − Γ0 )2 .
=
P i=1

(4.6)

On figure 7, we compare the results with and without variance reduction on Vibrato
at the final time step i.e. antithetic variables. The convergence history against the
number of simulation paths is displayed. Results show that variance reduction is
efficient on that test case. The standard error against the number of simulation paths
is also displayed. It is clear that a reduction variance is needed. It requires almost
ten times the number of simulation paths without the reduction variance technique to
obtain the same precision. The Gamma is computed for the same set of parameters
as given above.
On figures 8 we display the Vanna of an European Call option, computed with
VAD. And again, the convergence with respect to the number of simulation paths is
accelerated by more sampling of the final time step. Note that the Vanna requires
double the number of time steps
4.3.2. Third Order Derivatives. For third order derivatives, we compute second derivatives by Vibrato of Vibrato 2.6 and differentiate by AD (VVAD). The
sensitivity of the Gamma with respect to changes in X0 is ∂ 3 V /∂X03 . The sensitivity of the Vanna with respect to changes in the interest rate is ∂ 3 V /∂X0 ∂σ∂r. The
parameters of the European Call are the same but the Monte Carlo path number is
1, 000, 000 and 50 time steps for the discretization. The results are displayed on figure
9. The convergence is slow; we could not eliminate the small difference between the
analytical solution and the approximation by increasing the number of paths.
4.3.3. Ramp Function and High Order Derivatives. As mentioned in Section 3.4, it is possible to handle the non-differentiability of the function (x − K)+

17
✂✞✂✥

✥ ✥✁

❱❆✠ ✡☛☞ ❱✌ ✍✎❈✏
❱❆✠ ✡✑✒✓ ❱✌ ✍✎❈✏

❱❆✟ ✠✡☛ ❱☞ ✌✍❈✎
❱❆✟ ✠✏✑✒ ❱☞ ✌✍❈✎
✥ ✥✆☎

❆✓✔✕✖✑✏✗✔✕ ✘☛✕♦✑✏☛✓

✂✞✂✂✥
✥ ✥✆✄

✂✞✂✂✂✥

✥ ✥✆✂

✥ ✥✆✁

✥ ✁✂✝
✥ ✥✆

✥ ✁✂✆

✥ ✥✥☎

✥ ✥✥✄

✥ ✁✂☎
✥ ✥✥✂

✥ ✁✂✄

✥ ✥✥✁
✥

✆✥✥✥✥

✁✥✥✥✥

✝✥✥✥✥

✂✥✥✥✥

✞✥✥✥✥

✥✂

✥✂✂

✥✂✂✂

✥✂✂✂✂

✥✂✂✂✂✂

✥ ✟✂✆

Figure 7. On the left the Gamma versus the number of simulation paths is displayed when
computed by VAD with and without the variance reduction method on Z, the straight line is the
analytical solution at one point X0 = 120; On the right, the standard error of the two methods
versus the number of simulation paths with and without variance reduction.
✂✄

✲✡✁✄

❱❆✞ ✟✠✡☛
❆☞✌✍✎✏✑✒✌✍ ✓♦✍✔✏✑♦☞

✲✡✁✂
✲
✲ ✁
✲ ✁✠

✁✂✄

✲ ✁✟
✲ ✁✞

✁

✲ ✁✝
✲ ✁✆

✲✁✂✄

✲ ✁☎
❱❆☛ ☞✌✍ ✎☞✏✮
❱❆☛ ☞✌✍✠ ✎☞✏✮
❆✑✒✓✔✕✖✗✒✓ ✘♦✓✙✕✖♦✑

✲ ✁✄

✲

✁

✥✁

☎✁

✆✁

✝✁

✁✁

✥✁

☎✁

✆✁

✝✁

✥✁✁

✲ ✁✂

✡

✝✡✡✡

✡✡✡✡

✝✡✡✡

✠✡✡✡✡

✠✝✡✡✡

✟✡✡✡✡

✟✝✡✡✡

✞✡✡✡✡

✞✝✡✡✡

✝✡✡✡✡

Figure 8. On the left the Vanna versus Price is displayed when computed by VAD; the analytical
exact Vanna is also displayed; both curves overlap. On the right, the convergence history at one
point X0 = 120 is displayed with respect to the number of Monte Carlo samples MW . This is done
for two values of MZ , MZ = 1 (lower curve) and MZ = 2 (upper curve).
✁ ✥

☎

✁

✂

✆

✁

✄

✁

✁

☎

✂

✁

✆

✲✁

✲ ✁

✆

✲✆
✲☎

✲ ✁

☎

✲✄

✲ ✁

✄

✲ ✂

✲ ✁

✂

❱❱❆✝ ✞✟✠✡
❆☛☞✌✍✎✏✑☞✌ ✒♦✌✓✎✏♦☛
✆

☎

✄

✂

✥

✥✆

✥☎

✥✄

✥✂

✆

✲ ✁

❱❱❆✝ ✞✟✠✡
❆☛☞✌✍✎✏✑☞✌ ✒♦✌✓✎✏♦☛
✂

✁✂

✆✂

☎✂

✄✂

✂✂

✁✂

✆✂

☎✂

✄✂

✁✂✂

Figure 9. On the left ∂ 3 V /∂X03 versus Price is displayed when computed by VVAD; the analytical exact curve is also displayed; both curves practically overlap. On the right, the same for the
Vanna with respect to changes in interest rate (∂ 3 V /∂X0 ∂σ∂r).

18
at x = K by using distribution theory and program the ramp function explicitly
with a second derivative equal to an approximate Dirac function at K. We illustrate
this technique with a standard European Call option in the Black-Scholes model.
We computed the Gamma and the sixth derivative with respect to X0 . For the first
derivative, the parameter a does not play an important role but, as we evaluate higher
derivatives, the choice of the parameter a becomes crucial for the quality of a good
approximation and it requires more points to catch the Dirac approximation with
small a. Currently the choice of a is experimental.
We took the same parameters as previously for the standard European Call option
but the maturity for the Gamma now set at T = 5 years and T = 0.2 year for the
sixth derivative with respect to X0 . The initial asset price varies from 1 to 200. The
Monte Carlo parameters are also set to 100, 000 simulation paths and 25 time steps.
The results are displayed on figure 10.
✁ ✆✄

✁

❆✝ ✞✟✠✡
❆☛☞✌✍✎✏✑☞✌ ✒♦✌✓✎✏♦☛

✁ ✆✥

✂

❆✟ ✠✡☛☞ ✌✍✄ ✎✏✑✒
❆✟ ✠✡☛☞ ✓✔✕✌✓ ✌ ✎✏✑✒
❆✖✌✓✗☛✡✕✌✓ ✘✔✓♦☛✡✔✖

✄☎✲ ✄

✁ ✆✂
✁ ✆
✁ ☎
✁ ✄

✲✄☎✲ ✄

✁ ✥
✁

✂

✲ ✁

✂

✲ ✁

✂

✥

✄

☎

✆

✆✂

✆✥

✆✄

✆☎

✂

✲ ✁

✂

✂✄

✥

✆

✝

✞

✂

✂✥

✂✆

✂✝

✂✞

✥

Figure 10. On the left the Gamma versus Price is displayed when computed by AD with the
ramp function (with a = 1); the analytical exact Gamma is also displayed; both curves overlap. On
the right, the sixth derivative with respect to the parameter X0 is displayed when computed via the
same method; the analytical solution is also displayed. We computed the approximation with local
parameter a and with a = 5.

For the Gamma, the curves are overlapping but for the sixth derivative with
respect to the parameter X0 , we cannot take a constant parameter a anymore. When
we choose locally adapted parameter a, the curves are practically overlapping.
4.4. Baskets. A Basket option is a multidimensional derivative security whose
payoff depends on the value of a weighted sum of several risky underlying assets.
As before, Xt is given by (2.3). But now (Wt )t∈[0,T ] is a d-dimensional correlated
Brownian motion with E[dWti dWtj ] = ρi,j dt.
To simplify the presentation, we assume that r and σi are real constants and the
payoff is given by
d
X
ωi Xi T − K)+ ]
VT = e−rT E[(

(4.7)

i=1

Pd
where (ωi )i=1,...,d are positive weights with i=1 ωi = 1. Here, we choose to compare
three different methods. The reference values coming from an approximated momentmatching dynamics (Levy [30] and in Brigo et al. [4]), VAD and second order finite
difference (FD).

19
4.4.1. Algorithm to compute the Gamma of a Basket option. We make
use of the fact that r and σ are constant.
1. Generate M simulation paths using a one time step for the Euler scheme.
X̄ i T ±




d
d
X
X
√
1
Σij T Zj , i = 1, . . . , d,
|Σij |2 T ±
= X i T• exp −
2 j=1
j=1

with XT• = X0 exp (rT ), where Z denotes an N (0; Id ) random vector.
2. For each simulation path, with C = ΣΣT , compute (Vibrato)


T
∂µ
1
√ (VT+ − VT− )C −T Z
∂Xi0
2 h
∂Σ
1
: C −T (ZZ T − Id )C −1
(VT+ − 2VT• + VT− )
+
4h
∂Xi0

∆=

(4.8)

with VT. = (ω · X̄T. − K)+
3. Compute the mean of the resulting vector and discount the result.
4. Apply Automatic Differentiation to what precedes.
4.4.2. Numerical Test. In this numerical test d = 7 and the underlying asset
prices are:
X0 T = (1840, 1160, 3120, 4330.71, 9659.78, 14843.24, 10045.40).

(4.9)

The volatility vector is:
σ T = (0.146, 0.1925, 0.1712, 0.1679, 0.1688, 0.2192, 0.2068).
The correlation matrix is

1.0
0.9477
0.9477
1.0

0.8494 0.7558

0.8548 0.7919

0.8719 0.8209

0.6169 0.6277
0.7886 0.7354

0.8494
0.7558
1.0
0.9820
0.9505
0.6131
0.9303

0.8548
0.7919
0.9820
1.0
0.9378
0.6400
0.8902

0.8719
0.8209
0.9505
0.9378
1.0
0.6417
0.8424

0.6169
0.6277
0.6131
0.6400
0.6417
1.0
0.5927


0.7886
0.7354

0.9303

0.8902
.
0.8424

0.5927
1.0

(4.10)

(4.11)

The number of Monte Carlo paths varies from 1 to 106 with only one time step for
the time integration. Errors are calculated with reference to a solution computed by
approximate moment matching.
On figures 11 and 12, the plot of convergence for the computation of the Gamma
of a Basket made of the first 4 and 7 assets are displayed versus the number of
simulation paths Vibrato plus AD (direct mode) and for Finite differences applied
to a brute force Monte Carlo algorithm. The convergence speed of these methods is
almost the same (with a slight advantage for the Finite difference).
Table 3 displays results for a Basket with the 7 assets, in addition the table 4
displays the CPU time for Vibrato plus AD (direct mode); the finite difference method
is one third more expensive. Again, the method is very accurate.

20
✁

✵ ✵✵✵✆✵

✥✄

❋✠ ✡☛☞✌
❱❆✠ ✡☛☞✌
❆✍✎✏✑✒✓✔✎✏ ✕♦✏✖✒✓♦✍

❋✠ ✡☛☞✌
❱❆✠ ✡☛☞✌

✁

❘✞✍ ❱❛✎✏✞

✥

✵ ✵✵✵✁☎

✁

✂✄
✁

✵ ✵✵✵✁☎

✂

✄☎✲ ✄
✵ ✵✵✵✁✄

✲✄☎✲ ✄
✵ ✵✵✵✁✄

✲ ✁
✲ ✁

✵ ✵✵✵✁✂
✆✵✵✵✵✵

✥✵✵✵✵✵

✝✵✵✵✵✵

✄✵✵✵✵✵

✂
✂✄

✁✞✟✵✝

✥

Figure 11. d=4.

✆

✝

✞

✂☎✟ ✝

Figure 12. d=7.

Convergence of the computation of the Gamma of a Basket option when d = 4 and 7
via Vibrato plus Automatic Differentiation on Monte Carlo and via Finite differences,
versus the number of simulation paths. The parameters are for T = 0.1.
5. American Option. Recall that an American option is like a European option
which can be exercised at any time before maturity. The value Vt of an American
option requires the best exercise strategy. Let ϕ be the payoff, then
Vt := ess sup E[e−r(τ −t) ϕ(Xτ ) | Xt ]

(5.1)

τ ∈Tt

where Tt denotes the set of [t, T ]-valued stopping times (with respect to the (augmented) filtration of the process (Xs )s∈[0,T ] ).
Consider a time grid 0 < t1 < · · · < tn = T with time step h, i.e. tk = kh. To
discretize the problem we begin by assuming that the option can be exercised only at
tk , k = 0, .., n ; its value is defined recursively by

 V̄tn = e−rT ϕ(X̄T )

(5.2)
−rt
 V̄tk = max e k ϕ(X̄tk ), E[V̄tk+1 | X̄tk ] ,
0≤k≤n−1

5.1. Longstaff-Schwartz Algorithm . Following Longstaff et al. [31] let the
continuation value Ctk = E[e−rh V̄tk+1 | X̄tk ] as X is a Markov process. The holder
of the contract exercises only if the payoff at tk is higher than the continuation value
Ctk . The continuation value is approximated by a linear combination of a finite set
of R real basis functions:
Ck ≃

R
X

αk,i ψk,i (X̄tk ).

(5.3)

i=1

Typically, the (αk,i )i=1,...,R are computed by least squares,
 
!2 
R


X
αk,i ψk,i (X̄tk )  .
min E  E[e−rh V̄tk+1 | X̄tk ] −
α 

i=1

(5.4)

21
This leads to a Gram linear system
R
X
j=1


αk,i Gram ψk,i (X̄tk ), ψk,j (X̄tk ) = E[E[e−rh Vk+1 | Xtk ]ψk,i (X̄tk )], i = 1, . . . , R.

(5.5)

Remark 8. Once the optimal stopping time is known, the differentiation with
respect to θ of (5.2) can be done as for a European contract. The dependency of the
τ ∗ on θ is neglected; arguably this dependency is second order but this point needs to
be validated. Hence, the following algorithm is proposed.
5.2. Algorithm to compute the Gamma of an American option.
1. Generate M simulation paths of an Euler scheme with n time steps of size
h = Tn .
2. Compute the terminal value of each simulation path
VT = (K − X̄T )+

(5.6)

3. Compute the Gamma of the terminal condition using (4.3) in section (4.1)
for each simulation path.
4. Iterate from n − 1 to 1 and perform the following at the k-th time step.
(a) Solve the Gram linear system (5.5).
(b) Calculate the continuation value of each path.
Ck+1 (X̄tk ) =

R
X

αk,i ψi (X̄kn ).

(5.7)

i=1

(c) Compute the Gamma by differentiating the Vibrato formula from the
time step k − 1 with respect to X0

N
Zki
1
1 X ∂
n
√
Ȳk−1
(1 + rh) (Ṽki+ − Ṽki− )
Γ̃k =
N i=1 ∂X0
2
X0 σ h

√ 1 i
(Zki )2 − 1
i
i
n
√
.
+ Ȳk−1 σ h (Ṽk+ − 2Ṽk• + Ṽk− )
2
X̄0 σ h
(d) For i = 1, . . . , M
(
Vki = Ṽki , Γik = Γ̃ik

if Ṽki ≥ Ck+1 (X̄kn,i ),

i
Vki = e−rh Vk+1
, Γik = e−rh Γik+1 otherwise

(5.8)
(5.9)

(5.10)

n
with Ṽk+1 = (K − X̄k+1
)+ and

(

√
X̄k± = X̄k−1 + rhX̄k−1 ± σ X̄k−1 hZk
X̄k• = X̄k−1 + rhX̄k−1 .

(5.11)

5. Compute the mean of the vector V and Γ.
Remark 9. The differentiation with respect to X0 is implemented by automatic
differentiation of the computer program.

22
5.2.1. Numerical Test. We consider the following value : σ = 20% or σ = 40%,
X0 varying from 36 to 44, T = 1 or T = 2 year, K = 40 and r = 6%. The Monte
Carlo parameters are: 50, 000 simulation paths and 50 time steps for the time grid.
The basis in the Longstaff-Scharwtz algorithm is (xn )n=0,1,2 .
We compare with the solution of the Black-Scholes partial differential equation
discretized by an implicit Euler scheme in time, finite element in space and semismooth Newton for the inequalities [1]. A second order finite Difference approximation
is used to compute the Gamma. A large number of grid points are used to make it
a reference solution. The parameters of the method are 10, 000 and 50 time steps
per year. Convergence history for Longstaff Schwartz plus Vibrato plus AD is shown
on figure 13 with respect to the number of Monte Carlo paths (Finite Difference on
Monte Carlo is also displayed).
On figure 13, we display the history of convergence for the approximation of
the Gamma of an American Put option versus the number of simulation paths for Vibrato plus Automatic differentiation and for Finite Difference applied to the American
Monte Carlo, the straight line is the reference value computed by PDE+ semi-smooth
Newton. The convergence is faster for VAD than with second order Finite Difference
(the perturbation parameter is taken as 1% of the underlying asset price).
On table 5, the results are shown for different set of parameters taken from
Longstaff et al. [31]. The method provides a good precision when variance reduction (??) is used, for the different parameters, except when the underlying asset price
is low with a small volatility. As for the computation time, the method is faster than
Finite Difference applied to the American Monte Carlo which requires three evaluations of the pricing function whereas VAD is equivalent to two evaluations (in direct
mode).
✥ ✥✝✂
❋✟ ✠✡☛☞❆✌✍✎
❱❆✟ ✠✡☛☞❆✌✍✎

✥ ✥✝

❘✏✑

❱❛✒✓✏

✥ ✥✂✆

✥ ✥✂☎

✥ ✥✂✄

✥ ✥✂✂

✥ ✥✂

✥ ✥✁✆

✥ ✥✁☎

✥ ✥✁✄

✥ ✥✁✂
✥

✁✥✥✥✥

✂✥✥✥✥

✝✥✥✥✥

✄✥✥✥✥

✞✥✥✥✥

Figure 13. Convergence of the Gamma of an American option via Vibrato plus Automatic
Differentiation on the Longstaff-Schwartz algorithm and via Finite Difference, versus the number of
simulation paths. The parameters are σ = 40% and X0 = 40.

6. Second Derivatives of a Stochastic Volatility Model. The Heston model
[23] describes the evolution of an underlying asset (Xt )t∈[0,T ] with a stochastic volatility (Vt )t∈[0,T ] :
p
dWt1 ,
dXt = rXt dt + Vt Xtp
dVt = κ(η − Vt )dt + ξ Vt dWt2 , t ∈ [0, T ]; V0 , X0 given.
(6.1)

Here ξ is the volatility of the volatility, η denotes the long-run mean of Vt and κ the
mean reversion velocity. The standard Brownian process (Wt1 )t∈[0,T ] and (Wt2 )t∈[0,T ]

23
are correlated: E[dWt1 Wt2 ] = ρdt, ρ ∈ (−1, 1). If 2κη > ξ 2 , it can be shown that
Vt > 0 for every t ∈ [0, T ]. We consider the evaluation of a standard European Call
with payoff
VT = E[(XT − K)+ ].

(6.2)

6.1. Algorithm to Compute second derivatives in the Heston Model. To
compute the Gamma by Vibrato method for the first derivative coupled to automatic
differentiation for the second derivative one must do the following:
1. Generate M simulation paths for the underlying asset price (X̄, V̄) and its
X̄,V̄)
tangent process (Ȳ , Ū) = ∂(∂X
using an Euler scheme with n time steps of
0
size h = Tn ,
q

√ 1
n
n
n

V̄kn X̄kn hZ̃k+1
X̄
=
X̄
+
rh
X̄
+
,
X̄0n = X0 ,

k+1
k
k



q
√ 1
n
,
Ȳ0n = 1,
Ȳk+1
= Ȳkn + rhȲkn + V̄kn Ȳkn hZ̃k+1


q

√ 2

 n
V̄k+1 = V̄kn + κ(η − V̄kn )h + ξ V̄kn hZ̃k+1
, V̄0n = V0

with

 1 
1
Z̃
=
ρ
Z̃ 2

p 0
1 − ρ2

  1
Z
Z2

(6.3)

(6.4)

where (Zk1 , Zk2 )1≤k≤n denotes a sequence of N (0; I2 )-distributed random variables.
2. For each simulation path
(a) Compute the payoff
VT = (X̄nn − K)+ .

(6.5)

(b) Compute the Delta using Vibrato at maturity with the n − 1 time steps
and the following formula

with

Z1
1
n
¯ n = Ȳn−1
qn
∆
(1 + rh) (VT+ − VT− )
√
2
n
n
h
X̄n−1
V̄n−1
2
q
√ 1
Zn1 − 1
n
n
q
V̄n−1
h (VT+ − 2VT• + VT− )
+Ȳn−1
√
2
n
n
h
X̄n−1
V̄n−1

 X̄

 X̄

T±
T•

n
n
= X̄n−1
+ rhX̄n−1
±

n
n
= X̄n−1
+ rhX̄n−1
.

q
√ 1
n X̄ n
V̄n−1
n−1 hZ̃n ,

(6.6)

(6.7)

(6.8)

(c) Apply an Automatic Differentiation method on step (2b) to compute
the Gamma.
3. Compute the mean of the result and discount it.

24
6.1.1. Numerical Test. We have taken the following values: the underlying
asset price X0 ∈ [60, 130], the strike is K = 90, the risk-free rate r = 0.135% and the
maturity o is T = 1.
The initial volatility is V0 = 2.8087%, the volatility of volatility is ξ = 1%, the
mean reversion is κ = 2.931465 and the long-run mean is ν = 0.101. The correlation
between the two standard Brownian motions is ρ = 50%.
The number of Monte Carlo path is 500, 000 with 100 time steps each.
The results are displayed on figures 14, 15.
On figure 14 we compare the results obtained by Vibrato plus Automatic Differentiation (direct mode), with second order Finite Difference method applied to a
standard Monte Carlo simulation. On figures 15 we display the Vanna of an Eu✥ ✥✁

✁ ✆✂

❋✡ ☛☞✌✍
❱❆✡ ☛☞✌✍

✥ ✥✆✁

✁ ✆
✁ ☎✂

✥ ✥✆

✥ ✥☎✁

✁ ☎

✥ ✥☎

✁ ✄✂

✥ ✥✄✁

✁ ✄

✥ ✥✄

✁ ✥✂
✁ ✥

✥ ✥✂✁

✁

✥ ✥✂

✂
❋✝ ✞✟✠✡
❱❆✝ ✞✟✠✡
❘☛☞✁ ❱❛✌✍☛

✥ ✥✥✁

✲ ✁

✥
✝✥

✞✥

✟✥

✠✥

✂✥✥

✂✂✥

✂

✂✄✥

✂

✥

✥✂

✄

✄✂

☎

☎✂

✆

✆✂

✂

Figure 14. On the left the Gamma versus Price is displayed when computed by VAD; the
approximated Gamma via Finite Difference is also displayed; both curves overlap. On the right, the
convergence history at one point (X0 , V0 ) = (85, 2.8087) is displayed with respect to the number of
Monte Carlo samples.

☎

☎

❋✠ ✡☛☞✌
❱❆✠ ✡☛☞✌

❋✞ ✟✠✡☛
❱❆✞ ✟✠✡☛
✁ ✁

✁

❘☞✌ ❱❛✍✎☞

✁

✂
✄ ✁

✄
✄

✥
✂ ✁

✲✄
✂

✲✂
✲✁
✲

✥ ✁

✥

✆✥

✝✥

✞✥

✟✥

✄✥✥

✄✄✥

✄✂✥

✄✁✥

✆

✁✆✆✆

✝✆✆✆✆

✝✁✆✆✆

✥✆✆✆✆

✥✁✆✆✆

✂✆✆✆✆

✂✁✆✆✆

✄✆✆✆✆

✄✁✆✆✆

✁✆✆✆✆

Figure 15. On the left the Vanna versus Price is displayed when computed by VAD; the
approximated Vanna via Finite Difference is also displayed; both curves overlap. On the right, the
convergence history at one point (X0 , V0 ) = (85, 2.8087) is displayed with respect to the number of
Monte Carlo samples.

ropean Call option in the Heston model, and again, the convergence with respect to
the number of simulation paths. As for the Gamma, the method is quite precise.
provides a good precision for the approximation of the Vomma and the Vanna. Both
are computed at one point (X0 , V0 ) = (85, 2.8087) with the same set of parameters as

25
given above. The computation by VAD is 30% faster for the Gamma compared with
☎✂✂

✲✄✁✁

❋✟ ✠✡☛☞
❱❆✟ ✠✡☛☞

✆✂✂

✲✄✄✁

✁✂✂

✲✆✁✁

✂
✲✆✄✁

✲✁✂✂
✲☎✁✁

✲✆✂✂
✲☎✄✁

✲☎✂✂

✲✂✁✁

✲✄✂✂

✲ ✁✂✂

❋✠ ✡☛☞✌
❱❆✠ ✡☛☞✌
❘✍✎✏ ❱❛✑✒✍

✲✂✄✁

✲ ✂✂✂

☎✂

✥✂

✄✂

✝✂

✂✂

✂

✁✂

✞✂

✲ ✁✁

✁

✄✁✁✁

✥✁✁✁✁

✥✄✁✁✁

✝✁✁✁✁

✝✄✁✁✁

✞✁✁✁✁

✞✄✁✁✁

✟✁✁✁✁

✟✄✁✁✁

✄✁✁✁✁

Figure 16. On the left the Vomma versus Price is displayed when computed by VAD; the
approximated Vomma via Finite Difference is also displayed; both curves overlap. On the right, the
convergence history at one point (X0 , V0 ) = (85, 2.8087) is displayed with respect to the number of
Monte Carlo samples.

the Vanna. In the case of the Vomma and the Gamma, VAD is 30% faster. For the
Vanna Finite difference requires four times the evaluation of the pricing function so
VAD is twice times faster.
7. Vibrato plus Reverse AD (VRAD). If several greeks are requested at
once then it is better to use AD in reverse mode. To illustrate this point, we proceed
to compute all second and cross derivatives i.e. the following Hessian matrix for a
standard European Call option:

 2
∂ 2V
∂2V
∂2V
∂ V
 ∂X02
∂v∂X0 ∂r∂X0 ∂T ∂X0 

 2
 ∂ V
∂ 2V
∂2V 
∂2V


 ∂X0 ∂σ
∂σ 2
∂v∂r
∂T ∂v 
(7.1)
.
 2
∂ 2V
∂2V
∂2V 
 ∂ V


 ∂X0 ∂r
∂v∂r
∂r2
∂T ∂r 
 ∂ 2V
2
∂ V
∂2V
∂2V 
∂X0 ∂T
∂v∂T
∂r∂T
∂T 2

It is easily seen that a Finite Difference procedure will require 36 (at least 33)
evaluations of the original pricing function whereas we only call this function once if
AD is used in reverse mode. Furthermore, we have to handle 4 different perturbation
parameters.
The parameters are X0 = 90, K = 100, σ = 0.2, r = 0.05 and T = 1 year. The
parameters of Monte Carlo are set to 200, 000 simulation paths and 50 time steps. We
used the library adept 1.0 for the reverse mode. One great aspect here is that we only
have one formula in the computer program to compute all the greeks, consequently
one has just to specify which parameters are taken as variable for differentiation.
The results are shown in the table 1, clearly the reverse automatic differentiation
combined with Vibrato is almost 4 times faster than the finite difference procedures.
8. Malliavin Calculus and Likelihood Ratio Method . Here, we want to
point out that Malliavin calculus and LRM are excellent methods but they have their
own numerical issues especially with short maturities which may make VAD more
attractive for a general purpose software.

26
Mode
Time (sec)

FD (MC)
2.01

VRAD (MC)
0.47

Table 1
CPU time (in seconds) to compute the Hessian matrix of a standard European Call option
(considering X0 , σ, r, T as variables) in the Black-Scholes model.

Let us start by recalling briefly the foundations of Malliavin calculus (further
details are available in Nualart [34], Fournié et al.[10] and in Gobet et al. [18], for
instance). We recall the Bismut-Elworthy-Li formula (see [3], for example):
Proposition 8.1. (Bismut-Elworthy-Li formula) Let X be a diffusion process
given by (2.3) with d = 1, b and σ in C 1 . Let f : R → R be C 1 with E[f (XT )2 ]
and E[f ′ (XT )2 ] bounded. Let (Hth)t∈[0,T ] ani F -progressively measurable process in
RT
L2 ([0, T ] × Ω, dt ⊗ dP) such that E 0 Hs2 ds is finite. Then
#
"
#
"
Z T
Z T
σ(Xs )Hs
′
ds
(8.1)
Hs dWs = E f (XT )YT
E f (XT )
Ys
0
0
dXt
where Yt =
is the tangent process defined in (2.11). By choosing Ht = Yt /σ(Xt )
dx
the above yields


Z


1 T Ys
∂


E [f (XTx )] = E f (XTx )
dWs 
x
∂x
T 0 σ(Xs )


|
{z
}

(8.2)

Malliavin weight

2 #
Yt
is finite.
provided f has polynomial growth and E
σ(Xtx )
0
Second Derivative.. In the context of the Black-Scholes model, the Malliavin
weights, πΓ , for the Gamma is (see [2]):
 2

1
WT
1
πΓ = 2
(8.3)
− − WT .
X0 σT σT
σ
"Z

T



Hence

ΓMal = e

−rT



1
E (XT − K)
X02 σT
+



WT2
1
− − WT
σT
σ



.

The pure likelihood ratio method gives a similar formula (see Lemma 2.1)
!#
"
Z2 − 1
Z
−rT
+
.
ΓLR = e
E (XT − K)
− 2 √
X02 σ 2 T
X0 σ T

(8.4)

(8.5)

LRPW is an improvement of LRM obtained by combining it with a pathwise method
[15].



∂
K
Z
√
ΓLRPW =
= e−rT 2 √ E[Z1{XT >K} ].
e−rT E (XT − K)+
∂X0
X0 σ T
X0 σ T
(8.6)
LRPW is much cheaper than VAD, Malliavin or LRM and it is also less singular at
T = 0. However all these methods require new analytically derivations for each new
problem.

27
8.1. Numerical Tests. We compared VAD with LRPW and Malliavin calculus.
The results are shown on Table 2
T
1.00e+0
5.00e-1
1.00e-1
5.00e-2
1.00e-2
5.00e-3
1.00e-3
5.00e-4
1.00e-4
5.00e-5
1.00e-5

VAD (MC)
3.63e-5
8.55e-5
6.64e-4
1.49e-3
8.78e-3
1.86e-2
9.62e-2
1.85e-1
1.01e+0
1.98e+0
1.03e+1

FD (MC)
1.76e-4
3.11e-4
1.50e-3
2.80e-3
1.84e-2
3.95e-2
1.77e-1
3.34e-1
1.63e+0
3.46e+0
1.78e+1

LRPW (MC)
3.40e-4
7.79e-4
4.00e-3
7.51e-3
3.76e-2
7.55e-2
3.76e-1
7.56e-1
3.77e+0
7.54e+0
3.79e+1

Malliavin (MC)
9.19e-3
1.62e-2
6.54e-2
1.21e-1
5.44e-1
1.10e+0
5.74e+0
1.07e+1
5.26e+1
1.09e+2
5.40e+2

Table 2
Variance of the Gamma of a standard European Call with short maturities in the Black-Scholes
model. Gamma is computed with VAD, FD, LRPW and Malliavin. The computation are done on
the same samples.

The Gamma is computed with the same parameters as in the section 4.3. The
maturity is varying from T = 1 to 10−5 year. The Monte Carlo parameters are also
set to 100, 000 simulation paths and 25 time steps.
Notice the inefficiency of LRPW, Malliavin Calculus and to a lesser degree of
VAD and Finite Difference when T is small.
Note on CPU. Tests have been done on an Intel(R) Core(TM) i5-3210M Processor @ 2,50 GHz. The processor has turbo speed of 3.1 GHz and two cores. We did
not use parallelization in the code.
9. Conclusion. This article extends the work of Mike Giles and investigates the
Vibrato method for higher order derivatives in quantitative finance.
For a general purpose software Vibrato of Vibrato is too complex but we showed
that it is essentially similar to the analytical differentiation of Vibrato. Thus AD of
Vibrato is both general, simple and essentially similar to Vibrato of Vibrato of second
derivatives. We have also shown that Automatic differentiation can be enhanced
to handle the singularities of the payoff functions of finance. While AD for second
derivatives is certainly the easiest solution, it is not the safest and it requires an
appropriate choice for the approximation of the Dirac mass.
Finally we compared with Malliavin calculus and LRPW.
The framework proposed is easy to implement, efficient, faster and more stable
than its competitors and does not require analytical derivations if local volatilities or
payoffs are changed.
Further developments are in progress around nested Monte Carlo and MultilevelMultistep Richardson-Romberg extrapolation [29] (hence an extension to [6]).
Acknowledgment. This work has been done with the support of ANRT and
Global Market Solution inc. with special encouragements from Youssef Allaoui and
Laurent Marcoux.
REFERENCES

28
[1] Y. Achdou and O. Pironneau. Computation methods for option pricing. Frontiers in Applied
Mathematics. SIAM, Philadelphia, 2005. xviii+297 pp., ISBN 0-89871-573-3.
[2] E. Benhamou. Optimal Malliavin weighting function for the computation of the greeks. Mathematical Finance, 13:37–53, 2003.
[3] J. M. Bismut, K. D. Elworthy, and X. M. Li. Bismut type formulae for differential forms.
Probability Theory, 327:87–92, 1998.
[4] D. Brigo, F. Mercurio, F. Rapisarda, and R. Scotti. Approximated moment-matching dynamics
for basket options simulation. Product and Business Development Group, Banca IMI, 2002.
Working paper.
[5] M. Broadie and P. Glasserman. Estimating security price derivatives using simulation. Management Science, 42(2):269–285, 1996.
[6] S. Burgos and M. B. Giles. The computation of greeks with multilevel monte carlo. 2011.
arXiv:1102.1348.
[7] L. Capriotti. Fast greeks by algorithmic differentiation. Journal of Computational Finance,
14(3):3–35, 2011.
[8] L. Capriotti. Likelihood ratio method and algorithmic differentiation: fast second order greeks.
Preprint SSRN:1828503, 2014.
[9] P. S. Dywer and M. S. Macphail. Symbolic matrix derivatives. The Annals of Mathematical
Statistics, 19(4):517–534, 1948.
[10] E. Fournié, J. M. Lasry, J. Lebuchoux, and P. L. Lions. Application of Malliavin calculus to
Monte Carlo methods in finance. Finance and Stochastics, 2(5):201–236, 2001.
[11] M. Giles and P. Glasserman. Smoking adjoints: fast evaluation of greeks in Monte Carlo
calculations. NA-05/15, Numerical Analysis Group, Oxford University, July 2005.
[12] M. B. Giles. Vibrato Monte Carlo sensitivities. In P. L’Ecuyer and A. Owen, editors, Monte
Carlo and Quasi-Monte Carlo Methods 2008, pages 369–382. New York, Springer edition,
2009.
[13] M. B. Giles. Monte carlo evalutation of sensitivities in computational finance. September 20–22,
2007.
[14] P. Glasserman. Gradient estimation via pertubation analysis. Kluwer Academic Publishers,
Norwell, Mass, 1991.
[15] P. Glasserman. Monte Carlo methods in financial engineering, volume 53 of Application of
Mathematics. Springer, New York, 2003. xiii+598 pp., ISBN 0-387-00451-3.
[16] P. Glasserman and X. Zhao. Fast greeks by simulation in forward LIBOR models. Journal of
Computational Finance, 3(1):5–39, 1999.
[17] P. W. Glynn. Likelihood ratio gradient estimation: an overview. In Proceedings of the Winter
Simulation Conference, pages 366–374, New York, 1987. IEEE Press.
[18] E. Gobet and R. Munos. Sensitivity analysing using Itô-Malliavin calculus and martingales:
applications to stochastic optimal control. SIAM Journal on Control and Optimization,
43(5):1676–1713, 2005.
[19] A. Griewank. On automatic differentiation. In M. Iri and K. Tanabe, editors, Mathematical Programming: Recent Developments and Applications, pages 83–108. 1989. Kluwer
Akedemic Publishern Dordrecht.
[20] A. Griewank and A. Walther. Evaluating derivatives: principles and techniques of algorithmic
differentiation. Frontiers in Applied Mathematics. SIAM, Philadelphia, 2008. xxi+426
pp., ISBN 978-0-89871-659-7.
[21] A. Griewank and A. Walther. ADOL-C: A Package for the Automatic differentiation of algorithm written in C/C++. University of Paderborn, Germany, 2010.
[22] L. Hascoët and V. Pascual. The Tapenade automatic differentiation tool: principles, model,
and specification. ACM Transactions On Mathematical Software, 39(3), 2013.
[23] S. L. Heston. A closed-form solution for options with stochastic volatility with applications to
bond and currency options. Review of Financial Studies, 6:327–343, 1993.
[24] Y. C. Ho and X. R. Cao. Optimization and perturbation analysis of queuing networks. Journal
of Optimization Theory and Applications, 40:559–582, 1983.
[25] R. J. Hogan. Fast reverse-mode automatic differentiation using expression templates in C++.
Transactions on Mathematical Software, 40(26):1–26, 2014.
[26] C. Homescu. Adjoints and automatic (algorihmtic) differentiation in computational finance.
arXiv:1107.1831, 2011.
[27] H. Kunita. Stochastic Flows and Stochastic Differential Equations. Cambridge studies in
advanced mathematics. Cambridge University Press, Cambridge, 1990. xiv+361 pp., ISBN
0-521-35050-6.
[28] P. L’Ecuyer. A unified view of the ipa, sf and lr gradient estimation techniques. Management
Science, 36(11):1364–1383, 1990.

29
[29] V. Lemaire and G. Pagès.
Multistep Richardson-Romberg extrapolation.
Preprint
arXiv:1401.1177, 2014.
[30] E. Levy. Pricing European average rate and currency options. Journal of International Money
and Finance, 11:474–491, 1992.
[31] F. A. Longstaff and E. S. Schwartz. Valuing American options by simulation: a simple least
squares approach. Review of Financial Studies, 14:113–148, 2001.
[32] R. L. Mishkov. Generalization of the formula of fa‘a di Bruno for a composite function
with a vector argument. Internation Journal of Mathematics and Mathematical Sciences,
24(7):481–491, 2000.
[33] U. Naumann. The art of differentiating computer programs: an introduction to algorithmic differentiation. Software, Enivronments and Tools. SIAM, RWTH Aachen University, Aachen,
Germany, 2012. xviii+333 pp., ISBN 978-1-61197-206-1.
[34] D. Nualart. The Malliavin calculus and related topics. Probability and its Applications.
Springer-Verlag, Berlin, 2006. x+390 pp., ISBN 978-3-540-28328-7.
[35] O. Pironneau. Automatic differentiation for financial engineering. Université Pierre et Marie
Curie, Paris VI, 2008.
[36] M. Reiman and A. Weiss. Sensitivity analysis for simulations via likelihood ratios. Operations
Research, 37:830–844, 1989.
[37] R. Rubinstein. Sensitivity analysis and performance extrapolation for computer simulation
models. Operations Research, 37:72–81, 1989.
[38] W. Squire and G. Trapp. Using complex variables to estimate derivatives of real functions.
SIAM Review, 40(1):110–112, 1998.
[39] R. Suri and M. Zazanis. Pertubation analysis gives strongly consistent sensitivity estimates for
the m/g/1 queue. Management Science, 34:39–64, 1988.

Table 3
Results for the price, the Delta and the Gamma of a Basket Option priced with the moment-matching approximation (reference values), Finite Difference
on Monte Carlo and Vibrato plus Automatic Differentiation on Monte Carlo. The settings of Monte Carlo simulation are 1 time step and 1, 000, 000 simulation
paths.

d

T

Price
AMM

Price
(MC)

Delta
AMM

Delta
Delta
Vibrato (MC) FD (MC)

Gamma
AMM

Gamma
Gamma
VAD (MC) FD (MC)

1
2
3
4
5
6
7

0.1
0.1
0.1
0.1
0.1
0.1
0.1

38.4285
34.4401
46.0780
59.6741
92.8481
139.235
155.492

37.3823
34.1232
45.9829
58.7849
90.9001
141.766
153.392

0.55226
0.27452
0.18319
0.13750
0.10974
0.09128
0.07820

0.55146
0.27275
0.18220
0.13639
0.10889
0.09017
0.07744

0.55423
0.28467
0.18608
0.14147
0.10956
0.09048
0.07766

4.65557e-3
1.28903e-3
4.29144e-4
1.86107e-4
7.64516e-5
3.54213e-5
2.31624e-5

4.66167e-3
1.34918e-3
4.28572e-4
1.93238e-4
7.79678e-5
3.71834e-5
2.09012e-5

4.64998e-3
1.28193e-3
4.21012e-4
1.79094e-4
7.59901e-5
3.41114e-5
2.18123e-5

1
2
3
4
5
6
7

1
1
1
1
1
1
1

155.389
135.441
181.935
234.985
364.651
543.629
603.818

154.797
133.101
182.642
232.018
363.363
540.870
607.231

0.66111
0.32583
0.21775
0.16304
0.13023
0.10794
0.92420

0.66039
0.32186
0.21497
0.16055
0.12780
0.10477
0.08995

0.67277
0.32547
0.21619
0.01610
0.12804
0.10489
0.89945

1.30807e-3
3.80685e-4
1.26546e-4
5.49161e-5
2.25892e-5
1.04115e-5
6.87063e-6

1.30033e-3
3.86998e-4
1.34423e-4
5.62931e-5
2.38273e-5
8.99834e-6
7.70388e-6

1.32812e-3
3.83823e-4
1.24927e-4
5.50990e-5
2.19203e-5
1.13878e-5
7.22849e-6

30

Table 4
Time computing (in seconds) for the Gamma with Finite Difference on Monte Carlo and with Vibrato plus Automatic Differentiation on Monte Carlo
simulation, dimension of the problem are varying. The settings of Monte Carlo algorithm are the same as above.

Method (Computing Gamma)

d=1

d=2

d=3

d=4

d=5

d=6

d=7

FD (MC)
VAD (MC)

0.49
0.54

0.95
0.77

1.33
0.92

1.82
1.21

2.26
1.50

2.91
1.86

3.36
2.31

Table 5
Results of the price, the Delta and the Gamma of an American option. The reference values are obtained via the Semi-Newton method plus Finite Difference,
they are compared to Vibrato plus Automatic Differentiation on the Longstaff-Schwartz algorithm. We compute the standard error for each American Monte
Carlo results. The settings of the American Monte Carlo are 50 time steps and 50, 000 simulation paths.

Price
Price Standard
Ref. Value (AMC)
Error

Delta
Delta
Standard
Ref. Value Vibrato (AMC)
Error

Gamma
Gamma
Standard
Ref. Value VAD (AMC)
Error

S

σ

T

36
36
36
36

0.2
0.2
0.4
0.4

1
2
1
2

4.47919
4.83852
7.07132
8.44139

4.46289
4.81523
7.07985
8.45612

0.013
0.016
0.016
0.024

0.68559
0.61860
0.51019
0.44528

0.68123
0.59934
0.51187
0.44102

1.820e-3
1.813e-3
1.674e-3
1.488e-3

0.08732
0.07381
0.03305
0.02510

0.06745
0.06398
0.03546
0.02591

6.947e-5
6.846e-5
4.852e-5
5.023e-5

38
38
38
38

0.2
0.2
0.4
0.4

1
2
1
2

3.24164
3.74004
6.11553
7.59964

3.23324
3.72705
6.11209
7.61031

0.013
0.015
0.016
0.025

0.53781
0.48612
0.44726
0.39786

0.53063
0.46732
0.45079
0.39503

1.821e-3
1.669e-3
1.453e-3
1.922e-3

0.07349
0.05907
0.02989
0.02233

0.07219
0.05789
0.03081
0.02342

1.198e-4
1.111e-4
5.465e-5
4.827e-5

40
40
40
40

0.2
0.2
0.4
0.4

1
2
1
2

2.31021
2.87877
5.27933
6.84733

2.30565
2.86072
5.28741
6.85873

0.012
0.014
0.015
0.026

0.41106
0.38017
0.39051
0.35568

0.40780
0.39266
0.39485
0.35446

1.880e-3
1.747e-3
1.629e-3
1.416e-3

0.06014
0.04717
0.02689
0.01987

0.05954
0.04567
0.02798
0.02050

1.213e-4
5.175e-4
1.249e-5
3.989e-5

42
42
42
42

0.2
0.2
0.4
0.4

1
2
1
2

1.61364
2.20694
4.55055
6.17459

1.60788
2.19079
4.57191
6.18424

0.011
0.014
0.015
0.023

0.30614
0.29575
0.33973
0.31815

0.29712
0.28175
0.34385
0.29943

1.734e-3
1.601e-3
1.517e-3
1.347e-3

0.04764
0.03749
0.02391
0.01768

0.04563
0.03601
0.02426
0.01748

4.797e-5
5.560e-5
3.194e-5
2.961e-5

44
44
44
44

0.2
0.2
0.4
0.4

1
2
1
2

1.10813
1.68566
3.91751
5.57268

1.09648
1.66903
3.90838
5.58252

0.009
0.012
0.015
0.028

0.21302
0.22883
0.29466
0.28474

0.20571
0.21972
0.29764
0.28447

1.503e-3
1.487e-3
1.403e-3
1.325e-3

0.03653
0.02960
0.02116
0.01574

0.03438
0.02765
0.02086
0.01520

1.486e-4
2.363e-4
1.274e-4
2.162e-4

31

☎✂✁
✥✁✂✞✂✂✥
✥✄✁
✆✄
✄✂
✝✁✁

❱❆✟ ✠ ❆✝✞✟✠✡☛☞✞✟
✡❱❆✠
✠❡ ✡☛☞✝
☛☞✌
❱❆✆
❱❆✞
✌♦✟✍✡☛♦✝
✍▼✎✏
✌✡✍✮
✝✞✟✠
✟✠✡☛
❱❆✟
❱❆✞
❆✎✎✏♦①☛✑✞✡✒✓
❆✡☛☞✌✍✎✏☛☞
❆☞✌✍✎✏✑✒✌✍
✌✑✠
✌✑✟✠✒✓
❱❆✠
☛☞✌✓♦✍✔✏✑♦☞
☛☞✌✠
✌♦✟✍✡☛♦✝
✍▼✎✏
✌✡✍✮
✍▼✎✏
✒✓ ✡☛☞✞
✑♦☞✒✍✎♦✡
❆✎✏✑✒✓✔✕✏✑ ❡☞✖♦✓✑☞✔
❆✔✕✖✗✓✑✒✕✖
✖♦✑✗✓✔♦✎
❱❆✞ ✟✠✡✠✟❡✠ ☛☞✌✠ ✍▼✎✏

✆
☎✥✁
✁ ✆✥
✆✁✁ ✄
✂✞✂✂✂✥
✥ ✥✄
☎ ✁✁✆✂
☎
☎✁✁
✁ ✁✥✂
✆
☎✁
✥✥✥✂✁
✁✂✝
✁✄ ☎
✄✁
✁ ✥
✄✁
✂✁
✁ ✄
✥✥ ✥✂
✁✂✆
✂
✁ ✥✂
✥✁
✂✁
✁✁ ✂
✥✥✥✥✁
✁✂☎
✁
✁
✲✲✥✥✁✁ ✁✂✄
✥
✂
✁ ✥✥✂

✁✄✄✥
✂

✝✁✁✁✁
✥✁✥☎☎✥

✥✂✂✂✁✄✆✆✥

✞✁✁✁✁
✄✁☎✝✝✥

☎✁✁
✆✥✂✂✂ ✟✁✁✁✁
✥✂✥✥
☎✥✄
✆✂✁
✂✄✥

☎✥✁
✆✥ ✥✂✂✂✂
✥☎
✥✁✁✁✁
☎✂✁
✆✄
✥✆
✂☎✥
✂✆✥

☎✄✁
✆☎
✥✝
✂✝✥

✥✂✂✂✂✂
✄✄✥✥
✂✁✁
✁✁✁✁

