A numerical method to price European derivatives
based on the one factor LIBOR Market Model of interest rates
Luca Vincenzo Ballestra, Graziella Pacelli, Francesco Zirilli
7. Contact us
We consider the problem of pricing European interest rate derivatives based on the LIBOR Market Model (LMM) with one driving factor. We derive a closed-form approximation of the transition probability density functions associated to the stochastic dynamical systems that describe the behaviour of the forward LIBOR interest rates in the LMM. These approximate formulae are based on a truncated power series expansion of the solutions of the Fokker-Planck equations associated to the LMM. The approximate probability density functions obtained are used to price European interest rate derivatives using the method of discounted expectations. The resulting integrals are low-dimensional when the most commonly traded European interest rate derivatives are considered, and they can be computed efficiently using elementary numerical quadrature schemes (i.e. Simpson's rule). The algorithm obtained is very well suited for parallel computing and is tested on the problem of pricing several derivatives including an European swaption and an interest rate spread option. In both cases, the numerical method proposed appears to be accurate (i.e.: relative error of order 10-2, 10-3, or even 10-4) and approximatively between 278 and 63000 times faster than previous methods based on the Monte Carlo simulation of the LMM stochastic dynamical systems. This website distributes to the potential users the software LIBOR_SWAPTION, a Fortran code that implements the numerical method proposed to compute the price of European swaptions.
Given the current time value T0 = 0, we define the so-called “tenor structure” as a finite set of time values Ti, i = 1, 2, … ,n+1, such that T0 < T1 < T2 < … Tn+1. Moreover we define the so-called “years fractions” δi = Ti+1 - Ti , i = 1, 2, … , n.
For i = 1, 2, … , n+1, let us consider a bond which matures at time Ti ,whose price at time t, T0 ≤ t ≤ Ti , is denoted by
Pi(t). The forward LIBOR rate Li(t) is defined as:
The LIBOR Market Model (LMM), introduced by Brace, Gatarek and Musiela  , Miltersen, Sandmann, Sondermann  , and Jamshidian , is the following set of systems of stochastic differential equations:
with initial conditions:
where Wq+1(t) is a p-dimensional standard Wiener process, dWq+1(t) is its stochastic differential, is a p-dimensional vector valued function of time, that models the volatility of the LIBOR interest rate Li(t), i = 1, 2, … , n, · dWq+1 denotes the Euclidean scalar product of and dWq+1, Li,0 is a positive real number and is the value assumed by the random variable Li(T0) with probability one, i = 1, 2, … , n, q +1 is an integer such that 0 ≤ q ≤ n that indicates the bond Pq+1(t) that has been chosen as numeraire asset (see  ,  ), and
We consider the LMM model with one driving factor, that is the stochastic dynamical systems (2)-(3) written when p=1. In particular we derive an approximation of the probability density functions associated to the model (2)-(3) when p=1. When p=1 we write respectively dWq+1(t) and instead than dWq+1(t) and , i=1,2,…,n.
Let us consider the change of variables:
Using the Ito’s lemma and setting p=1, the set of stochastic differential equations (2)-(3) can be rewritten as:
As announced previously we derive a closed-form approximation of the probability density functions associated to the model (6)-(8)  .
Let h be a positive integer such that h≤Min(n,q+1), Zh(t) = (Zh(t), Zh+1(t),…, Zn(t)), zh(t) = ( zh(t), zh+1(t),…, zn(t)),
z´h(t) = (z´h(t), z´h+1(t),…, z´n(t)), τ = t-Th-1. Note that in the previous sentence the variable τ should depend on h. We omit this dependence to keep the notation simple.
Let ph(τ, zh , z´h), 0< τ <Th-Th-1 , denote the probability density function of having Zh(t) = z´h given Zh(Th-1) = zh ,
Th-1< t ≤Th, h = 1,2,…,Min(n,q+1). The probability density functions ph(τ, zh , z´h), h = 1,2,…, Min(n,q+1), are the solutions of the Fokker-Planck equations associated to the stochastic systems of differential equations (6)-(8) with the necessary initial and boundary conditions.
Let ηj,h be a suitable positive constant, we consider the problem (6)-(8) when we substitute to νj(τ) the function
ηj,h + εj,h (νj(τ) - ηj,h), where εj,h is a real parameter, j=1,2,…,n, h=1, 2,…, Min(n,q+1). We note that when εj,h =0, j=1,2,…,n, h=1,2,…,Min(n,q+1), we have an approximate problem with time independent volatilities and when εj,h =1, j=1,2,…,n, h=1,2,…,Min(n,q+1), we have the original problem (6)-(8).
We approximate the probability density function ph(τ, zh , z´h) using a perturbative power series in the parameters δi, εj,h, i = h,h+1,…,n, j = h,h+1,…,n, h=1,2,…, Min(n,q+1) and we truncate this perturbation series at first-order. That is we look for the coefficients of the asymptotic expansion:
In  we have obtained simple closed-form expressions involving elementary functions and the Dirac delta function and its derivatives for p0h(τ, zh, z´h), pδ,1h,i(τ, zh, z´h), pε,1h,j(τ, zh, z´h), i=1,2,…,n, j=1,2,…,n, h=1,2,…,Min(n,q+1). Moreover, due to the fact that the LMM has been considered with p=1, the elementary functions appearing in these closed-form expressions depend only on one variable, and therefore can be used efficiently to price European interest rate derivatives using the method of discounted expectations, since the integrals that solve the pricing problems are low-dimensional.
Let us consider the case of an European payer swaption, whose payoff, paid at time T1, is given by the following function:
and K>0 is the strike price.
The swaption price at time T0, that we denote with Vsw(T0), is computed using the method of discounted expectations, that is:
Due to the fact that the LMM model has been considered with one driving factor (p = 1 in (6)), the evaluation of the integral appearing in (12) leads to the evaluation of a one-dimensional integral, that can be done using an elementary quadrature rule.
We consider a (5 x10) European payer swaption based on the 6-month forward LIBOR rate, that is a swaption maturing at time T1 = 5 year on a swap spanning the tenor structure T1, T2, …, Tn+1 of length 10 year, where n = 20, Ti+1-Ti = 0.5 year, i = 1,2,…,n.
The volatility of the forward LIBOR rate Li(t) is chosen as suggested in  , that is:
where a = 0.2934 year-3/2, b = 1.251 year-1, c = 0.1315, d = 0 year-1/2. Moreover we consider three different cases: Li,0= 0.06 year-1, i = 1,2,…,n, and K= 0.04 year-1 (Test Case A), Li,0= 0.05 year-1, i = 1,2,…,n, and K= 0.05 year-1 (Test Case B), Li,0= 0.04 year-1, i = 1,2,…,n, and K= 0.06 year-1 (Test Case C), corresponding to a swaption in the money, at the money, and out of the money, respectively.
Let us begin with Test Case A. We have computed the price of the swaption at time T0 = 0 using formula (12). The discount factor P1(T0) in (12) has been chosen equal to 1. The (one-dimensional) integral contained in (12) is evaluated using Simpson's rule with 400 subdomains. This choice suffices to keep the relative error due to the numerical integration in formula (12) to be approximately of order 10-4. The result obtained using the analytical approximation (12), that we denote with Vsw,an , is Vsw,an = 0.1538. The computer time needed to evaluate formula (12) with the choices described previously using a computer with a Pentium 4 Processor 1700 MHZ 256 MB Ram is 0.17 second.
In order to test the accuracy of this
approximation an accurate estimation of the exact option value, that we denote
by Vsw,true , has been obtained with a
We have performed 100.000.000 simulations of the payoff function (10) where the Euler-Maruyama scheme has been applied with 500 time discretization steps of equal size to discretize the interval [T0, T1]. The “true” option value obtained is Vsw,true = 0.1539. Note that there is empirical evidence suggesting that the value computed for Vsw,true has four significant digits correct. Hence we evaluate the relative error of the approximate formula (12), that we denote by RelErr, as follows: RelErr = (0.1539 - 0.1538)/0.1539 = 6.5 x 10-4.
Numerical experiments reveal that the Monte Carlo simulation converges to Vsw,true, up to an error on the fourth significant decimal digit, also when the Euler-Maruyama scheme is applied with only 80 time discretization steps of equal size in the interval [T0, T1], that is with the previous choices the time discretization error of the stochastic differential equations (6) is negligible compared to the other errors present in the computation.
Figure 1 shows the (relative) numerical error
obtained employing both the approximation formula (12) and the
The computer time needed to perform 700.000
Monte Carlo simulations with 80 time
discretization steps in the simulation of (6), (8) is approximatively 583.2
seconds using the Pentium IV processor.
We conclude that the computer time needed to calculate the swaption
price with the accuracy discussed previously using the numerical method
described here and proposed in  is about 583.2/0.17 ≈ 3430 times
smaller than the computer time necessary to the
Figure 1, Test Case A, error of the
analytical formula (12) compared to the
error of the
The numerical simulation relative to Test Case
B and Test Case C is performed similarly to the simulation relative to Test
Case A. The results obtained are summarized in Table 1, where we show, for Test
Cases A, B, and C, the swaption price computed using the analytical
approximation (12), Vsw,an, the relative error RelErr
of the approximation (12) compared to a true value obtained with an
accurate Monte Carlo simulation, the computer time necessary to evaluate
formula (12), CPUTimean,
the computer time necessary to evaluate the swaption value using the Monte
Carlo simulation with approximately the same accuracy obtained using the
analytical method, CPUTimeMC.
We note that the relative error obtained using the approximation (12) is in
Test Case B, equal to 3.5 x 10-3, and, in Test Case C, equal to 1.1 x 10-2. Moreover, in Test
Case B, the analytical approximation
(12) is 62.5/0.17≈ 368 times faster than the
Test Case A
6.5 x 10-4
Test Case B
3.5 x 10-3
Test Case C
1.1 x 10-2
We have applied the numerical method described
above to price European swaptions considering tenor structures, volatility
functions, strike prices and initial data different from those used in the
simulation presented here. We have also employed the numerical method to the
problem of pricing interest rate spread options. The results obtained are
contained in  , and are very similar to those presented here, in fact in the
numerical experience reported in  the numerical method proposed appears to
be accurate (i.e.: relative error of order 10-2, 10-3, or
even 10-4) and approximately between 278 and 63000 times faster than
the Monte Carlo simulation when results of approximately the same accuracy are
produced. The results obtained demonstrate the possibility of pricing interest
rates derivatives under the LMM with one driving factor using numerical methods
based on an appropriate combination of mathematical analysis and elementary numerical
methods rather than on relying exclusively on the
 L. V. Ballestra, G. Pacelli, F.Zirilli, A Numerical Method to Price Interest Rate Derivatives based on the One Factor LIBOR Market Model of Interest Rates, Nonlinear Analysis: Hybrid Systems 2 (2008) 568-589.
LIBOR Market Model
The codes made available are: Fmain.f, eigeis.f.
LIBORS_WAPTION.zip contains the following files:
Fmain.f: a Fortran main code that computes Vsw(T0).
eigeis.f: a Fortran code that computes the eigenvalues and eigenvectors of a real symmetric matrix. This code is part of the EISPACK library routines, and has been taken from http://www.netlib.org/eispack .