AN EFFICIENT LATTICE ALGORITHM FOR THE LIBOR MARKET MODEL Tim Xiao Journal of Derivatives, 19 (1) 25-40, Fall 2011 ABSTRACT The LIBOR Market Model has become one of the most popular models for pricing interest rate products. It is commonly believed that Monte-Carlo simulation is the only viable method available for the LIBOR Market Model. In this article, however, we propose a lattice approach to price interest rate products within the LIBOR Market Model by introducing a shifted forward measure and several novel fast drift approximation methods. This model should achieve the best performance without losing much accuracy. Moreover, the calibration is almost automatic and it is simple and easy to implement. Adding this model to the valuation toolkit is actually quite useful; especially for risk management or in the case there is a need for a quick turnaround. Key Words: LIBOR Market Model, lattice model, tree model, shifted forward measure, drift approximation, risk management, calibration, callable exotics, callable bond, callable capped floater swap, callable inverse floater swap, callable range accrual swap. 1 The LIBOR Market Model (LMM) is an interest rate model based on evolving LIBOR market forward rates under a risk-neutral forward probability measure. In contrast to models that evolve the instantaneous short rates (e.g., Hull-White, BlackKarasinski models) or instantaneous forward rates (e.g., Heath-Jarrow-Morton (HJM) model), which are not directly observable in the market, the objects modeled using the LMM are market observable quantities. The explicit modeling of market forward rates allows for a natural formula for interest rate option volatility that is consistent with the market practice of using the formula of Black for caps. It is generally considered to have more desirable theoretical calibration properties than short rate or instantaneous forward rate models. In general, it is believed that Monte Carlo simulation is the only viable numerical method available for the LMM (see Piterbarg [2003]). The Monte Carlo simulation is computationally expensive, slowly converging, and notoriously difficult to use for calculating sensitivities and hedges. Another notable weakness is its inability to determine how far the solution is from optimality in any given problem. In this paper, we propose a lattice approach within the LMM. The model has similar accuracy to the current pricing models in the market, but is much faster. Some other merits of the model are that calibration is almost automatic and the approach is less complex and easier to implement than other current approaches. We introduce a shifted forward measure that uses a variable substitution to shift the center of a forward rate distribution to zero. This ensures that the distribution is symmetric and can be represented by a relatively small number of discrete points. The shift transformation is the key to achieve high accuracy in relatively few discrete finite nodes. In addition, we present several fast and novel drift approximation approaches. Other concepts used in the model are probability distribution structure exploitation, numerical integration and the long jump technique (we only position nodes at times when decisions need to be made). 2 This model is actually quite useful for risk management because normally fullrevaluations of an entire portfolio under hundreds of thousands of different future scenarios are required for a short time window (see FinPricing (2011)). Without an efficient algorithm, one cannot properly capture and manage the risk exposed by the portfolio. The rest of this paper is organized as follows: The LMM is discussed in Section I. In Section II, the lattice model is elaborated. The calibration is presented in Section III. The numerical implementation is detailed in Section IV, which will enhance the reader's understanding of the model and its practical implementation. The conclusions are provided in Section V. I. LIBOR MARKET MODEL Let ( , F ,  0ttF ,P ) be a filtered probability space satisfying the usual conditions, where  denotes a sample space, F denotes a  -algebra, P denotes a probability measure, and   0ttF denotes a filtration. Consider an increasing maturity structure NTTT = ...0 10 from which expiry-maturity pairs of dates ( 1−kT , kT ) for a family of spanning forward rates are taken. For any time 1− kTt , we define a rightcontinuous mapping function )(tn by )(1)( tntn TtT − . The simply compounded forward rate reset at t for forward period ( 1−kT , kT ) is defined by         −== −− 1 ),( ),(1 ),;(:)( 11 k k k kkk TtP TtP TTtFtF  (1) where ),( TtP denotes the time t price of a zero-coupon bond maturing at time T and ),(: 1 kkk TT −=  is the accrual factor or day count fraction for period ( 1−kT , kT ). Inverting this relationship (1), we can express a zero coupon bond price in terms of forward rates as: 3  = + = k tnj jj tnk tF TtPTtP )()( )(1 1 ),(),(  (2) LIBOR Market Model Dynamics Consider a zero coupon bond numeraire ),( iTP • whose maturity coincides with the maturity of the forward rate. The measure iQ associated with ),( iTP • is called iT forward measure. Terminal measure NQ is a forward measure where the maturity of the bond numeraire ),( NTP • matches the terminal date NT . For brevity, we discuss the one-factor LMM only. The one-factor LMM (Brace et al. [1997]) under forward measure iQ can be expressed as If iTtki  , , tkk k ij jj jjj kkk dXtFtdt tF tFt tFttdF )()( )(1 )()( )()()( 1     + + =  += (3a) If 1, −= kTtki , tkkk dXtFttdF )()()( = (3b) If 1, − kTtki , tkk i kj jj jjj kkk dXtFtdt tF tFt tFttdF )()( )(1 )()( )()()( 1     + + −=  += (3c) where tX is a Brownian motion. There is no requirement for what kind of instantaneous volatility structure should be chosen during the life of the caplet. All that is required is (see Hull-White [2000]): ( )  − − − == 1 0 2 1 2 1 2 )( 1 ),(:)( kT k k kkk duu T T   (4) where k  denotes the market Black caplet volatility and  denotes the strike. Given this equation, it is obviously not possible to uniquely pin down the instantaneous volatility function. In fact, this specification allows an infinite number of choices. People often assume that a forward rate has a piecewise constant instantaneous volatility. 4 Here we choose the forward rate )(tFk has constant instantaneous volatility regardless of t (see Brigo-Mercurio [2006]). Shifted Forward Measure The )(tFk is a Martingale or driftless under its own measure kQ . The solution to equation (3b) can be expressed as       +−=  t sk t kkk dXsdssFtF 00 2 )()( 2 1 exp)0()(  (5) where ),;0()0( 1 kkk TTFF −= is the current (spot) forward rate. Under the volatility assumption described above, equation (5) can be further expressed as         +−= tk k kk XtFtF   2 exp)0()( 2 (6) Alternatively, we can reach the same Martingale conclusion by directly deriving the expectation of the forward rate (6); that is ( ) )0( 2 exp 2 1 )0( 2 )( exp 2 1 )0( 2 exp 2 exp 2 1 )0()( 22 22 0 kt t kt kt k t t tk k kk FdY t Y t FdX t tX t F dX t X Xt t FtFE =         −=         − −=         −         +−=    −  −  −       (7) where tX , tY are both Brownian motions with a normal distribution (0, t) at time t, )|(:)( tt EE F•=• is the expectation conditional on the tF , and the variable substitution used for derivation is ktt tXY −= (8) This variable substitution that ensures that the distribution is centered on zero and symmetry is the key to achieve high accuracy when we express the LMM in discrete finite form and use numerical integration to calculate the expectation. As a matter of 5 fact, without this linear transformation, a lattice method in the LMM either does not exist or introduces too much error for longer maturities. After applying this variable substitution (8), equation (6) can be expressed as         +=         +−= tk k ktk k kk YtFXtFtF     2 exp)0( 2 exp)0()( 22 (9) Since the LMM models the complete forward curve directly, it is essential to bring everything under a common measure. The terminal measure is a good choice for this purpose, although this is by no means the only choice. The forward rate dynamic under terminal measure NQ is given by tkk N kj jj jjj kkk dXtFdt tF tF tFtdF )( )(1 )( )()( 1     + + −=  += (10) The solution to equation (10) can be expressed as         +−=         +−=  tk k kk t sk t k kkk XttFdXdstFtF      2 )(exp)0( 2 )(exp)0()( 2 00 2 (11a) where the drift is given by    +=+= + −=−= t N kj jk jj jj t N kj jkjk ds sF sF dsst 0 10 1 )(1 )( )()(     (11b) where  )(1/)()( sFsFs jjjjj  += is the drift term. Applying (8) to (11a), we have the forward rate dynamic under the shifted terminal measure as         ++= tk k kkk YttFtF    2 )(exp)0()( 2 (12) Drift Approximation Under terminal measure, the drifts of forward rate dynamics are statedependent, which gives rise to sufficiently complicated non-lognormal distributions. 6 This means that an explicit analytic solution to the forward rate stochastic differential equations cannot be obtained. Therefore, most work on the topic has focused on ways to approximate the drift, which is the fundamental trickiness in implementing the Market Model. Our model works backwards recursively from forward rate N down to forward rate k. The N-th forward rate )(tFN without drift can be determined exactly. By the time it takes to calculate the k-th forward rate )(tFk , all forward rates from )(1 tFk+ to )(tFN at time t are already known. Therefore, the drift calculation (11b) is to estimate the integrals containing forward rate dynamics )(sF j , for j=k+1,...,N, with known beginning and end points given by )0(jF and )(tF j . For completeness, we list all possible solutions below. Frozen Drift (FD). Replace the random forward rates in the drift by their deterministic initial values, i.e.,   +=+= + − + −= N kj jk jj jj t N kj jk jj jj k t F F ds sF sF t 10 1 )0(1 )0( )(1 )( )(        (13) Arithmetic Average of the Forward Rates (AAFR). Apply the midpoint rule (rectangle rule) to the random forward rates in the drift, i.e., ( ) ( ) += ++ + − N kj jk jjj jjj k t tFF tFF t 1 2 1 2 1 )()0(1 )()0( )(     (14) Arithmetic Average of the Drift Terms (AADT). Apply the midpoint rule to the random drift terms, i.e.,  +=         + + + − N kj kj jj jj jj jj k t tF tF F F t 1 )(1 )( )0(1 )0( 2 1 )(       (15) Geometric Average of the Forward Rates (GAFR). Replace the random forward rates in the drift by their geometric averages, i.e., 7  += +  − N kj kj jjj jjj k t tFF tFF t 1 )()0(1 )()0( )(     (16) Geometric Average of the Drift Terms (GADT). Replace the random drift terms by their geometric averages, i.e.,  += +  + − N kj kj jj jj jj jj k t tF tF F F t 1 )(1 )( )0(1 )0( )(       (17) Conditional Expectation of the Forward Rate (CEFR). In addition to the two endpoints, we can further enhance our estimate based on the dynamics of the forward rates. The forward rate )(sF j follows the dynamic (9) (The drift term is ignored). We can derive the expectation of the forward rate conditional on the two endpoints and replace the random forward rate in the drift by the conditional expectation of the forward rate. Proposition 1. Assume the forward rate )(sF j follows the dynamic (9), with the two known endpoints given by )0(jF and )(tF j . Based on the conditional expectation of the forward rate )(sF j , the drift of )(tFk can be expressed as  += + − N kj t kj tFFjj tFFjj k ds sFE sFE t jj jj 1 0 )(),0(0 )(),0(0 ]|)([1 ]|)([ )(     (18a) where the conditional expectation of the forward rate is given by         −         = t sts F tF FsFE j t s j j jtFFj jj 2 )( exp )0( )( )0(]|)([ 2 )(),0(0  (18b) Proof. See Appendix A. Conditional Expectation of the Drift Term (CEDT). Similarly, we can calculate the conditional expectation of the drift term and replace the random drift term by the conditional expectation. 8 Proposition 2. Assume the forward rate )(sF j follows the dynamic (9), with the two known endpoints given by )0(jF and )(tF j . Based on the conditional expectation of the drift term j , the drift of )(tFk can be expressed as  +=           + − N kj kj t tFF jj jj k ds sF sF Et jj 1 0 )(),0( 0 )(1 )( )(     (19a) where the conditional expectation of the drift term is given by ( ) )( )(/)(1 1 )(1 )( |)( 2 )(),0( 0)(),0(0 s ss sF sF EsE Cj CjCj tFF jj jj tFFj jj jj      + −=           + = (19b)         −         += t sts F tF Fs j t s j j jjCj 2 )( exp )0( )( )0(1)( 2  (19c)         −         −         −         = t sts t sts F tF Fs jj t s j j jjCj )( exp1 )( exp )0( )( )0()( 22 2 22   (19d) Proof. See Appendix A. The accuracy and performance of these drift approximation methods are discussed in section IV. II. THE LATTICE PROCEDURE IN THE LMM The "lattice" is the generic term for any graph we build for the pricing of financial products. Each lattice is a layered graph that attempts to transform a continuous-time and continuous-space underlying process into a discrete-time and discrete-space process, where the nodes at each level represent the possible values of the underlying process in that period. There are two primary types of lattices for pricing financial products: tree lattices and grid lattices (or rectangular lattices or Markov chain lattices). The tree lattices, e.g., traditional binomial tree, assume that the underlying process has two 9 possible outcomes at each stage. In contrast with the binomial tree lattice, the grid lattices (see Amin [1993], Gandhi-Hunt [1997], Martzoukos-Trigeorgis [2002], Hagan [2005], and Das [2011]) shown in Exhibit 1, which permit the underlying process to change by multiple states, are built in a rectangular finite difference grid (not to be confused with finite difference numerical methods for solving partial differential equations). The grid lattices are more realistic and convenient for the implementation of a Markov chain solution. This article presents a grid lattice model for the LMM. To illustrate the lattice algorithm, we use a callable exotic as an example. Callable exotics are a class of interest rate derivatives that have Bermudan style provisions that allow for early exercise into various underlying interest rate products. In general, a callable exotic can be decomposed into an underlying instrument and an embedded Bermudan option. We will simplify some of the definitions of the universe of instruments we will be dealing with for brevity. Assume the payoff of a generic underlying instrument is a stream of payments  iiiii CTFZ −= − )( 1 for i=1,...,N, where iC is the structured coupon. The callable exotic is a Bermudan style option to enter the underlying instrument on any of a sequence of notification dates ex M exex ttt ,...,, 21 . For any notification date ex jtt = , we define a right-continuous mapping function )(tn by )(1)( tntn TtT − . If the option is exercised at t, the reduced price of the underlying instrument, from the structured coupon payer's perspective, is given by ( )  = − =         − =        == N tni Ni iiii t N tni Ni i t N TTP CTF E TTP Z E TtP tI tI )( 1 )( ),( )( ),(),( )( :)( ~  (20) where the ratio )( ~ tI is usually called the reduced value of the underlying instrument or the reduced exercise value or the reduced intrinsic value. Lattice approaches are ideal for pricing early exercise products, given their "backward-in-time" nature. Bermudan pricing is usually done by building a lattice to 10 carry out a dynamic programming calculation via backward induction and is standard. The lattice model described below also uses backward induction but exploits the Gaussian structure to gain extra efficiencies. First we need to create the lattice. The random process we are going to model in the lattice is the LMM (12). Unlike traditional trees, we only position nodes at the determination dates (the payment and exercise dates). At each determination date, the continuous-time stochastic equation (12) shall be discretized into a discrete-time scheme. Such discretized schemes basically convert the Brownian motion into discrete variables. There is no restriction on discretization schemes. At any determination date t, for instance, we discretize the Brownian motion to be equally spaced as a grid of nodes tiy , , for i = 1,..., tS . The number of nodes tS and the space between nodes titit yy ,1, −−= at each determination date can vary depending on the length of time and the accuracy requirement. The nodes should cover a certain number of standard deviations of the Gaussian distribution to guarantee a certain level of accuracy. We have the discrete form of the forward rate as         ++= tik k tikktik ytytFytF , 2 ,, 2 ),(exp)0();(    (21) The zero-coupon bond (2) can be expressed in discrete form as  = + = k tnj tijj titntik ytF yTtPyTtP )( , ,)(, );(1 1 );,();,(  (22) We now have expressions for the forward rate (21) and discount bond (22), conditional on being in the state tiy , at time t, and from these we can perform valuation for the underlying instrument. At the maturity date, the value of the underlying instrument is equal to the payoff, i.e., )(),( ,, NN TiNTiN yZyTI = (23) 11 The underlying state process tX in the LMM (11) is a Brownian motion. The transition probability density from state ( tix , , t ) to state ( Tjx , , T ) is given by         − − − − = )(2 )( exp )(2 1 ),;,( 2 ,, ,, tT xx tT Txtxp tiTj Tjti  (24) Applying the variable substitution (8), equation (24) can be expressed as         − −+− − − = )(2 )( exp )(2 1 ),;,( 2 ,, ,, tT tTyy tT Tytyp tTtiTj Tjti   (25) Equation (20) can be further expressed as a conditional value on any state ( tiy , , t ) as: j jj j j T N tnj j tjTtiT TNj Tj jtiN ti dy tT tTyy yTTP yZ tTyTtP ytI  =         − −+− − − = )( 2 , , , )(2 )( exp );,( )( )(2 1 );,( );(   (26) This is a convolution integral. Some fast integration algorithms, e.g., Cubic Spline Integration, Fast Fourier Transform (FFT), etc., can be used for optimization. We use the Trapezoidal Rule Integration in this paper for ease of illustration. Incomplete information handling. Convolution is widely used in Electrical Engineering, particularly in signal processing. The important part is that the far left and far right parts of the output are based on incomplete information. Any models that try to compute the transition values using integration will be inaccurate if this problem is not solved, especially for longer maturities and multiple exercise dates. Our solution is to extend the input nodes by padding the far end values on each side and only take the original range of the output nodes. Next, we determine the option values in each final notification node. On the last exercise date, if we have not already exercised, the reduced option value in any state Miy , is given by 12         = 0, );,( );( max );,( ),( , , , , MiN ex M Mi ex M MiN ex M Mi ex M yTtP ytI yTtP ytV (27) Then, we conduct the backward induction process that is performed by iteratively rolling back a series of long jumps from the final exercise date ex Mt across notification dates and exercise opportunities until we reach the valuation date. Assume that in the previous rollback step ex jt , we calculated the reduced option value: );,(/),( ,, jiN ex jji ex j yTtPytV . Now, we go to ex jt 1− . The reduced option value at ex jt 1− is         = −− −− −− −− −− −− );,( ),( , );,( ),( max );,( ),( 1,1 1,1 1,1 1,1 1,1 1,1 jiN ex j ji ex j c jiN ex j ji ex j jiN ex j ji ex j yTtP ytV yTtP ytI yTtP ytV (28a) where the reduced continuation value is given by jex j ex j ex jj ex jjjij jN ex j j ex j ex j ex jjiN ex j ji ex j c dy tt ttyy yTtP ytV ttyTtP ytV          − −+− − − = − −−− −−− −− )(2 )( exp );,( ),( )(2 1 );,( ),( 1 2 111, 11,1 1,1   (28b) We repeat the rollback procedure and eventually work our way through the first exercise date. Then the present value of the Bermudan option is found by a final integration given by ( ) 1 1 2 111 11 11 1 2 exp );,( ),( 2 1 ),0()0( dy t ty yTtP ytV t TPpv ex ex N ex ex ex NBermudan          + −=   (29) The present value or the price of the callable exotic from the coupon payer's perspective is: )0()0()0( _ instrumentunderlyBermudanpayer pvpvpv −= (30) This framework can be used to price any interest rate products in the LMM setting and can be easily extended to the Swap Market Model (SMM). III. Calibration 13 First, if we choose the LMM as the central model, we need to price interest rate derivatives that depend on either or both of cap and swaption markets. Second, we will undoubtedly use various swaptions to hedge a callable exotic. It is a reasonable expectation that the calibrated model we intend to use to price our exotic, will at least correctly price the market instruments that we intend to hedge with. Therefore, in an exotic derivative pricing situation, recovery of both cap and swaption markets might be desired. The calibration of the LMM to caplet prices is quite straightforward. However, it is very difficult, if not impossible, to perfectly recover both cap and swaption markets. Fortunately for the LMM, there also exist extremely accurate approximate formulas for swaptions implied volatility, e.g., Rebonato's formula. We introduced a parameter  and set ii   = where i  denotes the market Black caplet volatility. One can choose different  for different i . For simplicity we describe one  situation here. By choosing 1= , we have perfectly calibrated the LMM to the caplet prices in the market. However, our goal is to select a  to minimize the sum of the squared differences of the volatilities derived from the market and the volatilities implied by our model for both caps and swaptions combined. In the optimization, we use Rebonato's formula for an efficient expression of the model swaption volatilities, given by ( ) ( )2Re , 2 1, 2 , 2 1, 02 , 2 , )0( )0()0()0()0( )()( )0( )0()0()0()0(1 bonato ji jjjiji ji T ji ijjijiLMM S FFww dttt S FFww T               == =    += +=  (31a) where ij =1 under one-factor LMM. The swap rate )0(,S is given by  +==   1, )0()0()0( i ii FwS (31b) 14 ( ) ( )   += += − += − + + =        1 1 1 1 1 )0(1 )0(1 )0( k k j jjk j jji i F F w (31c) Assume the calibration containing  caplets and G swaptions. The error minimization is given by ( ) ( ) = ++= −+− G j swn Nj bonato Nj M i ii 1 2 , Re ,1 2 min    (32) where swn Nj ,+  denotes the market Black swaption volatility. The optimization can be found at a stationary point where the first derivative is zero; that is,   = += = += + + = G j bonato Nj M i i G j swn Nj M i i 1 Re ,1 1 ,1        (33) In terms of forward volatilities, we use the time-homogeneity assumption of the volatility structure, where a forward volatility for an option is the same or close to the spot volatility of the option with the same time to expiry. The time-homogeneous volatility structure can avoid non-stationary behavior. In the LMM, forward swap rates are generally not lognormal. Such deviation from the lognormal paradigm however turns out to be extremely small. Rebonato [1999] shows that the pricing errors of swaptions caused by the lognormal approximation are well within the market bid/ask spread. For most short maturity interest rate products, we can use the lattice model without calibration (33). However, for longer maturity or deeply in the money (ITM) or out of the money (OTM) exotics we may need to use the calibration and even some specific skew/smile adjustment techniques to achieve high accuracy. IV. NUMERICAL IMPLEMENTATION In this section, we will elaborate on more details of the implementation. We will start with a simple callable bond for the purpose of an easy illustration and then move 15 on to some typical callable exotics, e.g., callable capped floater swap and callable range accrual swap. The reader should be able to implement and replicate the model after reading this section. Callable Bond A callable bond is a bond with an option that allows the issuer to retain the privilege of redeeming the bond at some points before the bond reaches the maturity date. For ease of illustration, we choose a very simple callable bond with a one-year maturity, a quarterly payment frequency, a $100 principal amount (A), and a 4% annual coupon rate (the quarterly coupon 1=C ). The call dates are 6 months, 9 months, and 12 months. The call price (H) is 100% of the principal. The bond spread (  ) is 0.002. Let the valuation date be 0. A detailed description of the callable bond and current (spot) market data is shown in Exhibit 2. For a short-term maturity callable bond, our lattice model can reach high accuracy even without calibration (33) and incomplete information handling. Therefore, we set 1= and ii   = . The valuation procedure for a callable bond consists of 4 steps: Step 1: Create the lattice. Based on the long jump technique, we position nodes only at the determination (payment/exercise) dates. The number of nodes and the space between nodes at each determination date may vary depending on the length of time and the accuracy requirement. To simplify the illustration, we choose the same settings across the lattice, with a grid space (space between nodes) 2/1= , and a number of nodes S=7. It covers 3)1( =−S standard deviations for a standard normal distribution. The nodes are equally spaced and symmetric, as shown in Exhibit 3. 16 Step 2: Find the option value at each final node. At the final maturity date 4T , the payoff of the callable bond in any state iy is given by ( )CAHyTVV ii +== ,min),(: 44, (34) where A denotes the principal amount, C denotes the bond coupon, and H denotes the call price. The option values at the maturity are equal to the payoffs as shown in Exhibit 3. Step 3: Find the option value at earlier nodes. Let us go to the penultimate notification date 3T . The option value in any state iy is given by ( )CVHyTVV ciii +== 3,33, ,min),(: (35) Equation (35) can be further expressed in the form of reduced value as         + == );,( , );,( min );,( : ~ 43 3, 4343 3, 3, i c i ii i i yTTP CV yTTP H yTTP V V (36a) where );,(/ 433, i C i yTTPV denotes the reduced continuation value in state iy at 3T given by ( ) ( ) ( ) ( ) ( ) ( ) ( )             − −+− −+             − −+− − − −− =       − −+− − − −− = − − =  34 2 33441 14 7 2 34 2 3344 4 34 34 34 2 3344 4 34 34 43 3, 2 )( exp),( 2 )( exp),( 22 )(exp 2 )( exp),( 2 )(exp );,( TT TTyy yTV TT TTyy yTV TT TT dY TT TTyY YTV TT TT yTTP V ij j j ij j i i c i        (36b) where  denotes the bond spread. Similarly we can compute the reduced callable bond values at 2T . All intermediate reduced values are shown in Exhibit 3. Step 4: Compute the final integration. The final integral at valuation date 0 is calculated as 17 ( ) ( ) 399.80 2 )( exp );,( ),( 2 )( exp );,( ),( 22 exp ),0( 2 )( exp );,( ),( 2 exp ),0()0( 2 2 221 142 12 7 2 2 2 22 42 2 2 2 4 2 2 22 42 2 2 2 4 =             + −+             + − − =       + − − = − − − =  T Ty yTTP yTV T Ty yTTP yTV T T TP dY T TY YTTP YTV T T TPV j j j j j j j        (37) Moreover, we need to add the present value of the coupon at 1T into the final price. The final callable bond value is given by 398.81),0()exp()0()0( 11 =−+= CTPTVV  (38) The pseudo-code is supplied in Appendix B for the implementation program. The convergence results shown in Exhibit 4 indicate what occurs for a given grid space  when we increase the number of nodes S. The speed of convergence is very fast, ensuring that a small number of grids are sufficient. All calculations are converged to 100.7518. One sanity check is that the callable bond price should be close to the straight bond price if the call prices become very high. Both of them are computed as 103.3536. Callable capped floater swap A callable capped floater swap has two legs: a regular floating leg and a structured coupon leg. The structured coupon rate of the j-th period ( jj TT ,1− ) is given by }],),(max{min[ 1 F j C jjjjjjjj KKTFAC −+=  (39) where jA is the notional amount, C j K is the rate cap, F j K is the rate floor, j is the spread and j is the scale factor. For j > 0, it is called a callable capped floater swap. For j < 0, it is called a callable inverse floater swap. 18 We choose a real middle life trade with more than 10 years remaining in its lifetime. The floating leg has a quarterly payment frequency with step-down notionals and step-up spreads. The structured coupon leg has a semi-annually payment frequency with varying notionals, spreads, scales, rate caps, and rate floors. The call schedule is semi-annual. Callable range accrual swap A callable range accrual swap has two legs: a regular floating leg and a structured coupon leg. The structured coupon rate of the j-th period ( jj TT ,1− ) is given by  += −= j ji T Tt j ijj j M RIA C 11  (40a) where    + = otherwise KtttFKif I jiiij i 0 ),;(1 maxmin  (40b) where R is the fixed rate, minjK and maxjK are the accrual range of the j-th period, ),;( +iii tttF is the LIBOR rate,  is the range accrual index term, jM is the total number of the business days in the j-th period. We choose a real 10 years maturity trade. The floating leg has a quarterly payment frequency and the structured coupon leg has a semi-annually payment frequency with varying accrual ranges. It starts with the first call opportunity being in 3 years from inception, and then every year until the last possibility being 9 years from inception. The range accrual index term is 6 months. The lattice implementation procedure for a callable capped floater swap or a callable range accrual swap is quite similar to the one for a callable bond except the valuation for the underlying instrument. 19 The convergence diagrams of pricing calculations are shown in Exhibits 5 and 6. Each curve in the diagrams represents the convergence behavior for a given grid space as nodes are increased. All of the lattice results are well converged. If the grid space is smaller, the algorithm has better convergence accuracy but a slower convergence rate, and vice verse. We benchmarked our model under different drift approximation methods with several standard market approaches, e.g., the regression-based Monte Carlo in the full LMM and the HJM trinomial tree. The model comparisons for the accuracy and speed are shown in Exhibits 7 and 8. With regards to accuracy, as expected, the FD performs very badly. AAFR and GAFR do a little better but errors go in different directions. The same conclusions can be drawn for AADT and GADT. Both CEFR and CEDT are the best. In terms of CPU times, FD, AAFR, AADT, GAFR and GADT are the same. But CEFR and CEDT are slower, especially in the callable range accrual swap case. V. CONCLUSION In this paper, we proposed a lattice model in the LMM to price interest rate products. Conclusions can be drawn, supported by the previous sections. First, the model is quite stable. The fast convergence behavior requires fewer discretization nodes. Second, this model has almost equivalent accuracy to the current pricing models in the market. Third, the implementation of the model is relatively easy. The calibration is very simple and straightforward. Finally, the performance of the model is probably the best among all known approaches at the time of writing. We use the following techniques in our model: shifted forward measure, drift approximation, probability distribution structure exploitation, long jump, numerical integration, incomplete information handling, and calibration. Combining these 20 techniques, the model achieves sufficient accuracy in relatively few time steps and discrete nodes, which makes it a very efficient method. For ease of illustration, we present the lattice model based on the Trapezoidal Rule integration. A better but slightly more complicated solution is to spline the payoff functions. The cubic spline of the option payoffs can achieve higher accuracy, especially for Greeks calculations, and higher speed. Although cubic spline takes some time, the lattice will require much fewer nodes (23 ~ 28 nodes are good enough) and can perform a much faster integration. In general, the spline method can provide a speedup factor around 3 ~ 5 times. We have implemented the lattice model to price a variety of interest rate exotics. The algorithm can always achieve a fast convergence rate. The accuracy, however, is a bit trickier, depending on many factors: drift approximation approaches, numerical integration schemes, volatility selections, and calibration, etc. Some work, such as calibration, is more of an art than a science. REFERENCE Amin, K. "Jump diffusion option valuation in discrete time." Journal of Finance, Vol. 48, No. 5 (1993), pp. 1833-1863. Brace, A., D. Gatarek, and M. Musiela. "The market model of interest rate dynamics." Mathematical Finance, Vol. 7, No. 4 (1997), pp. 127-155. Brigo, D., and F. Mercurio. "Interest Rate Models – Theory and Practice with Smiles, Inflation and Credit." Second Edition, Springer Finance, 2006. 21 Das, S. "Random lattices for option pricing problems in finance." Journal of Investment Management, Vol. 9, No.2 (2011), pp. 134-152. FinPricing, Capital market solution, https://finpricing.com/lib/FxVolIntroduction.html Gandhi, S. and P. Hunt. "Numerical option pricing using conditioned diffusions," Mathematics of Derivative Securities, Cambridge University Press, Cambridge, 1997. Hagan, P. "Accrual swaps and range notes." Bloomberg Technical Report, 2005. Hull. J., and A. White. "Forward rate volatilities, swap rate volatilities and the implementation of the Libor Market Model." Journal of Fixed Income, Vol. 10, No. 2 (2000), 46-62. Martzoukos, H., and L. Trigeorgis. "Real (investment) options with multiple sources of rare events." European Journal of Operational Research, 136 (2002), 696-706. Piterbarg, V. "A Practitioner's guide to pricing and hedging callable LIBOR exotics in LIBOR Market Models." SSRN Working paper, 2003. Rebonato, R. "Calibrating the BGM model." RISK, March (1999), 74-79. APPENDIX A: Proof of Proposition 1. We rewrite (9) as         −         = 2)0( )( ln 1 )( 2t F tF tY j j j j   (A1) 22 In the general Brownian Bridge case when the Wiener process )(tY has )( 1 tY =a and )( 2 tY =b, the distribution of )(tY at time ),( 21 ttt  is normal given by         − −− = − −− += )( ))(( )(, )( ))(( )(~)( 12 21 12 1 tt tttt t tt abtt atNtY YY  (A2) In our case: 01 =t , tt =2 , a=0, b= )(tY , ),0( ts , thus (A2) can be expressed as       − == t sts stY t s sNsY YY )( )(),()(~)(  (A3) Let 2/)()( 2ssYsA jjj  += . According to the linear transformation rule, )(sA j is a normal given by         − ==         =+= t sts ss F tF t ss sssA j YjAj j jj YjAjj )( )()(, )0( )( ln 2 )()(~)( 2 2 2     (A4) Let ( ))(exp)( sAsB jj = . By definition, )(sB j is a lognormal given by ( ))(),(~)( ssLogNsB AjAjj  . According to the characterizations of the lognormal distribution, the mean and variance of )(sB j are ( )         −         =         +== t sts F tFs sBEs j t s j jAj AjjBj 2 )( exp )0( )( 2 )( exp)()( 2 0   (A5a) ( )  ( )         −                 −         − =+−= t sts F tF t sts ssss j t s j jj AjAjAjBj )( exp )0( )( 1 )( exp)()(2exp1)(exp)( 2 2 2   (A5b) We have the conditional expectation of the forward rate )(sF j as ( ) ( )         −         == t sts F tF FsBEFsFE j t s j j jjjtFFj jj 2 )( exp )0( )( )0()()0(|)( 2 0)(),0(  (A6) 23 Proof of Proposition 2. Let )()0(1)(1)( sBFsFsC jjjjjj  +=+= where )(sB j is defined above. According to the linear transformation rule, )(sC j is a lognormal given by ( ))(),(~)( svsLogNsC j   . The mean and variance of )(sC j are         −         +=+= t sts F tF FsFs j t s j j jjBjjjCj 2 )( exp )0( )( )0(1)()0(1)( 2  (A7a)         −                 −         − == t sts F tF t sts FsFs j t s j jj jjBjjjCj )( exp )0( )( 1 )( exp)0()()0()( 2 2 2 2222   (A7b) On the other hand, according to the characterizations of the lognormal distribution, the mean and variance of )(sC j are       += 2 )( )(exp)( s ssCj     (A8a) ( )  ( ))()(2exp1)(exp)( ssssCj   +−= (A8b) Solving the equation (A8a) and (A8b), we get           + = )(/)(1 )( ln)( 2 ss s s CjCj Cj     (A9a)         += )( )( 1ln)( 2 s s s Cj Cj     (A9b) We know the first negative moment of the lognormal is ( ) ( )2/)()(exp)(1 sssCE j   +−=− and have the conditional expectation of the drift term as )( )(/)(1 1 2 )( )(exp1 )( 1 1 )(1 1 1 )(1 )( 2 00 )(),0( 0 s sss s sC E sF E sF sF E Cj CjCj jjj tFF jj jj jj      + −=      +−−=         −=         + −=           +   (A10) where )(sCj , )(sCj are given by (A7a) and (A7b). 24 APPENDIX B: The following pseudo-code (C++) demonstrates how to implement the model to price a callable bond. For the purpose of an easy illustration, we choose the same settings (the number of nodes and the grid space) across the lattice and use the Trapezoidal Rule for numerical integration. // 2*numNodes = 2*mNumNodes = the number of nodes (S); gap = mGap = the grid space (Phi) double priceCallableBond (BondTrade* bd, CallableBond* cb, int numNodes, double gap) { double pv; cb->fillLattice(); // The last exercise CallSchedule& cs = bd->callSch[numCallSch-1]; if (cs.term == bd->cFlow[numCashFlow-1].endDate) // The last exercise is at maturity for (int i= -numNodes; i <= numNodes; i++) cs.reducedValue[i+numNodes] = min (cs.callPrice, bd->cFlow[numCashFlow-1].reducedPayoff[i+numNodes]); else { // The last exercise is before maturity for (int i= -numNodes; i <= numNodes; i++) { pv = 0; for (int j = bd->numCF-1; (bd->cFlow[j].endDate >= cs.term) && (j >= 0); j--) { CashFlow& cf = bd->cFlow[j]; (cf.endDate == cs.term) ? pv += cf.reducedPayoff[i+numNodes] : pv += exp(-bondSpread*(cf.endDate-cs.term)) * cb->integral(i, cs.vol, cf.vol, cf.endDate, cs.term, cf.reducedPayoff); } cs.reducedValue[i+numNodes] = min (cs.callPrice/cs.df[i+numNodes], pv); } } if (numCallSch > 1) { // The remaining exercises for (int i = numCallSch 2; i>=0; i--) { CallSchedule& cs = bd->callSch[i]; CallSchedule& preCs = bd->callSch[i+1]; for (int j = -numNodes; j <= numNodes; j++) { pv = exp(-bondSpread * (preCs.term cs.term)) * cb->integral (j, cs.vol, preCs.vol, preCs.term, cs.term, preCs.reducedValue); for (int k=bd->numCF-1; k >= 0; k--) // Count intermediate coupons if ((bd->cFlow[k].endDate < preCs.term) && (bd->cFlow[k].endDate >= cs.term)) pv += bd->cFlow[k].reducedPayoff[j+numNodes] * exp (-bondSpread*(bd->cFlow[k].endDate cs.term)); cs.reducedValue[j+numNodes] = min (cs.callPrice/cs.df[j+numNodes], pv); } } } // The final integral CallSchedule& preCs = bd->callSch[0]; pv = cb->integral (0, 0, preCs.vol, preCs.term, 0, preCs.reducedValue) *exp(-bondSpread*(preCs.term)); pv *= bd->cFlow[bd->numCF-1].endDf; // endDf: discount factor from 0 to the end date for (int k=bd->numCF-1; k >= 0; k--) // Count intermediate coupons 25 if ((bd->cFlow[k].endDate < preCs.term)) pv += bd->cFlow[k].coupon * bd->cFlow[k].endDf * exp(-bondSpread * bd->cFlow[k].endDate); return pv; } void CallableBond::fillLattice() { for (int i = mTrade->numCF-1; i>=0; i--) { CashFlow& cf = mTrade->cFlow[i]; if (cf.endDate < mTrade->callSch[0].term) break; for (int j = -mNumNodes; j <= mNumNodes; j++) fillNode(i, j, cf.startDate, mDrift); } } void CallableBond::fillNode(int cI, int nI, double vT, DriftAppx flag) { int numCF = mTrade->numCF; double avgF, expon, fwdt, drift = 0; CashFlow& fl = mTrade->cFlow[cI]; if (cI == numCF-1) { // At maturity fl.df[nI + mNumNodes] = 1.0; fl.reducedPayoff[nI + mNumNodes] = fl.notional + fl.coupon; } else if (fl.startDate <= 0) // Starting before valuation date) fl.reducedPayoff[nI + mNumNodes] = fl.coupon * fl.endDf / mTrade->cFlow[numCF-1].endDf; else { fl.df[nI + mNumNodes] = 1.0; for (int i = numCF 1; i > cI; i--) { CashFlow& cf = mTrade->cFlow[i]; expon = (cf.vol * cf.vol * vT / 2) + cf.vol * nI * mGap; fwdt = cf.fwd0 * exp(-drift + expon); switch (flag) { // The other cases are similar to either AAFR or CEFR case AAFR: // Arithemic Average Fwd Rate avgF = 0.5 * (cf.fwd0 + fwdt); drift += vT * fl.vol * cf.vol * cf.delta * avgF / (1 + cf.delta * avgF); break; case CEFR: // Conditional Expectation of Fwd Rate drift += fl.vol * cf.vol * integralFwd(cf.fwd0, fwdt, 0, vT, cf.vol, cf.delta); break; default: break; } fl.df[nI + mNumNodes] /= (1 + fwdt * cf.delta); // df: discount factor maturing at maturity } fl.reducedPayoff[nI + mNumNodes] = fl.coupon / fl.df[nI + mNumNodes]; } } // Gauss-Legendre integration for drift const double xArray[] = {0, 0.1488743389, 0.4333953941, 0.6794095682, 0.8650633666, 0.9739065285}; const double wArray[] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513491, 0.0666713443}; double CallableBond::integralFwd(double F0, double Ft, double a, double b, double vol, double delta) { double xm = 0.5 * (b + a); double xr = 0.5 * (b a); double ss = 0, dx = 0; for (int j = 1; j <= 5; j++) { dx = xr * xArray[j]; ss += wArray[j] * (expectFwd(F0, Ft, (xm + dx), b, vol, delta) + expectFwd(F0, Ft, (xm dx), b, vol, delta)); } return ss * xr; } 26 double CallableBond::expectFwd(double F0, double Ft, double s, double t, double vol, double delta) { double mean = F0 * pow ((Ft / F0), (s / t)) * exp(0.5 * vol * vol * s * (t s) / t); return delta * mean / (1 + delta * mean);- } // Trapezoidal Rule Integration double CallableBond::integral (int curPos, double curVol, double preVol, double preTerm, double curTerm, double* value){ double diffPos, tmpV, sum = 0; for (int k = -mNumNodes; k <= mNumNodes; k++) { diffPos = k*mGap curPos*mGap + preVol * preTerm curVol * curTerm; tmpV = value[k+mNumNodes] * exp (-diffPos * diffPos/(2*(preTerm curTerm))); ((k == -mNumNodes) || (k == mNumNodes)) ? sum += 0.5 * tmpV : sum += tmpV; } return sum * mGap / sqrt(2 * PI * (preTerm curTerm)); } EXHIBIT 1. The Grid/Rectangular Lattice This exhibit defines the state space for the underlying process tY over the first two discrete time periods. The starting state 0y at valuation date 0 is the single root of the lattice. At each date it the underlying process itY is discretized into a number of vertical nodes/states indexed by j. The value itj y , denotes the underlying process in state j at date it . The node 1,1 ty , for instance, can evolve to any discrete state in 2t with certain transition probabilities. For a Brownian motion, the transition probability can be easily determined by (25). EXHIBIT 2: The Callable Bond and Associated Spot Market Data 1t 1,2 t y 1,3 t y 1,4 t y 1,5 t y 2,1 t y 2,2 t y 2,3 t y 2,4 t y 2,5 t y 0y 1,1 t y 2t 27 The callable bond has a one-year maturity, a $100 principal, a quarterly payment frequency, and a 4% annual coupon rate. Delta = (end date – start date)/365 (day count: ACT/365). The discount bond ),0( iTP matures at the end date iT . The call dates are 6, 9, and 12 months. Cash flow index 1 2 3 4 Start date (days) 0 92 181 273 End date (days) 92 ( 1T ) 181 ( 2T ) 273 ( 3T ) 365 ( 4T ) Delta (years) 0.252055 0.243836 0.252055 0.252055 Payoff ($) 1 1 1 101 Call Schedule (days) 181 273 365 Discount bond ),0( iTP 0.999313 0.998557 0.997293 0.995667 Black Volatility i  0.337631 0.344218 0.350878 EXHIBIT 3: The LMM Lattice Structure of the Callable Bond. The callable bond is defined in Exhibit 2. ),( ~ : ~ , ijji yTVV = denotes the reduced value of the callable bond at any node (i, j). 1V denotes the coupon at 1T . )0(V is the value calculated by the final integration. )0(V is the final callable bond value that is equal to )0(V plus the present value of 1V . The grid space is 5.0= and the number of nodes is 7=S . This lattice has 3 steps and 7 nodes. 0 )6(2 mT )12(4 mT)9(3 mT 5.114,13,12,1 −==== yyyy 124,23,22,2 −==== yyyy )3(1 mT 07.29 ~ 2,1 =V 99.66 ~ 2,2 =V 398.81)0( =V 5.034,33,32,3 −==== yyyy 044,43,42,4 ==== yyyy 5.054,53,52,5 ==== yyyy 1004,1 =V 1004,2 =V 1004,3 =V 1004,4 =V 1004,5 =V 24.44 ~ 3,1 =V 05.78 ~ 3,2 =V 24.96 ~ 3,3 =V 17.100 ~ 3,4 =V 82.98 ~ 3,5 =V 41.86 ~ 2,3 =V 72.96 ~ 2,4 =V 55.94 ~ 2,5 =V 7 5.0 = = S  11 =V 399.80)0( =V 1004,6 =V 1004,7 =V 164,63,62,6 ==== yyyy 5.174,73,72,7 ==== yyyy 11.87 ~ 3,6 =V 72.57 ~ 3,7 =V 36.77 ~ 2,6 =V 46.45 ~ 2,7 =V 28 EXHIBIT 4: The Convergence Results for the Callable Bond. The callable bond is defined in Exhibit 2. 1= and drift approximation is AADT. Each curve represents the convergence behavior for a given grid space (phi) as nodes are added. All calculations are converged to 100.7518. . EXHIBIT 5: The Convergence Results for the Callable Capped Floater Swap The callable capped floater swap has more than 10 years remaining in its lifetime. The floating leg has a quarterly payment frequency. The structural leg has a semi-annually payment frequency. The call schedule is semi-annual. =1 and drift approximation is CEDT. Each curve represents the convergence behavior for a given grid space (phi) as nodes (N) are added. Convergence of a callable bond 70 80 90 100 110 11 15 19 23 27 31 35 39 43 47 Number of nodes N P ri c e s phi=1/2 phi=1/3 phi=1/4 29 EXHIBIT 6: The Convergence Results for the Callable Range Accrual Swap The callable range accrual swap has 10 years maturity. The floating leg has a quarterly payment frequency. The structural leg has a semi-annually payment frequency. There are 7 call opportunities.  =1 and drift approximation is CEDT. Each curve represents the convergence behavior for a given grid space (phi) as nodes are added. EXHIBIT 7: The Benchmark Results for the Callable Capped Floater Swap Convergence of a callable capped floater swap -9.0% -8.0% -7.0% -6.0% -5.0% -4.0% -3.0% -2.0% -1.0% 0.0% 1.0% 2.0% 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 Number of nodes N R e la ti v e p ri c e e rr o r phi=1/4 phi=1/6 phi=1/8 phi=1/10 Convergence of a callable range accrual swap -5.0% -4.0% -3.0% -2.0% -1.0% 0.0% 1.0% 2.0% 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 Number of nodes N R e la ti v e p ri c e e rr o rs phi=1/4 phi=1/6 phi=1/8 phi=1/10 30 This exhibit presents the results for model comparison. We benchmark the lattice model under different drift approximation methods with several standard market approaches, e.g., the regression-based Monte Carlo in the full LMM and the HJM trinomial tree, for both accuracy and speed. The trade is the same as the one in Exhibit 5. The grid space is  =1/8 and the number of nodes is S=200. PC denotes Predictor-Corrector. The column 'Dif from MC' = 1 – (current row price) / (price of MC in LMM). All computational times are denoted in seconds on a computer with a 2.33 GHz Duo Core CPU. Model  Drift Steps n Calls Nodes/Paths Price Err from MC Run time MC in LMM PC 40 20 1 million 4,546,863.3 0 290.32 HJM tri-tree - 1979 20 2n+1 4,602,136.3 1.22% 15.01 Our Model 1 FD 40 20 200 4,822,728.4 6.07% 0.32 1 AAFR 40 20 200 4,637,263.2 1.99% 0.32 1 AADT 40 20 200 4,637,718.1 2.00% 0.32 1 GAFR 40 20 200 4,698,215.6 3.33% 0.32 1 GADT 40 20 200 4,698,441.3 3.33% 0.32 1 CEFR 40 20 200 4,665,210.3 2.60% 0.38 1 CEDT 40 20 200 4,665,552.4 2.61% 0.39 0.99 FD 40 20 200 4,708,768.9 3.56% 0.32 0.99 AAFR 40 20 200 4,504,989.2 -0.92% 0.32 0.99 AADT 40 20 200 4,505,426.3 -0.91% 0.32 0.99 GAFR 40 20 200 4,609,779.5 1.38% 0.32 0.99 GADT 40 20 200 4,609,996.6 1.39% 0.32 0.99 CEFR 40 20 200 4,563,689.2 0.37% 0.38 0.99 CEDT 40 20 200 4,563,730.9 0.37% 0.39 EXHIBIT 8: The Benchmark Results for the Callable Range Accrual Swap This exhibit presents the results for model comparison. We benchmark the lattice model under different drift approximation methods with several standard market approaches, e.g., the regression-based Monte Carlo in the full LMM and the HJM trinomial tree, for both accuracy and speed. The trade is the same as the one in Exhibit 6. The grid space is  =1/8 and the number of nodes is S=200. The column 'Dif from MC' = 1 – (current row price) / (price of MC in LMM). All computational times are denoted in seconds on a computer with a 2.33 GHz Duo Core CPU. Model  Drift Steps n Calls Nodes/Paths Price Dif from MC Run time MC in LMM Euler 1801 7 1 million 585793.2 0.00% 2372.21 HJM tri-tree - 1801 7 2n+1 582167.8 -0.62% 15.62 31 Our Model 1 FD 1801 7 200 648365.4 10.68% 0.21 1 AAFR 1801 7 200 602482.2 2.85% 0.21 1 AADT 1801 7 200 602742.1 2.89% 0.21 1 GAFR 1801 7 200 616318.6 5.21% 0.21 1 GADT 1801 7 200 616425.3 5.23% 0.21 1 CEFR 1801 7 200 598253.3 2.13% 2.21 1 CEDT 1801 7 200 598372.4 2.15% 2.35 0.99 FD 1801 7 200 609373.9 4.03% 0.21 0.99 AAFR 1801 7 200 579337.2 -1.10% 0.21 0.99 AADT 1801 7 200 579386.3 -1.09% 0.21 0.99 GAFR 1801 7 200 591981.5 1.06% 0.21 0.99 GADT 1801 7 200 591917.6 1.05% 0.21 0.99 CEFR 1801 7 200 588918.9 0.53% 2.21 0.99 CEDT 1801 7 200 588935.7 0.54% 2.