A Comparison of Control Variate Methods for Pricing Interest Rate Derivatives in the LIBOR Market Model

1 Introduction

In the past decades, many exotic interest rate derivatives have been developed to meet the needs of investors, such as Caps, Floors, Swaptions, and so forth. In these derivatives, interest rate is used not only for discounting as well as for defining the payoffs of the products. Therefore, choosing a suitable interest rate model is extremely significant. The LIBOR market model (LMM), which was developed by Brace, Gatarek and Musiela^[1], has become a standard model for the pricing of interest rate derivatives in recent years. Meanwhile, various details of its implementation are still being worked out, just like how to approximate the drift of the LMM given by a stochastic differential equation (SDE), for better approximations will lead to more accurate and efficient Monte Carlo simulations. Miltersen, Sandmann and Sondermann^[2] drove a unified structure of interest rates model and provided the closed form solutions for some rate derivatives. Brace, Gatarek and Musiela^[1] introduced the Frozen drift approximation. Glasserman and Zhao^[3] studied the arbitrage-free discretizational approximations of the LMM. Hull and White^[4], Jackel and Rebonato^[5] developed the methods of approximation for the LMM by high dimensional lognormal processes. Hunter, Jackel and Joshi^[6] approximated the drift of interest rate model over the mentioned long time steps to reduce the calculation considerably. Kurbanmuradov, Sabelfeld and Schoenmakers^[7] studied the neglecting the drift term approximation, the Log-normal approximation and others. Joshi and Stacey^[8] studied the tpc, mpc, capc and cami drift term of the LMM approximation methods. Siopacha and Teichmann^[9], Papapantoleon and Siopacha^[10] developed the strong Taylor approximation for the drift term of the LMM while Papapantoleon and Skovmand^[11] studied the Picard approximation for it.

On the other hand, it is not easy to solve explicitly the SDE corresponding to the LMM in pricing the interest rate derivatives, therefore it requires the use of some numerical methods to solve them. Joshi and Stacey^[8] and some others mentioned above used Monte Carlo Simulations in the pricing of LIBOR rate derivatives. Furthermore, Pelsser, Jong, Driessen^[12] used Monte Carlo method to price Caps within the LMM. But they had no discussion about acceleration method to improve the Monte Carlo simulation efficiency and reduce the variance in the simulations, the latter plays key role in increasing the efficiency of simulations.

The method of control variate is one of the most widely used variance reduction techniques and so it is among the most effective applicable techniques for improving the efficiency of Monte Carlo simulations. It exploits information about the errors in estimates of known quantities to reduce the errors in an estimate of an unknown quantity. To illustrate the method of control variate, we consider the problem of estimating the expectation of a random variable V , the discounted payoff of a derivative, such as caplet, Caps, Swaption, etc. Let V₁, V₂, · · · , V_m be the outputs from m replications and suppose that each V_i is independent and identically distributed (i.i.d.). The estimator of the price is of the average

Suppose that on each simulation there is another output X_j along with V_j and the pairs (V_j, X_j) are i.i.d.. Meanwhile, the expectation E[X_j] = E[X] is assumed known. Then for any fixed coefficient b we can calculate

from the jth simulation. Therefore the control variate estimator of the price is given by

It is obvious that V (b) = V - b(X - E[X]) is an unbiased estimator of V and its variance is

To minimize the variance, the optimal coefficient b∗=Cov[V,X]Var[X] is selected and then

Consequently, we should choose proper control variate X which has high correlation with V to reduce the variance.

In this paper, we will use the control variate method to improve the efficiency of Monte Carlo simulation in evaluating the value of some interest rate derivatives in the LIBOR market model. Firstly, we will present three approximation methods which were widely used in practice and introduce two new ones based on direct simulations. Secondly, we will combine these approximations with control variate method to reduce the variance of simulation estimates and increase the efficiency of Monte Carlo simulations. We focus on the LIBOR market model for a discrete set of tenors given by a stochastic differential equation (SDE) in the terminal bond measure as developed in [13]. We will compare the efficiencies of simulations in the valuation of the LIBOR derivatives for different simulation methods based on the different control variates.

The rest of the paper is organized as follows. In Section 2, we recall three approximation methods of the LIBOR market model, the neglecting the drift term approximation, the frozen drift approximation and the Log-normal approximation, which had been studied in some literatures and propose two new ones based on direct simulations, the optimal linear function approximation and the optimal piecewise function approximation. In Section 3, we construct different control variate approximates, based on the approximation results in Section 2. As an example, we will demonstrate the simulation efficiency for evaluating the value of Caplets by the method of control variate acceleration. We will also analyze the simulation results with different control variates as well as various parameters in LIBOR market model. The idea in this paper may also applicable for pricing other interest rate derivatives base on the LIBOR market model.

