Weak and Strong Taylor methods for numerical solutions of stochastic
  differential equations
October 31, 2018
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL
SOLUTIONS OF STOCHASTIC DIFFERENTIAL EQUATIONS
MARIA SIOPACHA AND JOSEF TEICHMANN
Abstract. We apply results of Malliavin-Thalmaier-Watanabe for strong and
weak Taylor expansions of solutions of perturbed stochastic differential equa-
tions (SDEs). In particular, we work out weight expressions for the Taylor
coefficients of the expansion. The results are applied to LIBOR market models
in order to deal with the typical stochastic drift and with stochastic volatility.
In contrast to other accurate methods like numerical schemes for the full SDE,
we obtain easily tractable expressions for accurate pricing. In particular, we
present an easily tractable alternative to “freezing the drift” in LIBOR market
models, which has an accuracy similar to the full numerical scheme. Numerical
examples underline the results.
1. Introduction and Setting
Let (Ω,F ,P) be a probability space carrying an N -dimensional Brownian motion
(Wt)t≥0 with a d × d correlation matrix. We consider smooth curves Fǫ : R →
L2(Ω;RN ) of random variables, where ǫ ∈ R is a parameter. We apply Taylor
theorems to obtain strong approximations of the curve Fǫ at ǫ = 0 and we apply
partial integration on Wiener space to obtain weak approximations of the law of Fǫ
for small values of ǫ.
We choose the notion Taylor expansion instead of asymptotic expansion in order
to point out that the strong method is indeed a classical Taylor expansion with usual
conditions for convergence. The weak method represents a truncated converging
power series in the parameter ǫ if – for instance – the payoff f : RN → R stems
from a real analytic function and some distributional properties are satisfied.
2. Weak and strong Taylor methods - Structure Theorems
We introduce in this section two concepts of approximation. Consider a curve
ǫ 7→ Fǫ, where ǫ ∈ R and Fǫ ∈ L2(Ω;RN ).
Definition 1. A strong Taylor approximation of order n ≥ 0 is a (truncated)
power series
(2.1) Tnǫ (Fǫ) :=
such that
(2.2) E
|Fǫ −Tnǫ (Fǫ)|
= o(ǫn),
Financial support from the Austrian Science Fund (FWF) under grant P 15889 and the START-
prize-grant Y328-N13 is gratefully acknowledged. Furthermore this work was financially supported
by the Christian Doppler Research Association (CDG). The authors gratefully acknowledge a
fruitful collaboration and continued support by Bank Austria and the Austrian Federal Financing
Agency (ÖBFA) through CDG.
http://arxiv.org/abs/0704.0745v1
2 MARIA SIOPACHA AND JOSEF TEICHMANN
holds true as ǫ→ 0.
Remark 1. In our setting a strong Taylor approximation of any order n ≥ 0 of
the curve Fǫ can always be obtained, see for instance [KM97].
Let f : RN → R be a Lipschitz function with Lipschitz constant K, then we
obtain
(2.3) E
|f(Fǫ)− f
ǫ (Fǫ)
‖Fǫ −Tnǫ (Fǫ)‖
= Ko(ǫn).
Equation (2.3) does not hold anymore if f is not globally Lipschitz continuous.
In particular, we observe the dependence of the right hand side on the Lipschitz
constant K. Hence, truncating an a-priori known Taylor expansion leads to an
error term, which contains the Lipschitz constant and is therefore not useful for
non-Lipschitz claims. The weak method navigates around this feature by partial
integration.
Definition 2. A weak Taylor approximation of order n ≥ 0 is a power series
for each bounded, measurable f : RN → R,
ǫ (f, Fǫ) :=
E(f(F0)πi),
where πi ∈ L1(Ω) denote real valued, integrable random variables, such that
f(Fǫ)
−Wnǫ (f, Fǫ)| = o(ǫn).
Remark 2. The weights πi for i ≥ 1 are called Malliavin weights.
Remark 3. If the law of Fǫ is real analytic at ǫ = 0 in the weak sense, i.e. if there
exist (signed) measures µi such that for all bounded, measurable f : R
N → R the
following series converges and the equality
f(Fǫ)
f(x)µi(dx),
holds true, precisely then we do have a converging weak Taylor expansion. We aim
for constructing stochastic representations of the following type, for i ≥ 0:
f(x)µi(dx) = E(f(F0)πi).
For the definition of the weak Taylor approximation to make sense, existence of
the Malliavin weights has to hold. The following theorem can be found in a slightly
different version in [MT06] and goes back to S. Watanabe. For the definition and
notion of D∞(RN ) see [Mal97] or [Nua06].
Theorem 1. Let Fǫ : R → D∞(RN ) be smooth and assume that the Malliavin co-
variance matrix γ(Fǫ) is invertible with p-integrable inverse for every p ≥ 1 around
ǫ = 0 (i.e. on an open interval containing ǫ = 0). Then there is a weak Taylor
approximation of any order n ≥ 0 and there are explicit formulas for the weights
πi. If we only know that the Malliavin covariance matrix γ(F0) is invertible with
p-integrable inverse, then we can also calculate the Malliavin weights, since they
depend only on γ(F0).
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 3
Proof. Fix n ≥ 0 and take a smooth test function f : RN → R and assume that
γ−1(Fǫ) exists as a smooth curve in D∞ on a open ǫ-interval containing ǫ = 0. By
standard arguments we can prove the following formula
f(Fǫ)
f(Fǫ)δ
s 7→ (DsFǫ)Tγ−1(Fǫ)
More precisely, by the integration by parts [Nua06, Definition 1.3.1-(1.42)], the
chain rule [Nua06, Proposition 1.2.3] and the definition of the Malliavin covariance
matrix, [Nua06, page 92], we obtain from the right hand side the desired left hand
side.
Notice that the ǫ-dependence of the Skorohod integral is smooth due to basic
properties of D∞. Hence, we can calculate higher derivatives of the left hand side
by iterating the above procedure and differentiating the Skorohod integral. We
denote
(2.4) π1 := δ
s 7→ (DsFǫ)Tγ−1(Fǫ)
We write then, pars pro toto, the formula for the second derivative
E(f(Fǫ)) = E
f(Fǫ)δ
s 7→ π1(DsFǫ)Tγ−1(Fǫ)
f(Fǫ)δ
s 7→ (Ds
γ−1(Fǫ)
f(Fǫ)δ
s 7→ (DsFǫ)Tγ−1(Fǫ)
dγ(Fǫ)
γ−1(Fǫ)
f(Fǫ)δ
s 7→ (DsFǫ)Tγ−1(Fǫ)
This formula makes perfect sense at ǫ = 0 and – by induction – we see that we can
perform this step for any derivative. The general, recursive result is the following:
as := (DsFǫ)
γ−1(Fǫ)
for 0 ≤ s ≤ T,
πn := δ(s 7→ asπn−1) +
πn−1,
π0 := 1.
Here we understand the weights πn as ǫ-dependent, whereas in the final formulas
we put ǫ = 0. This proves the result for smooth test functions f and under the
assumption that the Malliavin covariance matrix is invertible around ǫ = 0. If we
approximate a bounded, measurable function f by smooth test functions we obtain
the desired assertion by standard arguments, since the weights are integrable. �
Remark 4. By Taylor’s theorem and the Faà-di-Bruno-formula we obtain
dnf(Fǫ)
|α|≤n
f (α)(Fǫ)pα,
where pα is a well-defined polynomial in derivatives of the curve ǫ 7→ Fǫ, for a
multi-index α. Since D∞ is an algebra, see [Mal97], the above expression lies in
4 MARIA SIOPACHA AND JOSEF TEICHMANN
Lp(Ω) for each p ≥ 0. The previous result provides a representation of the partial
integration result for
|α|≤n
f (α)(Fǫ)pα) = E(f(Fǫ)πn).
The structure of the weights is seen from above. The result can be considered as
a dual version of the Faà-di-Bruno-formula. However, the structure of this dual
formula is much simpler.
We provide an example to demonstrate the strong and weak method of approx-
imation. The method works in order to replace time-consuming iteration schemes,
like the Euler-scheme, by simulations of “simple” Itô integrals.
Example 1. We deal with a generic, real-valued random variable over a one-
dimensional Gaussian space, see [Nua06], i.e.
where the F i lie in the (i+1)st Wiener chaos Hi+1(Ω) (one can think of a Hermite
expansion for instance) and the sum is understood in the L2-sense. From the strong
expansion we obtain immediately – for a given Lipschitz function f : R → R – that
f(Fǫ)
f(F 0 + ǫF 1)
| ≤ Ko(ǫ),
as ǫ→ 0, where K denotes the Lipschitz constant of f . This simple approximation
can be sometimes quite useful.
We assume now that F 0 =
h(s)dWs has non-vanishing variance in order to
calculate the weights, which do depend only on γ(F0). The strong Taylor approxi-
mation is given by definition, the weak Taylor expansion can be constructed by the
previous recursive formulas and the specifications
0 = h(s),
γ(F 0) =
h(s)2ds,
h(s)2ds
In order to obtain a first-order approximation for bounded, measurable random
variables we therefore have to calculate
f(F 0)
f(F 0)π1
where
π1 = δ
s 7→ asF 1
This amounts to an integration of f times a polynomial with respect to a Gaussian
density, since:
f(F 0)π1
f(F 0)F 1
asdWs
f(F 0)DsF
Notice that the strong approximation does not yield such a result for bounded,
measurable random variables. Notice also that in the given case the approximation
can be calculated in a deterministic way, since we deal with Gaussian integrations.
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 5
The second-order weak Taylor approximation is given by
f(F 0)
f(F 0)π1
+ ǫ2E
f(F 0)π2
where
π2 = δ
s 7→ π1asF 1
s 7→ asF 2
3. Applications from Financial Mathematics
For applications we want to deal with strong and weak Taylor approximations
of a given curve of random variables. We are particulary interested in cases, where
the first derivative dFǫ
|ǫ=0 is of simple form or – even more important – where the
Malliavin covariance matrix γ(F0) is of simple form. In these cases it is easy to
obtain first or second order approximations of the respective quantities in the weak
or strong sense.
In what follows, first we will present one of the most applied interest rate models,
namely the LIBOR market model (LMM). Then, we will introduce the commonly
used technique of freezing the drift. We will show how to embed the”freezing the
drift” technique into our framework of Taylor approximations. We understand
freezing the drift as a strong Taylor approximation of order zero in the drift term
of the LIBOR SDE. Our goal is to put this technique into a method, where we
can in particular improve the order of approximation. We will finally extend the
assumption of log normality and develop a stochastic volatility LMM, where we will
show how to obtain tractable option prices via our weak Taylor approximations.
3.1. The LIBOR Market Model. We apply our concepts to the LMM, initially
constructed by [BGM97], [MSS97] and [Jam97]. Let T denote a strictly positive
fixed time horizon and (Ω,FT ,P, (Ft)0≤t≤T ) be a complete probability space, sup-
porting an N -dimensional Brownian motion Wt = (W
t , ...W
t )0≤t≤T . The factors
are correlated with dW it dW
t = ρijdt. Let 0 = T0 < T1 < T2 < . . . < TN <
TN+1 =: T be a discrete tenor structure and α := Ti+1 − Ti the accrual factor for
the time period [Ti, Ti+1], i = 0, . . . , N . Let P (t, Ti) denote the value at time t
of a zero coupon bond with maturity Ti ∈ [0, T ]. The measure P is the terminal
forward measure, which corresponds to taking the final bond P (t, T ) as numéraire.
The forward LIBOR rate Lit := Lt(Ti, Ti+1) at time t ≤ Ti for the period [Ti, Ti+1]
is given by:
Lit = Lt(Ti, Ti+1) =
P (t, Ti)
P (t, Ti+1)
We assume that for any maturity Ti there exists a bounded, continuous, determin-
istic function σi(t) : [0, Ti] → R, which represents the volatility of the LIBOR Lit,
i = 1, ..., N . The log normal LIBOR market model can be expressed under the
measure P as:
(3.1) dLit = σ
i(t)Lit
j=i+1
1 + αL
dt+ σi(t)LitdW
t , i = 1, ..., N.
3.2. Freezing the Drift. The dynamics of forward LIBORs for i = 1, ..., N − 1
depend on the stochastic drift term
, i ≤ j ≤ N , which is determined by LI-
BOR rates with longer maturities. This random drift prohibits analytic tractability
when pricing products that depend on more that one LIBOR rate, since there is
no unifying measure under which all LIBOR rates are simultaneously log normal.
6 MARIA SIOPACHA AND JOSEF TEICHMANN
In addition, it encumbers the numerical implementation of the model. Common
practice is to approximate this term by its starting value
or as it is widely
referred to as freezing the drift, i.e.
1 + αL
1 + αL
It was first implemented in the original paper [BGM97] for the pricing of swaptions
based on the LMM. [BW00] and [Sch02] argue that freezing the drift is justified due
to the fact that this term has small variance. However, by freezing the drift there
is a difference in option prices with the real and the frozen drift. It has not been
examined how big the error is or for which assets it works well or not. Our aim is
to investigate such a phenomenon and improve the performance by providing with
correction terms of order one.
3.3. Correcting the Frozen Drift. The purpose of this section is to embed the
well-known and often applied technique of freezing the drift into the strong and
weak Taylor approximations, in order to develop a method to improve the order of
accuracy. Specifically for the strong Taylor approximation, the method works well,
since we always deal with a globally Lipschitz drift term x 7→ αx+
1+αx+
with small
Lipschitz constant α.
Remark 5. As it will be clear later, the strong Taylor correction method can be
accommodated with any extension of the log normal LMM, for example with the
Lévy LIBOR model by Eberlein and Özkan [EÖ05].
3.3.1. Strong Taylor Approximation. We first state a useful lemma, asserting that
we can indeed freeze the drift under special model formulation and choice parame-
ters.
Lemma 1. Let ǫ1 ∈ R and consider for i = 1, . . . , N the following stochastic
differential equation:
(i,ǫ1)
t = ǫ1
σi(t)X
(i,ǫ1)
j=i+1
(j,ǫ1)
1 + αX
(j,ǫ1)
ρijdt+ dW
,(3.2)
defined on the complete probability space (Ω,FT ,P, (Ft)0≤t≤T ) where Wt is an N -
dimensional Brownian motion under the measure P with dW it dW
t = ρijdt. Then
the first-order strong Taylor approximation for X
(i,ǫ1)
t is given by:
(3.3) T1ǫ1(X
(i,ǫ1)
t ) = X
(i,0)
t + ǫ1
(i,ǫ1)
Proof. By (1) we obtain for n = 1:
(i,ǫ1)
t ) ≃ X
(i,ǫ1)
t = X
(i,0)
0 + ǫ1Y
t + o(ǫ1),
since X
(i,0)
t = X
(i,0)
0 and where Y
|ǫ1=0X
(i,ǫ1)
t is the first-order correction
term. By differentiating (3.2) with respect to ǫ1, we calculate:
(i,ǫ1)
= σi(t)X
(i,0)
j=i+1
(j,0)
1 + αX
(j,0)
ρijdt+ dW
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 7
and derive Y it as the solution to the above linear SDE:
(3.4) Y it =
−σi(s)X(i,0)0
j=i+1
(j,0)
1 + αX
(j,0)
σi(s)X
(i,0)
with Y i0 = 0. �
Remark 6. We parametrise the LIBOR market model in terms of the parameter
ǫ1 as follows:
(i,ǫ1)
t = σ
i(t)L
(i,ǫ1)
j=i+1
(j,ǫ1)
1 + αX
(j,ǫ1)
ρijdt+ dW
and assume at t = 0 that L
(i,ǫ1)
0 = X
(i,ǫ1)
0 for all ǫ1 and all i = 1, ..., N . If ǫ1 = 1,
what we obtain is the standard LIBOR market model formulation and in particular
(i,1)
t = X
(i,1)
t . For ǫ1 = 0, X
(i,0)
t equals its starting value and thus the drift term
in the following SDE is no longer stochastic:
(3.5) dL
(i,0)
t = σ
i(t)L
(i,0)
j=i+1
(j,0)
1 + αX
(j,0)
ρijdt+ dW
The next proposition provides a way for a pathwise approximation of L
(i,ǫ1)
t , by
means of adjusting its SDE. This is achieved by adding Tnǫ1(X
(j,ǫ1)
t ) in the frozen
drift part.
Proposition 1. Assume the setup of Lemma 1 and assume further at t = 0 that
(i,ǫ1)
0 = X
(i,ǫ1)
0 for all ǫ1 and all i = 1, ..., N . Then the stochastic differential
equation for L
(i,ǫ1)
t with the unfrozen drift:
(3.6) dL
(i,ǫ1)
t = σ
i(t)L
(i,ǫ1)
j=i+1
(j,ǫ1)
1 + αX
(j,ǫ1)
ρijdt+ dW
can be strongly approximated as ǫ1 ↓ 0 by
(3.7) dL̂
(i,ǫ1)
t = σ
i(t)L̂
(i,ǫ1)
j=i+1
Tnǫ1(X
(j,ǫ1)
σj(t)
1 + α
Tnǫ1(X
(j,ǫ1)
ρijdt+ dW
Remark 7. For n = 0, we derive the ”freezing the drift” case. For n = 1, we
already obtain an improvement.
Proof. First step is to interchange X
(j,ǫ1)
t with (X
(j,ǫ1)
t )+ in (3.6) to obtain:
(i,ǫ1)
t = σ
i(t)L
(i,ǫ1)
j=i+1
(j,ǫ1)
t )+σ
1 + α(X
(j,ǫ1)
ρijdt+ dW
This yields no change for the dynamics of L
(i,ǫ1)
t , since X
(j,ǫ1)
t = (X
(j,ǫ1)
t )+.
8 MARIA SIOPACHA AND JOSEF TEICHMANN
By Taylor’s expansion, we know that as ǫ1 ↓ 0, L̂(i,ǫ1)t → L
(i,ǫ1)
t P-a.s. The
estimate for the error term is given by
log L̂
(i,ǫ1)
t − logL
(i,ǫ1)
σi(s)
j=i+1
Tnǫ1(X
(j,ǫ1)
σj(s)
1 + α
Tnǫ1(X
(j,ǫ1)
ρij +
j=i+1
(j,ǫ1)
t )+σ
1 + α(X
(j,ǫ1)
α|X(j,ǫ1)s − (Tnǫ1(X
(j,ǫ1)
s ))+|ds.
Remark 8. The SDE for the approximated L̂
(i,ǫ1)
t is easier and faster to simulate
than (3.1), as it is exhibited by the following example. Notice additionally that
(i,ǫ1)
t is a continuous functional of the process Y
t (3.4) and of the Brownian path
W it . Eventually, by using L̂
(i,ǫ1)
t as the LIBOR rates, the computational complexity
of the drift and thus of the model can be reduced substantially, while maintaining
accuracy of prices.
Example 2. In this example, we examine the performance of the strong Taylor
correction method. Let N = 3 and consider pricing a caplet on the LIBOR rate L1
with strike K. Its price is given by:
0 = αEP
L1T1 −K
Assume that the volatility functions σi(t) : [0, Ti] → R for i = 1, 2, 3 are given by
(cf. Brigo and Mercurio [BM01], formulation (6.12)):
σi(t) =
a(Ti − t) + d
− b(Ti − t)
where the constants a, b, d, e are the same for all three LIBOR rates and are equal
to a = −0.113035, b = 0.22911, d = −a, e = 0.684784. Thus, we can write the
model under the terminal measure P as:
(1,ǫ1)
t = σ
1(t)L
(1,ǫ1)
(2,ǫ1)
2(t)ρ12
1 + αX
(2,ǫ1)
(3,ǫ1)
3(t)ρ13
1 + αX
(3,ǫ1)
dt+ σ1(t)L
(1,ǫ1)
(2,ǫ1)
t = σ
2(t)L
(2,ǫ1)
(3,ǫ1)
3(t)ρ23
1 + αX
(3,ǫ1)
dt+ σ2(t)L
(2,ǫ1)
dL3t = σ
3(t)L3t dW
(1,ǫ1)
t = ǫ1
σ1(t)X
(1,ǫ1)
(2,ǫ1)
2(t)ρ12
1 + αX
(2,ǫ1)
(3,ǫ1)
3(t)ρ13
1 + αX
(3,ǫ1)
dt+ σ1(t)X
(1,ǫ1)
(2,ǫ1)
t = ǫ1
σ2(t)X
(2,ǫ1)
(3,ǫ1)
3(t)ρ23
1 + αX
(3,ǫ1)
dt+ σ2(t)X
(2,ǫ1)
(3,ǫ1)
t = ǫ1
σ3(t)X
(3,ǫ1)
with initial values L
(i,ǫ1)
0 = X
(i,ǫ1)
0 = ci, for i = 1, 2, 3 and for all ǫ1. The Brownian
motion vector (W 1t ,W
t ) is correlated with correlation coefficient ρij given by:
ρij = 0.49 + (1− 0.49) exp (−0.13|i− j|), i, j = 1, 2, 3.
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 9
The SDEs for the approximated LIBOR rates L̂
(1,ǫ1)
t and L̂
(2,ǫ1)
t are given by:
(1,ǫ1)
t = σ
1(t)L̂
(1,ǫ1)
c2 + ǫ1Y
σ2(t)ρ12
1 + α
c2 + ǫ1Y
c3 + ǫ1Y
σ3(t)ρ13
1 + α
c3 + ǫ1Y
+ σ1(t)L̂
(1,ǫ1)
(2,ǫ1)
t = σ
2(t)L̂
(2,ǫ1)
c3 + ǫ1Y
σ3(t)ρ23
1 + α
c3 + ǫ1Y
dt+ σ2(t)L̂
(2,ǫ1)
The partial derivative terms Y 2t and Y
t are equal to:
Y 2t = c2
( ∫ t
σ2(s)dW 2s −
αc3ρ23
1 + αc3
σ2(s)σ3(s)ds
Y 3t = c3
σ3(s)dW 3s .
We compare three caplet prices:
• benchmark price, underlying L(1,ǫ1)t ;
• strong Taylor price, underlying L̂(1,ǫ1)t ;
• frozen drift price, underlying L(1,0)t .
Numerical results in basis points (bps) are displayed in Table 2 for parameters
ǫ1 = 1, N = 3, α = 0.50137, c1 = 3.86777%, c2 = 3.7574%, c3 = 3.8631%,
T1 = 1.53151, Ti = T1 + iα, i = 2, 3, 4. We characteristically observe the difference
in prices between the benchmark and frozen drift price, whilst our strong Taylor
correction method performs very well and is computationally simpler and faster.
strikes K=3% K=3.5% K=4% K=5.75% K=6.25% K=8%
benchmark 11.1831 8.5897 6.5503 3.0349 2.4423 1.2969
strong Taylor 11.0687 8.5691 6.5867 3.1448 2.5513 1.3926
frozen drift 13.9551 11.1822 8.8803 4.6313 3.8506 2.2524
Table 1: Caplet values in bps for parameters ǫ1 = 1, α = 0.50137, c1 = 3.86777%, c2 = 3.7574%,
c3 = 3.8631% and T1 = 1.53151.
3.3.2. Weak Taylor Approximation. In what follows, we provide some results on
how to correct option prices obtained by the SDE with the frozen drift (3.5) by
adding a correction term involving the appropriate Malliavin weight. Let L
i,k,ǫ1
denote the vector of the LIBOR rates (L
(i,ǫ1)
, . . . , L
(k,ǫ1)
Proposition 2. Assume the setup of Lemma 1, where the ith LIBOR rate is given
(3.8) dL
(i,ǫ1)
t = σ
i(t)L
(i,ǫ1)
j=i+1
(j,ǫ1)
1 + αX
(j,ǫ1)
ρijdt+ dW
with L
(i,ǫ1)
0 = X
(i,ǫ1)
0 for all ǫ1 and all i = 1, ..., N . Assume furthermore that the
Malliavin covariance matrix γ(L
i,k,0
) is invertible. Then the price of an option with
10 MARIA SIOPACHA AND JOSEF TEICHMANN
payoff g(L
i,k,ǫ1
), for i ≤ k ≤ N and g bounded measurable, can be approximated by
the weak Taylor approximation of order one:
a(g,L
i,k,ǫ1
) = P (0, T )
i,k,0
+ ǫ1EP
i,k,0
,(3.9)
where the Malliavin weight ζTi is given by:
ζTi = δ
i,k,0
)Tγ−1(L
i,k,0
i,k,ǫ1
,(3.10)
for t ≤ Ti.
Proof. The weight ζTi is obtained by (2.4). Notice that we can write:
i,k,ǫ1
i,k,0
and hence the result (3.9) by Definition 2 for n = 1. �
Example 3. In this example we let N = 3 and we price a payers swaption with
strike price K and maturity T1, where the underlying swap is entered at T1 and has
payment dates T2 and T3. We assume that the volatility functions σ
i(t) : [0, Ti] → R
for i = 1, 2, 3 are constant:
σ1(t) = σ1, σ
2(t) = σ2, σ
3(t) = σ3,
such that we obtain under the terminal measure P:
(1,ǫ1)
t = σ1L
(1,ǫ1)
(2,ǫ1)
1 + αX
(2,ǫ1)
(3,ǫ1)
1 + αX
(3,ǫ1)
dt+ dW 1t
(2,ǫ1)
t = σ2L
(2,ǫ1)
(3,ǫ1)
1 + αX
(3,ǫ1)
dt+ dW 2t
dL3t = σ3L
t ,(3.11)
(1,ǫ1)
t = ǫ1
(1,ǫ1)
(2,ǫ1)
1 + αX
(2,ǫ1)
(3,ǫ1)
1 + αX
(3,ǫ1)
dt+ dW 1t
(2,ǫ1)
t = ǫ1
(2,ǫ1)
(3,ǫ1)
1 + αX
(3,ǫ1)
dt+ dW 2t
(3,ǫ1)
t = ǫ1
(3,ǫ1)
with initial values L
(i,ǫ1)
0 = X
(i,ǫ1)
0 = ci, for i = 1, 2, 3 and for all ǫ1. W
t and
W 2t are correlated with correlation coefficient ρ12. We freeze the drifts in the above
equations to obtain:
(1,0)
t = c1 exp
( αc2σ2
1 + αc2
αc3σ3
1 + αc3
(2,0)
t = c2 exp
( αc3σ3
1 + αc3
L3t = c3 exp
Similarly to the previous example, we compare four option prices:
• benchmark price;
• frozen drift;
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 11
• strong Taylor price;
• weak Taylor price.
The weak correction formula (3.9) adds a correction term to the closed form price
of the option. The swaption payoff at Ti can be found for example in [MR98]:
swptn
k=i+1
(1 + αL
if the underlying swap is entered at time Ti and has payment dates Ti+1, ..., T . αk
is given by:
Kα, k = i+ 1, . . . , N,
1 +Kα, k = N + 1.
The payers swaption value at time t = 0 can be written as:
swptn
0 = P (0, Ti)EPi
swptn
= P (0, T )EP
(1 + αL
)− (1 +Kα)
,(3.12)
where αi := −1 and Pi denotes the forward measure corresponding to the bond
P (t, Ti) as numéraire. Therefore, its benchmark price is given by the above formula
with N = 2 and i = 1:
swptn
0 = P (0, T )
αL1T1 + αL
+ α2L1T1L
−Kα2L2T1 − 2Kα
Its weak Taylor price is given by (3.9) with i = 1 and k = N = 2.
swptn
0 = P (0, T )
(1,0)
(2,0)
+ α2L
(1,0)
(2,0)
−Kα2L(2,0)
− 2Kα
+ ǫ1EP
(1,0)
(2,0)
+ α2L
(1,0)
(2,0)
−Kα2L(2,0)
− 2Kα
The weight ζT1 is given by (3.10). The partial derivative terms C
|ǫ1=0L
(1,ǫ1)
and C2T1 :=
|ǫ1=0L
(2,ǫ1)
are given by:
C1T1 = L
(1,0)
σ1ρ12
( σ3αc3β2
(1 + αc3)
t− (β2 + β3)W 2t
C2T1 = L
(2,0)
−σ2β3W 2t dt,
with C10 = C
0 = 0 and β2 :=
(1+αc2)2
, β3 :=
(1+αc3)2
. The Malliavin covariance
matrix of the vector (L
(1,0)
(2,0)
) is equal to:
(1,0)
(2,0)
(1 + ρ212)(L
(1,0)
)2T1σ
1 2ρ12(L
(1,0)
(2,0)
)T1σ1σ2
2ρ12(L
(1,0)
(2,0)
)T1σ1σ2 (1 + ρ
12)(L
(2,0)
)2T1σ
⇒ det
(1,0)
(2,0)
(1,0)
(2,0)
)2T 21 σ
2(1− ρ212).
12 MARIA SIOPACHA AND JOSEF TEICHMANN
The determinant is not zero as long as ρ12 6= 1, which is a natural assumption.
Hence under this condition, its inverse is given by:
(1,0)
(2,0)
(1− ρ212)
1+ρ212
(1,0)
)2T1σ
− 2ρ12
(1,0)
(2,0)
T1σ1σ2
− 2ρ12
(1,0)
(2,0)
T1σ1σ2
1+ρ212
(2,0)
)2T1σ
Write the weight ζT1 = ζ
+ ζ2T1 , where the first weight ζ
is obtained as:
ζ1T1 =
(1,0)
C1T1γ
11 + C
γ−112
+D1tL
(2,0)
C1T1γ
21 + C
γ−122
δW 1t ,
and ζ2T1 similarly:
ζ2T1 =
(1,0)
C1T1γ
11 + C
γ−112
+D2tL
(2,0)
C1T1γ
21 + C
γ−122
δW 2t .
Performing all necessary calculations, we conclude that:
ζ1T1 = ρ12
W 1T1
(σ3αc3β2T1
2(1 + αc3)
− (β2 + β3)
W 2t dt
ρ12(β2 + β3)T1
− ρ12
(ρ12β3T1
W 2t dt
Analogously we obtain ζ2T1 as:
ζ2T1 = ρ
W 2T1
(σ3αc3β2T1
2(1 + αc3)
(β2 + β3)
W 2t dt
(β2 + β3)T1
(β3T1
W 2t dt
Notice that the weights are functions of normal variables and thus the calculation of
the weak Taylor price amounts just to computation of deterministic integrals. Table
3 gives the swaption prices in bps for parameters N = 3, α = 0.25, σ1 = 18%,
σ2 = 15%, σ3 = 12%, c0 = 5.28875%, c1 = 5.37375%, c2 = 5.40%, c3 = 5.40125%
and ρ12 = 0.75.
strikes K=4% K=4.5% K=4.75% K=5% K=5.15% K=5.25%
benchmark 10.2240 6.5386 4.7454 3.1060 2.2599 1.7758
frozen drift 10.2132 6.5326 4.7419 3.1028 2.2582 1.7618
strong Taylor 10.2240 6.5386 4.7454 3.1060 2.2599 1.7758
weak Taylor 10.2266 6.5407 4.7485 3.1064 2.2593 1.7626
Table 2: Swaption values in bps for parameters ǫ1 = 1, α = 0.25, σ1 = 18%, σ2 = 15%,
σ3 = 12%, c0 = 5.28875%, c1 = 5.37375%, c2 = 5.40%, c3 = 5.40125% and ρ12 = 0.75.
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 13
3.4. The Stochastic Volatility LIBOR Market Model. In this section, we
develop a stochastic volatility LMM. The stochastic volatility parameter vt follows
a square root process, like in the extensively applied Heston model [Hes93]. The
resulting model, called hereafter the stochastic volatility LMM (SVLMM), has the
following dynamics under the terminal measure:
dLit = σ
i(t)Lit
j=i+1
1 + αL
vtdt+ dW
, i = 1, ..., N,(3.13)
dvt = κ(θ − vt)dt+ ǫ2
vtdBt,
where κ, θ, ǫ2 ∈ R+. The Brownian motions Wt = (W 1t , ...,WNt ) and Bt are ex-
pressed under the terminal measure with correlations dW it dBt = ρidt and dW
ρijdt for i, j = 1, ...N . We assume additionally that the filtration (Ft)0≤t≤T is gen-
erated by both Brownian motions. Observe that the process vt is a time-changed
squared Bessel process with dimension δ = 4κθ/ǫ22. If δ ≥ 2, then the point zero is
unattainable. So we require 2κθ ≥ ǫ22 for the process vt not to reach zero.
3.4.1. Pricing a multi-LIBOR option. In this section, we aim at approximating the
price of an option with payoff depending on the vector L
i,k,ǫ1,ǫ2
i,ǫ1,ǫ2
, . . . , L
k,ǫ1,ǫ2
We interpret the volatility of the volatility parameter ǫ2 as a parameter on which
the LIBOR rates depend. Overall, we parametrise the SVLMM by both ǫ1 and ǫ2
and correct prices in a weak sense introducing Malliavin weights.
Proposition 3. Consider the SVLMM (3.13) and assume that the Malliavin co-
variance matrix γ(L
i,k,0,0
) is invertible. Then the price of an option with payoff
i,k,ǫ1,ǫ2
), i ≤ k ≤ N , where ψ is a bounded measurable function, can be approx-
imated by the weak Taylor approximation of order one:
(ǫ1,ǫ2)
i,k,ǫ1,ǫ2
)) = P (0, T )
i,k,0,0
+ ǫ1EP
i,k,0,0
+ ǫ2EP
i,k,0,0
,(3.14)
where the Malliavin weights ζTi , πTi are given by:
ζTi = δ
i,k,0,0
)Tγ−1(L
i,k,0,0
i,k,ǫ1,0
,(3.15)
πTi = δ
i,k,0,0
)Tγ−1(L
i,k,0,0
i,k,0,ǫ2
),(3.16)
for t ≤ Ti.
Proof. The weights ζTi and πTi are obtained by (2.4). We derive (3.14) by noticing
that:
i,k,ǫ1,ǫ2
i,k,0,0
i,k,ǫ1,0
i,k,0,ǫ2
i,k,0,0
+ ǫ1E
i,k,0,0
+ ǫ2E
i,k,0,0
from Definition 2 for n = 1. �
14 MARIA SIOPACHA AND JOSEF TEICHMANN
Example 4. Let N = 2 and consider the SVLMM where the volatility functions
σi(t) : [0, Ti] → R for i = 1, 2 are assumed to be constant and in particular σ1(t) =
σ1, σ
2(t) = σ2. We derive an approximative formula for the price of a payers
swaption with maturity T1 and strike price K. The underlying swap is entered at
T1 and has payment dates T2, T3. Under the terminal measure P we can write the
SDEs for the LIBOR rates and stochastic volatility as:
dvǫ2t = κ
θ − vǫ2t
dt+ ǫ2
vǫ2t dBt,
(1,ǫ1,ǫ2)
t = −L
(1,ǫ1,ǫ2)
t ρ12
(2,ǫ1,ǫ2)
1 + αX
(2,ǫ1,ǫ2)
t dt+ σ1L
(1,ǫ1,ǫ2)
vǫ2t dW
(2,ǫ2)
t = σ2L
(2,ǫ2)
vǫ2t dW
(2,ǫ1)
t = ǫ1
(2,ǫ2)
vǫ2t dW
W 1t and W
t are assumed to be correlated, so correlations are as dW
t dBt = ρidt
and dW 1t dW
t = ρ12 for i = 1, 2. The (0, 0)-model is given by:
v0t = exp (−κt)(v00 − θ) + θ,
(1,0,0)
= c1 exp
v0t dW
( αc2ρ12
1 + αc2
(2,0)
= c2 exp
v0t dW
(2,0,0)
t = c2,
with c :=
v0t dt = θT1 −
(exp (−κT1) − 1). As in the previous example, we
compare the following option prices:
• benchmark price;
• frozen drift;
• weak Taylor price (3.14).
The benchmark price is given by (3.12) with N = 2 and i = 1:
swptn
0 = P (0, T )EP
(1,ǫ1,ǫ2)
(2,ǫ2)
(1,ǫ1,ǫ2)
(2,ǫ2)
−Kα2L(2,ǫ2)
The weak Taylor price is obtained by (3.14):
swptn
0 = P (0, T )
(1,0,0)
(2,0)
+ α2L
(1,0,0)
(2,0)
−Kα2L(2,0)
− 2Kα
+ ǫ1EP
(1,0,0)
(2,0)
+ α2L
(1,0,0)
(2,0)
−Kα2·
· L(2,0)
− 2Kα
+ ǫ2EP
(1,0,0)
(2,0)
+ α2L
(1,0,0)
(2,0)
−Kα2L(2,0)T1 − 2Kα
We calculate the Malliavin weights ζT1 , πT1 as given by (3.15) and (3.16) corre-
spondingly. We can express the weight ζT1 as:
ζT1 = ζ
+ ζ2T1 ,
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 15
with:
ζ1T1 =
(1,0,0)
(1,ǫ1,0)
γ−1(L
(1,0,0)
(2,0,0)
(2,0,0)
(1,ǫ1,0)
γ−1(L
(1,0,0)
(2,0,0)
δW 1t ,
ζ2T1 =
(1,0,0)
(1,ǫ1,0)
γ−1(L
(1,0,0)
(2,0,0)
(2,0,0)
(1,ǫ1,0)
γ−1(L
(1,0,0)
(2,0,0)
δW 2t .
since ∂
|ǫ1=0 L
(2,ǫ1,0)
= 0. The partial derivative term with respect to ǫ1 for L
given by:
(1,ǫ1,0)
(1,0,0)
(1 + αc2)2
v0sdW
= −σ1ρ12β2L(1,0,0)T1
θ(T1 − t)−
v00 − θ
(exp (−κT1)− exp (−κt))
dW 2t ,
where β2 =
(1+αc2)2
. Similarly the weight πT1 is given by:
πT1 = π
+ π2T1 ,
with:
π1T1 =
(l,0,0)
(j,0,ǫ2)
(1,0,0)
(2,0,0)
δW 1t ,
π2T1 =
(l,0,0)
(j,0,ǫ2)
(1,0,0)
(2,0,0)
δW 2t .
Partial derivative terms are equal to:
(1,0,ǫ2)
(1,0,0)
exp (−κt)
exp (κs)
v0sdBsdW
αc2σ1σ2ρ12
1 + αc2
exp (κs)
exp (−κT1)− exp (−κs)
Doing similar calculations, we derive the second partial derivative:
(2,ǫ2)
(2,0)
σ1σ2Vtdt
where Vt = exp (−κt)
exp (κs)
v0sdBs.
16 MARIA SIOPACHA AND JOSEF TEICHMANN
We calculate the Malliavin covariance matrix γ
(1,0,0)
(2,0,0)
and its in-
verse.
(1 + ρ212)(L
(1,0,0)
)2σ21
v0t dt
︸ ︷︷ ︸
2ρ12L
(1,0,0)
(2,0,0)
σ1σ2c
2ρ12L
(1,0,0)
(2,0,0)
σ1σ2c (1 + ρ
12)(L
(2,0,0)
)2σ22c
⇒ det
(1,0,0)
(2,0,0)
(1,0,0)
(2,0,0)
)2σ21σ
2(1− ρ212).
Hence its inverse is given by, for ρ12 6= 1:
γ−1 =
(1− ρ212)
1+ρ212
(1,0,0)
)2σ21c
− 2ρ12
(1,0,0)
(2,0,0)
σ1σ2c
− 2ρ12
(1,0,0)
(2,0,0)
σ1σ2c
1+ρ212
(2,0,0)
)2σ22c
If we define Xi =
v0t dW
t , i = 1, 2 and Y =
θ(T1−t)− v
(exp (−κT1)−
exp (−κt))
dW 2t , we finally obtain the weights as:
ζ1T1 = −
ρ12β2
X1Y − Cov(X1, Y )
ζ2T1 =
ρ212β2
X2Y − Cov(X2, Y )
Moreover, for the weight πT1 we define:
exp (−κT1)− exp (−κt)
and random variables Di, Zi for i = 1, 2:
g(s)dW isdW
g(s)dZisdW
where the Brownian motions Zit are independent from W
t and f(t) =
exp (−κt)√
g(s) = exp (κs)
v0s . Therefore, we obtain the weights as:
π1T1 =
X1(ρ1D1 +
1− ρ21Z1) +
(αc2(2ρ12σ2 + σ1) + σ1
1 + αc2
(αc2(2ρ12σ2 + σ1) + σ1
1 + αc2
X1(ρ2D2 +
1− ρ22Z2)+
σ1BX1
− σ1ρ1E
WEAK AND STRONG TAYLOR METHODS FOR NUMERICAL SOLUTIONS OF SDES 17
where E equals to
1− exp (−κT1)
. Similarly we get π2T1 as:
π2T1 =
X2(ρ2D2 +
1− ρ22Z2) +
σ2X2 + 1
− σ2ρ2E
− ρ12
X2(ρ1D1 +
1− ρ21Z1) +
σ1BX2
(αc2(2ρ12σ2 + σ1) + σ1
1 + αc2
− σ1ρ2E
(αc2(2ρ12σ2 + σ1) + σ1
1 + αc2
In this example, the weights are functions of normal variables and double sto-
chastic integrals, which are computed via simulation. Table 4 reports the swaption
prices in bps with parameters N = 2, α = 1.5, σ1 = 25%, σ2 = 15%, c0 = 5.28875%,
c1 = 5.4%, c2 = 5.39%, v0 = 1, ρ1 = −0.75, ρ2 = −0.6, κ = 2.3767, θ = 0.2143,
ǫ2 = 25%, ρ12 = 0.63.
strikes K=3.5% K=4% K=5% K=6% K=7% K=8%
benchmark 3.8984 2.9221 1.2588 0.3858 0.1019 0.0216
(0, 0)-model 3.8951 2.9053 1.2705 0.3966 0.0942 0.0185
weak Taylor 3.8990 2.9159 1.2694 0.3791 0.1042 0.0210
Table 3: Stochastic volatility swaption values in bps for parameters ǫ1 = 1, α = 1.5, σ1 = 25%,
σ2 = 15%, c0 = 5.28875%, c1 = 5.4%, c2 = 5.39%, v0 = 1, ρ1 = −0.75, ρ2 = −0.6,
κ = 2.3767, θ = 0.2143, ǫ2 = 25%, ρ12 = 0.63.
References
[BGM97] A. Brace, D. Gatarek, and M. Musiela, The Market Model of Interest Rate Dynamics,
Mathematical Finance 7 (1997), no. 2, 127–155.
[BM01] D. Brigo and F. Mercurio, Interest Rate Models: Theory and Practice, Springer Finance,
Springer, 2001.
[BW00] A. Brace and R.S. Womersley, Exact Fit to the Swaption Volatility Matrix using Semi-
definite Programming, Working paper, presented at ICBI Global Derivatives Confer-
ence, Paris, April 2000, 2000.
[EÖ05] E. Eberlein and F. Özkan, The Lévy LIBOR model, Finance and Stochastics 9 (2005),
327348.
[Hes93] S. Heston, A Closed-Form Solution for Options with Stochastic Volatility with Applica-
tions to Bond and Currency Options, The Review of Financial Studies 6 (1993), no. 2,
327–343.
[Jam97] F. Jamshidian, LIBOR and Swap Market Models and Measures, Finance and Stochas-
tics 1 (1997), 293–330.
[KM97] A. Kriegl and P. W. Michor, The Convenient Setting of Global Analysis, American
Mathematical Society, 1997.
[Mal97] P. Malliavin, Stochastic Analysis, Springer, 1997.
[MR98] M. Musiela and M. Rutkowski, Martingale Methods in Financial Modelling, second ed.,
Springer, 1998.
[MSS97] K. Miltersen, K. Sandmann, and D. Sondermann, Closed Form Solutions for Term
Structures Derivatives with Log-Normal Interest Rates, Journal of Finance 52 (1997),
409–430.
[MT06] P. Malliavin and A. Thalmaier, Stochastic Calculus of Variations in Mathematical Fi-
nance, Springer, 2006.
18 MARIA SIOPACHA AND JOSEF TEICHMANN
[Nua06] D. Nualart, The Malliavin Calculus and Related Topics, second ed., Springer Verlag,
2006.
[Sch02] E. Schlögl, A Multicurrency Extension of the Lognormal Interest Rate Market Models,
Finance and Stochastics 6 (2002), 173196.
Department of Mathematical Methods in Economics, Vienna University of Technol-
ogy, Wiedner Hauptstrasse 8–10/105–1, A-1040 Vienna, Austria.
E-mail address: [josef.teichmann,siopacha]@fam.tuwien.ac.at
	1. Introduction and Setting
	2. Weak and strong Taylor methods - Structure Theorems
	3. Applications from Financial Mathematics
	3.1. The LIBOR Market Model
	3.2. Freezing the Drift
	3.3. Correcting the Frozen Drift
	3.4. The Stochastic Volatility LIBOR Market Model
	References
