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

1.

2.

3.

4.

5.

7.

1. Abstract

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.

2. Mathematical formulation

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: (1)

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: (2)

with initial conditions: (3)

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 (4)

where e=Min(i,q+1).

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: (5)

Using the Ito’s lemma and setting p=1, the set of stochastic differential equations (2)-(3) can be rewritten as: (6) (7)

where (8)

As announced previously we derive a closed-form approximation of the probability density functions associated to the model (6)-(8)  .

3. The numerical method

Let h be a positive integer such that hMin(n,q+1),  Zh(t) = (Zh(t), Zh+1(t),…, Zn(t)), zh(t) = ( zh(t), zh+1(t),…, zn(t)),

z´h(t) = (h(t),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: (9)

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: (10)

where ,      (11)

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: ,      (12)

where (13)

and (14)

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.

4. Numerical results

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: (15)

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 Monte Carlo simulation. The equations of the dynamics of the LIBOR interest rates (equations (6)-(7) with h=1, q=0) are discretized in time using the Euler-Maruyama scheme  .

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 Monte Carlo simulation with 80 Euler-Maruyama time discretization steps. We note that the numerical method proposed in this paper is more accurate than the Monte Carlo simulation if the size of the Monte Carlo sample is smaller than (about) 700.000.

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 Monte Carlo method to obtain a result of similar accuracy. Figure 1, Test Case A, error of the analytical  formula (12) compared to the error of the Monte Carlo simulation

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 Monte Carlo simulation, and, in Test Case C, the analytical approximation (12) is 65.8/0.17≈ 387 times faster than the Monte Carlo simulation.

Table 1

 Vsw,an RelErr CPUTimean CPUTimeMC Test Case A 0.1538 6.5 x 10-4 0.17 s 583.2 s Test Case B 0.04816 3.5 x 10-3 0.17 s 62.5 s Test Case C 0.005668 1.1 x 10-2 0.17 s 65.8 s

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 Monte Carlo simulation.

5. References

     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.

     F. Black, The Pricing of Commodity Contracts, Journal of Financial Economics 3 (1976) 167-179.

     A. Brace, M. Gatarek, M. Musiela, The Market Model of Interest Rate Dynamics, Mathematical Finance 7 (1997) 127-155.

     D. Brigo, F. Mercurio, Interest Rate Models: Theory and Practice, Springer-Verlag, Berlin (2001).

     D. Brigo, F. Mercurio, M. Morini, The LIBOR Model Dynamics: Approximations, Calibration and Diagnostics, European Journal of Operational Research 163 (2005) 30-51.

     F. Jamshidian, LIBOR and Swap Market Models and Measures, Finance and Stochastics 1 (1997) 293-330.

     P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag, New York (1992).

     K. Miltersen, K. Sandmann, D. Sondermann, Closed Form Solutions for Term Structure Derivatives with Lognormal Interest Rates, Journal of Finance 52 (1997) 409-430.

6.  LIBOR_SWAPTION: A Fortran code to compute the price of an European swaption using the one factor

LIBOR Market Model

The codes made available are: Fmain.f, eigeis.f.

The archive 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 . 