2 Approximations of LIBOR Market Model

For a given tenor structure 0 = T₀ < T₁ < ··· < T_n < T_n+1, we consider a LIBOR market model for the forward LIBOR process L_i(t) in the terminal bond numeraire Qⁿ⁺¹ (see, e.g., [21])

where

L_i(t)is defined in the interval [T₀, T_i], δ_i = T_i+1 – T_i (i = 1, 2, · · · , n) and σ_i,j(t) (j = 1, 2,· · · ,d) are given deterministic functions also defined in [T₀, T_i], σ_i(t) = (σ_i,₁(t),σ_i,₂(t), · · · ,σ_i,_d(t))^T. W₁(t), W₂(t), · · · , Wd(t)are Brown motions independent with each other, W(ⁿ⁺¹)(t) = (W₁(t), W₂(t),··· ,W_d(t))^T (t ∈ [T₀,T_n]) is a standard d-dimensional Wiener process under Qⁿ⁺¹. When d = 1, (1) becomes a one-factor model with σ_i(t) = σ_i,1(t), W(ⁿ⁺¹)(t) = W₁(t).

Next, we present several approximation methods. Among them, the neglecting the drift term approximation is studied by lots of authors, the frozen drift approximation was first introduced by Brace, Gatarek and Musiela^[1], the Log-normal approximation was developed by Kurbanmuradov, Sabelfeld and Schoenmakers^[7].

Method 1: Neglecting the drift term approximation

Neglecting the random drift term in (1), we have that

with μi(1)(t)=0. So the solution of (1) can be approximated by

Method 2: Frozen drift approximation

Replacing the random drift term in (1) by their deterministic initial values, (1) can be approximated by

where

Then the solution of (5) can be given by

The above approximation was first proposed by Brace et al^[1].

Method 3: Log-normal approximation

Let f be a function defined by f(x)=δix1+δix and Y_i(t) = f(L_i(t))(i = 1, 2, · · · ,n). By (1) and using Ito formula, we find that Y_i(t) satisfies the following SDE

The first term οf the Picard Iteration to the solution of above SDE is Yi(0)(t)≈Yi(0)=δiLi(0)1+δiLi(0), which can be regarded as the approximation in the method 2. The second term o the Picard Iteration to the solution of above SDE is

Replacing δlLl(t)1+δlLl(t)in expression of μ_i (t)in 2 by Y_l^(t)(t), we get

Thus the solution of (1) can be approximated by

Next, we will propose two new approximations for the random drift function μ_i(t) by μ⁽⁴⁾(t) and μ⁵(t) respectively, based on direct simulations.

Method 4: Optimal linear function approximation based on direct simulations

Divide each interval [T_i,T_i+1] into K equal parts with mesh size Δt = t_k+1 — t_k. Based on the Euler discretization for SDE (1), generate a standard normal vector Z = (Z₁, Z₂, · · · , Z_n) each Z_i (i = 1, 2, · · · , n) is a standard normal random variable and independent with each other. For the jth simulation, let

where

Let m be the simulation times and

Choosing a_i, b_i (i = 1, 2, · · · , n), such that (a1,b1,⋯,an,bn)=arg⁡minai,biJ(a1,b1,⋯,an,bn). Then

lead to

So (1) can be approximated by equation

and the corresponding solution of (1) can by approximated by

Method 5: Optimal piecewise function approximation based on direct simulations

We define piecewise functions μi(5)(t), such that μi(5)(tk)=ai,k for t ∈ [t_k-1, t_k], and let a = (a₁₁, a₂₁, a₂₂ · · · ,a_nn),

Choosing a_i,_k such that a=arg⁡minai,kJ(a), then

lead to

Thus the solution of (1) can be approximated by

3 Numerical Results for Pricing of Caplets

Due to a Cap is a collection of caplets and that each caplet may be viewed as a call option on a forward rate, so we primarily analyze simulation results for each caplet next. In order to evaluate the value of caplet, different Monte Carlo simulation methods are used based on the different control variates.

The value of a caplet_i, defined in the accrual period [T_i, T_i+1] (i = 1, 2, · · · , n – 1), can be represented in the numeraire Qⁿ⁺¹ by (see, e.g., [14])

where B_n+1(t) is the time t price of zero-coupon bond with maturing at T_n+1 (0 ≤ t ≤ T_n+1),

The algorithm for the valuation of a caplet with control variate technique is then given as follows.

