Weak and Strong Taylor methods for numerical solutions of stochastic differential equations

TL;DR

This paper develops weak and strong Taylor expansion methods for solving stochastic differential equations, providing tractable formulas for LIBOR market models that improve pricing accuracy without complex numerical schemes.

Contribution

It introduces weight expressions for Taylor coefficients of SDE solutions and applies them to LIBOR models, offering a practical alternative to traditional drift freezing methods.

Findings

Accurate pricing formulas for LIBOR models derived from Taylor expansions.

Comparable accuracy to full numerical schemes with simpler expressions.

Numerical examples confirm the effectiveness of the proposed methods.

Abstract

We apply results of Malliavin-Thalmaier-Watanabe for strong and weak Taylor expansions of solutions of perturbed stochastic differential equations (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.

Tables

Table 1. Table 1. Caplet values in bps for parameters ϵ 1 = 1 subscript italic-ϵ 1 1 \epsilon_{1}=1 , α = 0.50137 𝛼 0.50137 \alpha=0.50137 , c 1 = 3.86777 % subscript 𝑐 1 percent 3.86777 c_{1}=3.86777\% , c 2 = 3.7574 % subscript 𝑐 2 percent 3.7574 c_{2}=3.7574\% , c 3 = 3.8631 % subscript 𝑐 3 percent 3.8631 c_{3}=3.8631\% and T 1 = 1.53151 subscript 𝑇 1 1.53151 T_{1}=1.53151 .

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 2. Table 2. Swaption values in bps for parameters ϵ 1 = 1 subscript italic-ϵ 1 1 \epsilon_{1}=1 , α = 0.25 𝛼 0.25 \alpha=0.25 , σ 1 = 18 % subscript 𝜎 1 percent 18 \sigma_{1}=18\% , σ 2 = 15 % subscript 𝜎 2 percent 15 \sigma_{2}=15\% , σ 3 = 12 % subscript 𝜎 3 percent 12 \sigma_{3}=12\% , c 0 = 5.28875 % subscript 𝑐 0 percent 5.28875 c_{0}=5.28875\% , c 1 = 5.37375 % subscript 𝑐 1 percent 5.37375 c_{1}=5.37375\% , c 2 = 5.40 % subscript 𝑐 2 percent 5.40 c_{2}=5.40\% , c 3 = 5.40125 % subscript 𝑐 3 percent 5.40125 c_{3}=5.40125\% and ρ 12 = 0.75 subscript 𝜌 12 0.75 \rho_{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 3. Table 3. Stochastic volatility swaption values in bps for parameters ϵ 1 = 1 subscript italic-ϵ 1 1 \epsilon_{1}=1 , α = 1.5 𝛼 1.5 \alpha=1.5 , σ 1 = 25 % subscript 𝜎 1 percent 25 \sigma_{1}=25\% , σ 2 = 15 % subscript 𝜎 2 percent 15 \sigma_{2}=15\% , c 0 = 5.28875 % subscript 𝑐 0 percent 5.28875 c_{0}=5.28875\% , c 1 = 5.4 % subscript 𝑐 1 percent 5.4 c_{1}=5.4\% , c 2 = 5.39 % subscript 𝑐 2 percent 5.39 c_{2}=5.39\% , v 0 = 1 subscript 𝑣 0 1 v_{0}=1 , ρ 1 = − 0.75 subscript 𝜌 1 0.75 \rho_{1}=-0.75 , ρ 2 = − 0.6 subscript 𝜌 2 0.6 \rho_{2}=-0.6 , κ = 2.3767 𝜅 2.3767 \kappa=2.3767 , θ = 0.2143 𝜃 0.2143 \theta=0.2143 , ϵ 2 = 25 % subscript italic-ϵ 2 percent 25 \epsilon_{2}=25\% , ρ 12 = 0.63 subscript 𝜌 12 0.63 \rho_{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

Equations

$$\mathbf{T}^{n}_{\epsilon}(F_{\epsilon}):=\sum_{i=0}^{n}\frac{{\epsilon}^{i}}{i!}\frac{\partial^{i}}{\partial{\epsilon}^{i}}\Big{|}_{{\epsilon}=0}F_{\epsilon},$$

$$\mathbb{E}\big{(}\lvert F_{\epsilon}-\mathbf{T}^{n}_{\epsilon}(F_{\epsilon})\rvert\big{)}=o({\epsilon}^{n}),$$

$$\mathbb{E}\big{(}\lvert f(F_{\epsilon})-f\big{(}\mathbf{T}^{n}_{\epsilon}(F_{\epsilon})\big{)}\rvert\big{)}\leq K\mathbb{E}\big{(}\lVert F_{\epsilon}-\mathbf{T}^{n}_{\epsilon}(F_{\epsilon})\rVert\big{)}=Ko({\epsilon}^{n}).$$

$${\mathbf{W}}_{ϵ}^n(f,F_{ϵ}):=\sum _{i=0}^n\frac{{ϵ}^i}{i!}\mathbb{E}(f(F_0){\pi }_i),$$

$$\lvert\mathbb{E}\big{(}f(F_{\epsilon})\big{)}-\mathbf{W}^{n}_{\epsilon}(f,F_{\epsilon})\rvert=o({\epsilon}^{n}).$$

$$\mathbb{E}\big{(}f(F_{\epsilon})\big{)}=\sum_{i\geq 0}\frac{{\epsilon}^{i}}{i!}\int_{\mathbb{R}^{N}}f(x)\mu_{i}(dx),$$

$${\int }_{{\mathbb{R}}^N}f(x){\mu }_i(dx)=\mathbb{E}(f(F_0){\pi }_i).$$

$$\frac{d}{d\epsilon}\mathbb{E}\big{(}f(F_{\epsilon})\big{)}=\mathbb{E}\Bigg{(}f(F_{\epsilon})\delta\Big{(}s\mapsto{(D_{s}F_{\epsilon})}^{\mathrm{T}}\gamma^{-1}(F_{\epsilon})\frac{d}{d\epsilon}F_{\epsilon}\Big{)}\Bigg{)}.$$

$$\pi_{1}:=\delta\Big{(}s\mapsto{(D_{s}F_{\epsilon})}^{\mathrm{T}}\gamma^{-1}(F_{\epsilon})\frac{d}{d\epsilon}F_{\epsilon}\Big{)}.$$

$$\frac{d^2}{d{ϵ}^2}\mathbb{E}(f(F_{ϵ}))$$

$$\displaystyle+\mathbb{E}\Bigg{(}f(F_{\epsilon})\delta\Big{(}s\mapsto{(D_{s}\frac{dF_{\epsilon}}{d\epsilon})}^{\mathrm{T}}\gamma^{-1}(F_{\epsilon})\frac{dF_{\epsilon}}{d\epsilon}\Big{)}\Bigg{)}-$$

$$\displaystyle-\mathbb{E}\Bigg{(}f(F_{\epsilon})\delta\Big{(}s\mapsto{(D_{s}F_{\epsilon})}^{\mathrm{T}}\gamma^{-1}(F_{\epsilon})\frac{d\gamma(F_{\epsilon})}{d\epsilon}\gamma^{-1}(F_{\epsilon})\frac{dF_{\epsilon}}{d\epsilon}\Big{)}\Bigg{)}+$$

$$\displaystyle+\mathbb{E}\Bigg{(}f(F_{\epsilon})\delta\Big{(}s\mapsto{(D_{s}F_{\epsilon})}^{\mathrm{T}}\gamma^{-1}(F_{\epsilon})\frac{d^{2}F_{\epsilon}}{d{\epsilon}^{2}}\Big{)}\Bigg{)}.$$

$$a_s$$

$${\pi }_n$$

$${\pi }_0$$

$$\frac{d^{n}f(F_{ϵ})}{d{ϵ}^n}=\sum _{|\alpha |\le n}{f}^{(\alpha )}(F_{ϵ})p_{\alpha },$$

$$\mathbb{E}(\sum _{|\alpha |\le n}{f}^{(\alpha )}(F_{ϵ})p_{\alpha })=E(f(F_{ϵ}){\pi }_n).$$

$$F_{ϵ}=\sum _{i=0}^{\infty }\frac{{ϵ}^i}{i!}{F}^i,$$

$$\lvert\mathbb{E}\big{(}f(F_{\epsilon})\big{)}-\mathbb{E}\big{(}f(F^{0}+\epsilon F^{1})\big{)}\rvert\leq Ko(\epsilon),$$

$$D_s{F}^0=h(s),$$

$$\gamma (F^0)={\int }_0^{\infty }h(s{)}^{2}ds,$$

$$a_s=\frac{h(s)}{{\int }_0^{\infty }h(s{)}^{2}ds}.$$

$$\mathbb{E}\big{(}f(F^{0})\big{)}+\epsilon\mathbb{E}\big{(}f(F^{0})\pi_{1}\big{)},$$

$$\pi_{1}=\delta\big{(}s\mapsto a_{s}F^{1}\big{)}.$$

$$\mathbb{E}\big{(}f(F^{0})\pi_{1}\big{)}=\mathbb{E}\big{(}f(F^{0})F^{1}\int_{0}^{\infty}a_{s}dW_{s}\big{)}-\int_{0}^{\infty}\mathbb{E}\big{(}f(F^{0})D_{s}F^{1}a_{s}\big{)}ds.$$

$$\mathbb{E}\big{(}f(F^{0})\big{)}+\epsilon\mathbb{E}\big{(}f(F^{0})\pi_{1}\big{)}+\epsilon^{2}\mathbb{E}\big{(}f(F^{0})\pi_{2}\big{)},$$

$$\pi_{2}=\delta\big{(}s\mapsto\pi_{1}a_{s}F^{1}\big{)}+\delta\big{(}s\mapsto a_{s}F^{2}\big{)}.$$

$$L_{t}^{i}=L_{t}(T_{i},T_{i+1})=\frac{1}{\alpha}\bigg{(}\frac{P(t,T_{i})}{P(t,T_{i+1})}-1\bigg{)}.$$

$$dL_{t}^{i}=\sigma^{i}(t)L_{t}^{i}\Big{(}-\sum_{j=i+1}^{N}\frac{\alpha L_{t}^{j}\sigma^{j}(t)}{1+\alpha L_{t}^{j}}\rho_{ij}\Big{)}dt+\sigma^{i}(t)L_{t}^{i}dW_{t}^{i},i=1,...,N.$$

$$\frac{\alpha L_t^j}{1+\alpha L_t^j}\approx \frac{\alpha L_0^j}{1+\alpha L_0^j}.$$

$$d{X}_t^{(i,{ϵ}_1)}$$

$$\mathbf{T}_{\epsilon_{1}}^{1}(X_{t}^{(i,\epsilon_{1})})=X_{t}^{(i,0)}+\epsilon_{1}\frac{\partial}{\partial\epsilon_{1}}\Big{|}_{\epsilon_{1}=0}X_{t}^{(i,\epsilon_{1})}.$$

$${\mathbf{T}}_{{ϵ}_1}^1(X_t^{(i,{ϵ}_1)})\simeq X_t^{(i,{ϵ}_1)}=X_0^{(i,0)}+{ϵ}_1{Y}_t^i+o({ϵ}_1),$$

Full text

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 (ÖBFA) through CDG

1. Introduction and Setting

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space carrying an $N$-dimensional Brownian motion ${(W_{t})}_{t\geq 0}$ with a $d\times d$ correlation matrix. We consider smooth curves $F_{\epsilon}:\mathbb{R}\to L^{2}(\Omega;\mathbb{R}^{N})$ of random variables, where $\epsilon\in\mathbb{R}$ is a parameter. We apply Taylor theorems to obtain strong approximations of the curve $F_{\epsilon}$ at $\epsilon=0$ and we apply partial integration on Wiener space to obtain weak approximations of the law of $F_{\epsilon}$ for small values of $\epsilon$.

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 $\epsilon$ if – for instance – the payoff $f:\mathbb{R}^{N}\to\mathbb{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 $\epsilon\mapsto F_{\epsilon}$, where $\epsilon\in\mathbb{R}$ and $F_{\epsilon}\in L^{2}(\Omega;\mathbb{R}^{N})$.

Definition 1**.**

A strong Taylor approximation of order $n\geq 0$ is a (truncated) power series

[TABLE]

such that

[TABLE]

holds true as $\epsilon\to 0$.

Remark 1**.**

In our setting a strong Taylor approximation of any order $n\geq 0$ of the curve $F_{\epsilon}$ can always be obtained, see for instance [KM97].

Let $f:\mathbb{R}^{N}\to\mathbb{R}$ be a Lipschitz function with Lipschitz constant $K$, then we obtain

[TABLE]

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\geq 0$ is a power series for each bounded, measurable $f:\mathbb{R}^{N}\to\mathbb{R}$,

[TABLE]

where $\pi_{i}\in L^{1}(\Omega)$ denote real valued, integrable random variables, such that

[TABLE]

Remark 2**.**

The weights $\pi_{i}$ for $i\geq 1$ are called Malliavin weights.

Remark 3**.**

If the law of $F_{\epsilon}$ is real analytic at $\epsilon=0$ in the weak sense, i.e.  if there exist (signed) measures $\mu_{i}$ such that for all bounded, measurable $f:\mathbb{R}^{N}\to\mathbb{R}$ the following series converges and the equality

[TABLE]

holds true, precisely then we do have a converging weak Taylor expansion. We aim for constructing stochastic representations of the following type, for $i\geq 0$:

[TABLE]

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 $\mathcal{D}^{\infty}(\mathbb{R}^{N})$ see [Mal97] or [Nua06].

Theorem 1**.**

Let $F_{\epsilon}:\mathbb{R}\to\mathcal{D}^{\infty}(\mathbb{R}^{N})$ be smooth and assume that the Malliavin covariance matrix $\gamma(F_{\epsilon})$ is invertible with $p$-integrable inverse for every $p\geq 1$ around $\epsilon=0$ (i.e. on an open interval containing $\epsilon=0$). Then there is a weak Taylor approximation of any order $n\geq 0$ and there are explicit formulas for the weights $\pi_{i}$. If we only know that the Malliavin covariance matrix $\gamma(F_{0})$ is invertible with $p$-integrable inverse, then we can also calculate the Malliavin weights, since they depend only on $\gamma(F_{0})$.

Proof.

Fix $n\geq 0$ and take a smooth test function $f:\mathbb{R}^{N}\to\mathbb{R}$ and assume that $\gamma^{-1}(F_{\epsilon})$ exists as a smooth curve in $\mathcal{D}^{\infty}$ on a open $\epsilon$-interval containing $\epsilon=0$. By standard arguments we can prove the following formula

[TABLE]

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 $\epsilon$-dependence of the Skorohod integral is smooth due to basic properties of $\mathcal{D}^{\infty}$. Hence, we can calculate higher derivatives of the left hand side by iterating the above procedure and differentiating the Skorohod integral. We denote

[TABLE]

We write then, pars pro toto, the formula for the second derivative

[TABLE]

This formula makes perfect sense at $\epsilon=0$ and – by induction – we see that we can perform this step for any derivative. The general, recursive result is the following:

[TABLE]

Here we understand the weights $\pi_{n}$ as $\epsilon$-dependent, whereas in the final formulas we put $\epsilon=0$. This proves the result for smooth test functions $f$ and under the assumption that the Malliavin covariance matrix is invertible around $\epsilon=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 Faà-di-Bruno-formula we obtain

[TABLE]

where $p_{\alpha}$ is a well-defined polynomial in derivatives of the curve $\epsilon\mapsto F_{\epsilon}$, for a multi-index $\alpha$. Since $\mathcal{D}^{\infty}$ is an algebra, see [Mal97], the above expression lies in $L^{p}(\Omega)$ for each $p\geq 0$. The previous result provides a representation of the partial integration result for

[TABLE]

The structure of the weights is seen from above. The result can be considered as a dual version of the Faà-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 approximation. The method works in order to replace time-consuming iteration schemes, like the Euler-scheme, by simulations of “simple” Itô integrals.

Example 1**.**

We deal with a generic, real-valued random variable over a one-dimensional Gaussian space, see [Nua06], i.e.

[TABLE]

where the $F^{i}$ lie in the $(i+1)^{st}$ Wiener chaos $\mathcal{H}_{i+1}(\Omega)$ (one can think of a Hermite expansion for instance) and the sum is understood in the $L^{2}$-sense. From the strong expansion we obtain immediately – for a given Lipschitz function $f:\mathbb{R}\to\mathbb{R}$ – that

[TABLE]

as $\epsilon\to 0$, where $K$ denotes the Lipschitz constant of $f$. This simple approximation can be sometimes quite useful.

We assume now that $F^{0}=\int_{0}^{\infty}h(s)dW_{s}$ has non-vanishing variance in order to calculate the weights, which do depend only on $\gamma(F_{0})$. The strong Taylor approximation is given by definition, the weak Taylor expansion can be constructed by the previous recursive formulas and the specifications

[TABLE]

In order to obtain a first-order approximation for bounded, measurable random variables we therefore have to calculate

[TABLE]

where

[TABLE]

This amounts to an integration of $f$ times a polynomial with respect to a Gaussian density, since:

[TABLE]

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.

The second-order weak Taylor approximation is given by

[TABLE]

where

[TABLE]

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 $\frac{dF_{\epsilon}}{d\epsilon}|_{\epsilon=0}$ is of simple form or – even more important – where the Malliavin covariance matrix $\gamma(F_{0})$ 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 $(\Omega,\mathcal{F}_{T},\mathbb{P},(\mathcal{F}_{t})_{0\leq t\leq T})$ be a complete probability space, supporting an $N$-dimensional Brownian motion $W_{t}=(W_{t}^{1},...W_{t}^{N})_{0\leq t\leq T}$. The factors are correlated with $dW_{t}^{i}$ $dW_{t}^{j}=\rho_{ij}dt$. Let $0=T_{0}<T_{1}<T_{2}<\ldots<T_{N}<T_{N+1}=:T$ be a discrete tenor structure and $\alpha:=T_{i+1}-T_{i}$ the accrual factor for the time period $[T_{i},T_{i+1}]$, $i=0,\ldots,N$. Let $P(t,T_{i})$ denote the value at time $t$ of a zero coupon bond with maturity $T_{i}\in[0,T]$. The measure $\mathbb{P}$ is the terminal forward measure, which corresponds to taking the final bond $P(t,T)$ as numéraire. The forward LIBOR rate $L_{t}^{i}:=L_{t}(T_{i},T_{i+1})$ at time $t\leq T_{i}$ for the period $[T_{i},T_{i+1}]$ is given by:

[TABLE]

We assume that for any maturity $T_{i}$ there exists a bounded, continuous, deterministic function $\sigma^{i}(t):[0,T_{i}]\to\mathbb{R}$, which represents the volatility of the LIBOR $L_{t}^{i}$, $i=1,...,N$. The log normal LIBOR market model can be expressed under the measure $\mathbb{P}$ as:

[TABLE]

3.2. Freezing the Drift

The dynamics of forward LIBORs for $i=1,...,N-1$ depend on the stochastic drift term $\frac{\alpha L_{t}^{j}}{1+\alpha L_{t}^{j}}$, $i\leq j\leq N$, which is determined by LIBOR 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. In addition, it encumbers the numerical implementation of the model. Common practice is to approximate this term by its starting value $\frac{\alpha L_{0}^{j}}{1+\alpha L_{0}^{j}}$ or as it is widely referred to as freezing the drift, i.e.

[TABLE]

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\mapsto\frac{\alpha x_{+}}{1+\alpha x_{+}}$ with small Lipschitz constant $\alpha$.

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 Lévy LIBOR model by Eberlein and Özkan [EÖ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 parameters.

Lemma 1**.**

Let $\epsilon_{1}\in\mathbb{R}$ and consider for $i=1,\ldots,N$ the following stochastic differential equation:

[TABLE]

defined on the complete probability space $(\Omega,\mathcal{F}_{T},\mathbb{P},(\mathcal{F}_{t})_{0\leq t\leq T})$ where $W_{t}$ is an $N$-dimensional Brownian motion under the measure $\mathbb{P}$ with $dW_{t}^{i}dW_{t}^{j}=\rho_{ij}dt$. Then the first-order strong Taylor approximation for $X_{t}^{(i,\epsilon_{1})}$ is given by:

[TABLE]

Proof.

By (1) we obtain for $n=1$:

[TABLE]

since $X_{t}^{(i,0)}=X_{0}^{(i,0)}$ and where $Y^{i}_{t}:=\frac{\partial}{\partial\epsilon_{1}}|_{\epsilon_{1}=0}X_{t}^{(i,\epsilon_{1})}$ is the first-order correction term. By differentiating (3.2) with respect to $\epsilon_{1}$, we calculate:

[TABLE]

and derive $Y^{i}_{t}$ as the solution to the above linear SDE:

[TABLE]

with $Y^{i}_{0}=0$. ∎

Remark 6**.**

We parametrise the LIBOR market model in terms of the parameter $\epsilon_{1}$ as follows:

[TABLE]

and assume at $t=0$ that $L_{0}^{(i,\epsilon_{1})}=X_{0}^{(i,\epsilon_{1})}$ for all $\epsilon_{1}$ and all $i=1,...,N$. If $\epsilon_{1}=1$, what we obtain is the standard LIBOR market model formulation and in particular $L_{t}^{(i,1)}=X_{t}^{(i,1)}$. For $\epsilon_{1}=0$, $X_{t}^{(i,0)}$ equals its starting value and thus the drift term in the following SDE is no longer stochastic:

[TABLE]

The next proposition provides a way for a pathwise approximation of $L_{t}^{(i,\epsilon_{1})}$, by means of adjusting its SDE. This is achieved by adding $\mathbf{T}_{\epsilon_{1}}^{n}(X_{t}^{(j,\epsilon_{1})})$ in the frozen drift part.

Proposition 1**.**

Assume the setup of Lemma 1 and assume further at $t=0$ that $L_{0}^{(i,\epsilon_{1})}=X_{0}^{(i,\epsilon_{1})}$ for all $\epsilon_{1}$ and all $i=1,...,N$. Then the stochastic differential equation for $L_{t}^{(i,\epsilon_{1})}$ with the unfrozen drift:

[TABLE]

can be strongly approximated as $\epsilon_{1}\downarrow 0$ by

[TABLE]

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_{t}^{(j,\epsilon_{1})}$ with $(X_{t}^{(j,\epsilon_{1})})_{+}$ in (3.6) to obtain:

[TABLE]

This yields no change for the dynamics of $L_{t}^{(i,\epsilon_{1})}$, since $X_{t}^{(j,\epsilon_{1})}=(X_{t}^{(j,\epsilon_{1})})_{+}$.

By Taylor’s expansion, we know that as $\epsilon_{1}\downarrow 0$, $\hat{L}_{t}^{(i,\epsilon_{1})}\rightarrow L_{t}^{(i,\epsilon_{1})}$ $\mathbb{P}$-a.s. The estimate for the error term is given by

[TABLE]

∎

Remark 8**.**

The SDE for the approximated $\hat{L}_{t}^{(i,\epsilon_{1})}$ is easier and faster to simulate than (3.1), as it is exhibited by the following example. Notice additionally that $\hat{L}_{t}^{(i,\epsilon_{1})}$ is a continuous functional of the process $Y^{j}_{t}$ (3.4) and of the Brownian path $W_{t}^{i}$. Eventually, by using $\hat{L}_{t}^{(i,\epsilon_{1})}$ 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 $L^{1}$ with strike $K$. Its price is given by:

[TABLE]

Assume that the volatility functions $\sigma^{i}(t):[0,T_{i}]\rightarrow\mathbb{R}$ for $i=1,2,3$ are given by (cf. Brigo and Mercurio [BM01], formulation (6.12)):

[TABLE]

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 $\mathbb{P}$ as:

[TABLE]

with initial values $L_{0}^{(i,\epsilon_{1})}=X_{0}^{(i,\epsilon_{1})}=c_{i}$, for $i=1,2,3$ and for all $\epsilon_{1}$. The Brownian motion vector $(W_{t}^{1},W_{t}^{2},W_{t}^{3})$ is correlated with correlation coefficient $\rho_{ij}$ given by:

[TABLE]

The SDEs for the approximated LIBOR rates $\hat{L}_{t}^{(1,\epsilon_{1})}$ and $\hat{L}_{t}^{(2,\epsilon_{1})}$ are given by:

[TABLE]

The partial derivative terms $Y_{t}^{2}$ and $Y_{t}^{3}$ are equal to:

[TABLE]

We compare three caplet prices:

• •

benchmark price, underlying $L_{t}^{(1,\epsilon_{1})}$;

• •

strong Taylor price, underlying $\hat{L}_{t}^{(1,\epsilon_{1})}$;

• •

frozen drift price, underlying $L_{t}^{(1,0)}$.

Numerical results in basis points (bps) are displayed in Table 1 for parameters $\epsilon_{1}=1$, $N=3$, $\alpha=0.50137$, $c_{1}=3.86777\%$, $c_{2}=3.7574\%$, $c_{3}=3.8631\%$, $T_{1}=1.53151$, $T_{i}=T_{1}+i\alpha$, $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.

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 $\mathbf{L}_{T_{i}}^{i,k,\epsilon_{1}}$ denote the vector of the LIBOR rates $(L_{T_{i}}^{(i,\epsilon_{1})},\ldots,L_{T_{i}}^{(k,\epsilon_{1})})$.

Proposition 2**.**

Assume the setup of Lemma 1, where the $i^{th}$ LIBOR rate is given by:

[TABLE]

with $L_{0}^{(i,\epsilon_{1})}=X_{0}^{(i,\epsilon_{1})}$ for all $\epsilon_{1}$ and all $i=1,...,N$. Assume furthermore that the Malliavin covariance matrix $\gamma(\mathbf{L}_{T_{i}}^{i,k,0})$ is invertible. Then the price of an option with payoff $g(\mathbf{L}_{T_{i}}^{i,k,\epsilon_{1}})$, for $i\leq k\leq N$ and $g$ bounded measurable, can be approximated by the weak Taylor approximation of order one:

[TABLE]

where the Malliavin weight $\zeta_{T_{i}}$ is given by:

[TABLE]

for $t\leq T_{i}$.

Proof.

The weight $\zeta_{T_{i}}$ is obtained by (2.4). Notice that we can write:

[TABLE]

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 $T_{1}$, where the underlying swap is entered at $T_{1}$ and has payment dates $T_{2}$ and $T_{3}$. We assume that the volatility functions $\sigma^{i}(t):[0,T_{i}]\rightarrow\mathbb{R}$ for $i=1,2,3$ are constant:

[TABLE]

such that we obtain under the terminal measure $\mathbb{P}$:

[TABLE]

with initial values $L_{0}^{(i,\epsilon_{1})}=X_{0}^{(i,\epsilon_{1})}=c_{i}$, for $i=1,2,3$ and for all $\epsilon_{1}$. $W_{t}^{1}$ and $W_{t}^{2}$ are correlated with correlation coefficient $\rho_{12}$. We freeze the drifts in the above equations to obtain:

[TABLE]

Similarly to the previous example, we compare four option prices:

• •

benchmark price;

• •

frozen drift;

• •

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 $T_{i}$ can be found for example in [MR98]:

[TABLE]

if the underlying swap is entered at time $T_{i}$ and has payment dates $T_{i+1},...,$ $T$. $\alpha_{k}$ is given by:

[TABLE]

The payers swaption value at time $t=0$ can be written as:

[TABLE]

where $\alpha_{i}:=-1$ and $\mathbb{P}^{{i}}$ denotes the forward measure corresponding to the bond $P(t,T_{i})$ as numéraire. Therefore, its benchmark price is given by the above formula with $N=2$ and $i=1$:

[TABLE]

Its weak Taylor price is given by (3.9) with $i=1$ and $k=N=2$.

[TABLE]

The weight $\zeta_{T_{1}}$ is given by (3.10). The partial derivative terms $C_{T_{1}}^{1}:=\frac{\partial}{\partial\epsilon_{1}}|_{\epsilon_{1}=0}L_{T_{1}}^{(1,\epsilon_{1})}$ and $C_{T_{1}}^{2}:=\frac{\partial}{\partial\epsilon_{1}}|_{\epsilon_{1}=0}L_{T_{1}}^{(2,\epsilon_{1})}$ are given by:

[TABLE]

and:

[TABLE]

with $C_{0}^{1}=C_{0}^{2}=0$ and $\beta_{2}:=\frac{\alpha c_{2}^{2}\sigma_{2}^{2}}{(1+\alpha c_{2})^{2}}$, $\beta_{3}:=\frac{\alpha c_{3}^{2}\sigma_{3}^{2}}{(1+\alpha c_{3})^{2}}$. The Malliavin covariance matrix of the vector $(L_{T_{1}}^{(1,0)},L_{T_{1}}^{(2,0)})$ is equal to:

[TABLE]

The determinant is not zero as long as $\rho_{12}\neq 1$, which is a natural assumption. Hence under this condition, its inverse is given by:

[TABLE]

Write the weight $\zeta_{T_{1}}=\zeta^{1}_{T_{1}}+\zeta^{2}_{T_{1}}$, where the first weight $\zeta^{1}_{T_{1}}$ is obtained as:

[TABLE]

and $\zeta_{T_{1}}^{2}$ similarly:

[TABLE]

Performing all necessary calculations, we conclude that:

[TABLE]

Analogously we obtain $\zeta_{T_{1}}^{2}$ as:

[TABLE]

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 2 gives the swaption prices in bps for parameters $N=3$, $\alpha=0.25$, $\sigma_{1}=18\%$, $\sigma_{2}=15\%$, $\sigma_{3}=12\%$, $c_{0}=5.28875\%$, $c_{1}=5.37375\%$, $c_{2}=5.40\%$, $c_{3}=5.40125\%$ and $\rho_{12}=0.75$.

3.4. The Stochastic Volatility LIBOR Market

Model

In this section, we develop a stochastic volatility LMM. The stochastic volatility parameter $v_{t}$ 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:

[TABLE]

where $\kappa,\theta,\epsilon_{2}\in\mathbb{R}_{+}$. The Brownian motions $W_{t}=(W_{t}^{1},...,W_{t}^{N})$ and $B_{t}$ are expressed under the terminal measure with correlations $dW^{i}_{t}dB_{t}=\rho_{i}dt$ and $dW^{i}_{t}dW_{t}^{j}=\rho_{ij}dt$ for $i,j=1,...N$. We assume additionally that the filtration $(\mathcal{F}_{t})_{0\leq t\leq T}$ is generated by both Brownian motions. Observe that the process $v_{t}$ is a time-changed squared Bessel process with dimension $\delta=4\kappa\theta/\epsilon_{2}^{2}$. If $\delta\geq 2$, then the point zero is unattainable. So we require $2\kappa\theta\geq\epsilon_{2}^{2}$ for the process $v_{t}$ 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 $\mathbf{L}_{T_{i}}^{i,k,\epsilon_{1},\epsilon_{2}}=(L_{T_{i}}^{i,\epsilon_{1},\epsilon_{2}},\ldots,L_{T_{i}}^{k,\epsilon_{1},\epsilon_{2}})$. We interpret the volatility of the volatility parameter $\epsilon_{2}$ as a parameter on which the LIBOR rates depend. Overall, we parametrise the SVLMM by both $\epsilon_{1}$ and $\epsilon_{2}$ and correct prices in a weak sense introducing Malliavin weights.

Proposition 3**.**

Consider the SVLMM (3.13) and assume that the Malliavin covariance matrix $\gamma(\mathbf{L}_{T_{i}}^{i,k,0,0})$ is invertible. Then the price of an option with payoff $\psi(\mathbf{L}_{T_{i}}^{i,k,\epsilon_{1},\epsilon_{2}})$, $i\leq k\leq N$, where $\psi$ is a bounded measurable function, can be approximated by the weak Taylor approximation of order one:

[TABLE]

where the Malliavin weights $\zeta_{T_{i}},\pi_{T_{i}}$ are given by:

[TABLE]

for $t\leq T_{i}$.

Proof.

The weights $\zeta_{T_{i}}$ and $\pi_{T_{i}}$ are obtained by (2.4). We derive (3) by noticing that:

[TABLE]

from Definition 2 for $n=1$. ∎

Example 4**.**

Let $N=2$ and consider the SVLMM where the volatility functions $\sigma^{i}(t):[0,T_{i}]\rightarrow\mathbb{R}$ for $i=1,2$ are assumed to be constant and in particular $\sigma^{1}(t)=\sigma_{1}$, $\sigma^{2}(t)=\sigma_{2}$. We derive an approximative formula for the price of a payers swaption with maturity $T_{1}$ and strike price $K$. The underlying swap is entered at $T_{1}$ and has payment dates $T_{2},T_{3}$. Under the terminal measure $\mathbb{P}$ we can write the SDEs for the LIBOR rates and stochastic volatility as:

[TABLE]

$W_{t}^{1}$* and $W_{t}^{2}$ are assumed to be correlated, so correlations are as $dW_{t}^{i}dB_{t}=\rho_{i}dt$ and $dW_{t}^{1}dW_{t}^{2}=\rho_{12}$ for $i=1,2$. The $(0,0)$-model is given by:*

[TABLE]

with $c:=\int_{0}^{T_{1}}v_{t}^{0}dt=\theta T_{1}-\frac{v_{0}^{0}-\theta}{\kappa}(\exp{(-\kappa T_{1})}-1)$. As in the previous example, we compare the following option prices:

• •

benchmark price;

• •

frozen drift;

• •

weak Taylor price (3).

The benchmark price is given by (3) with $N=2$ and $i=1$:

[TABLE]

The weak Taylor price is obtained by (3):

[TABLE]

We calculate the Malliavin weights $\zeta_{T_{1}}$, $\pi_{T_{1}}$ as given by (3.15) and (3.16) correspondingly. We can express the weight $\zeta_{T_{1}}$ as:

[TABLE]

with:

[TABLE]

and:

[TABLE]

since $\frac{\partial}{\partial\epsilon_{1}}|_{\epsilon_{1}=0}$ $L_{T_{1}}^{(2,\epsilon_{1},0)}=0$. The partial derivative term with respect to $\epsilon_{1}$ for $L^{1}$ is given by:

[TABLE]

where $\beta_{2}=\frac{\alpha\sigma_{2}^{2}c_{2}^{2}}{(1+\alpha c_{2})^{2}}$. Similarly the weight $\pi_{T_{1}}$ is given by:

[TABLE]

with:

[TABLE]

and:

[TABLE]

Partial derivative terms are equal to:

[TABLE]

Doing similar calculations, we derive the second partial derivative:

[TABLE]

where $V_{t}=\exp{(-\kappa t)}\int_{0}^{t}\exp{(\kappa s)}\sqrt{v_{s}^{0}}dB_{s}$.

We calculate the Malliavin covariance matrix \gamma\big{(}(L_{T_{1}}^{(1,0,0)},L_{T_{1}}^{(2,0,0)})\big{)} and its inverse.

[TABLE]

Hence its inverse is given by, for $\rho_{12}\neq 1$:

[TABLE]

If we define $X_{i}=\int_{0}^{T_{1}}\sqrt{v_{t}^{0}}dW_{t}^{i}$, $i=1,2$ and Y=\int_{0}^{T_{1}}\sqrt{v_{t}^{0}}\Big{(}\theta(T_{1}-t)-\frac{v_{0}^{0}-\theta}{\kappa}(\exp{(-\kappa T_{1})}-\exp{(-\kappa t)})\Big{)}dW_{t}^{2}, we finally obtain the weights as:

[TABLE]

and:

[TABLE]

Moreover, for the weight $\pi_{T_{1}}$ we define:

[TABLE]

and random variables $D_{i}$, $Z_{i}$ for $i=1,2$:

[TABLE]

where the Brownian motions $Z_{t}^{i}$ are independent from $W_{t}^{i}$ and $f(t)=\frac{\exp{(-\kappa t)}}{\sqrt{v_{t}^{0}}}$, $g(s)=\exp{(\kappa s)}\sqrt{v_{s}^{0}}$. Therefore, we obtain the weights as:

[TABLE]

where $E$ equals to \Big{(}\frac{1}{\kappa}\big{(}1-\exp{(-\kappa T_{1})}\big{)}-T_{1}\Big{)}. Similarly we get $\pi_{T_{1}}^{2}$ as:

[TABLE]

In this example, the weights are functions of normal variables and double stochastic integrals, which are computed via simulation. Table 3 reports the swaption prices in bps with parameters $N=2$, $\alpha=1.5$, $\sigma_{1}=25\%$, $\sigma_{2}=15\%$, $c_{0}=5.28875\%$, $c_{1}=5.4\%$, $c_{2}=5.39\%$, $v_{0}=1$, $\rho_{1}=-0.75$, $\rho_{2}=-0.6$, $\kappa=2.3767$, $\theta=0.2143$, $\epsilon_{2}=25\%$, $\rho_{12}=0.63$.