i) Divide each interval [T_i, T_i+1] into K equal parts with mesh size Δt = t_k+1 – t_k. Based on the Euler discretization for SDE (1), generate a standard normal vector Z = (Z₁, Z₂, · · · , Z_d), Z_i (i = 1, 2, · · · , d) are standard normal random variables and independent with each other,

then the jth (j = 1, 2, · · · , m) replication of the LIBOR value L_i(t) is attained, and the price of caplet_i is given by

ii) Let

with deterministic drift functions μ^i(t), which will replaced by μi(1)(t),⋯,μi(5)(t)respectively. Then for any k ≤ i – 1,

and

where

iii) Set

then

It is easy to see that

where

iv) The control variate estimate of the price of the caplet_i is finally obtained from the mean of m replications

In practice, b=cov(Vcapleti,Vcapl^eti)var(Vcapl^eti) is given by

where VcapletijandVcapl^etij are calculated by (19) and (20), respectively.

Now, we present some computation results based on the control variate method illustrated above under the multi-factor model. Taking the initial LIBOR value as L_i(0) = 0.05, the strike price as K = 0.05 and σ_i,1(t) = · · · = σ_i,_d(t) = σ (d = n, i = 1,2, · · · ,n).

Figure 1

Scatter plots of values of caplet and control variate for n = 2, L_i(0) = 0.05, K = 0.05, σ = 0.2, ΔT = 0.5, △t=130, m = 10⁴

Firstly, we present the simulation results for one caplet and the tenor structure is 0 = T₀ < T₁ < T₂ < T₃, the terminal numeraire is Q³. Figure 1 shows scatter plots of simulated values of B3(0)δ1(L1(T1)−K)+B3(T2) with approximated values B^3(0)δ1(L1(s)(T1)−K)+B^3(T2)(serving as the corresponding control variates for s = 1,2, ⋯,5) in the case of n = 2, σ = 0.2, mesh size △t=130,δi = ΔT = T_i+1 – T_i = 0.5 (i = 1, 2, 3), simulation times m = 10⁴. The plot shows the strong correlation between the simulated caplet values and the approximated caplet values of different control variate and their correlations are all more than 0.9995.

Let Error=∑j=1m(Vcapletj−V¯caplet)2/m(m−1)be the simulations error with classical Monte Carlo method and Errors=∑j=1m(Vcapletj(b)−V¯caplet(b))2/m(m−1) (s = 1, 2, 3, 4, 5) be the sth control variate Monte Carlo method as described above. Then Tables 1~3 show that the control variate simulations have smaller errors than the classical Monte Carlo simulation errors. Moreover, the errors decrease when Δt → 0 and m → ∞.

Secondly, we present the simulation results for three caplets, the tenor structure is 0 = T₀ < Τ₁ < T₂ < T₃ < T₄ < T₅, and the terminal numeraire is Q⁵. Let Rs=ErrorErrors(s=1,2,3,4,5) be the variance reduction ratios. Then Table 4 shows the computation results of the variance reduction ratios for different control variates method described above and various simulation paths. Each R_i is a 3-dimension vector for there are 3 caplets.

Table 5 shows the computation results for different control variate as described above and different volatilities σ. We can see that the less σ is, the greater variance reduction ratios are.

Table 6 shows the simulation results for different control variates, and different mesh sizes Δt. We can see that the less Δt is, the greater variance reduction ratios are.

Table 7 shows the simulation results for different control variates and different tenor sizes ΔT with σ = 0.2, m = 10³, △t=130.

From Tables 4~7, we find that the control variate methods show obvious variate reduction effects. While the less volatility σ or mesh size Δt or tenor size ΔT we take, the bigger variance reduction ratios R we get. They also show that the Log-normal approximation, the linear optimal approximation exhibits a little better acceleration effect than other three.

We also compare the computation results for one-factor model with multi-factor model in Table 8, which shows that the variance reduction ratios of caplets based on one-factor model of LIBOR rate are greater than that based on multi-factor model of LIBOR.

Finally, we present the computational results in Table 9 for different numbers of caplets, which shows that the variance reduction ratios of caplets are evidently reduced when we take more caplets based on the approximation of LMM.

4 Conclusion

In this paper, we have proposed different control variates for Monte Carlo simulations, we illustrate the effects for different control variates of the underlying process for Monte Carlo simulation. The idea in this paper can also be extended to pricing other interest rate derivatives based on LIBOR market model, such as different kinds of Swaps, and so forth, by constructing control variate H(L^(T)), where H(L(T)) is the payoff of the interest rate derivative based on the LIBOR market model.