Yield Curve Modelling with Skews and Stochastic Volatility by Leif Andersen and Jesper Andreasen Bank of America Securities April 2002 Modified August 2002 Abstract This paper discusses a variety of techniques for modeling the evolution of interest rates in the presence of stochastic volatily. We cover both Libor Market models and low-dimensional Markovian HJM models. To facilitate model calibration, special emphasis is put on efficient pricing of plain-vanilla instruments. Effects of stochastic volatility on CMS structures as well as Bermudan swaptions are discussed. 1. Introduction. There is currently great interest in improving fixed income models to better capture observed volatility skews and smiles in cap and swaption markets. Approaches suggested in the literature include state-dependent diffusions (see e.g. Andersen and Andreasen, 2000) and geometric Brownian motion overlaid with a jump process (as in Glasserman and Kou, 1999). While these approaches are both reasonable, they are not without problems and it is unlikely that they tell the whole story. For instance, as discussed in Rebonato (2001) and in a number of empirical studies, inclusion of stochastic volatility into fixed income models can significantly improve the realism of the models. This paper continues this line of research by focusing on flexible, yet practical techniques for the financial engineer to construct interest rate models with stochastic volatility. We work with both multi-factor Libor Market (LM) models and low-dimensional Markovian HJM models, and demonstrate that the proposed models are reasonable from an empirical perspective. To aid calibration and fast mark-to- market of simple derivatives, we pay particular attention to the development a variety of exact and approximative techniques to compute the prices of basic European derivatives such as caps, swaptions, and CMS options. Further, the paper introduces numerical methods for complex derivatives and gives examples of the effects of stochastic volatility on Bermudan swaptions. 2. Notation. Libor Market model. Let {T}N be a discrete tenor structure and let F (t) denote the time t value of the i i=1 k discretely compounded forward rate spanning [,T T] . Also, let R (t) denote the par k k+1 s,e rate of a plain-vanilla swap exchanging fixed for floating payments at points in the tenor structure TT ,,..., T . With d ” T - T and P(,t T) being the time t discount factor to ss+1 +2 e ii +1i time T, we have FtP()(t,T)=/P-(d,)-t11T,(t£ T ; ) kk kk k +1 P(,t)T(P, t-)T (cid:229) e Rste,A,se(ti)P,(=)t(T=,), it eAT1 .(ts) s £ d- s,e i=s+1 We now postulate model arbitrage-free forward rate dynamics of the form dFtF()t(V)=(t)tjV()lt((t)d(t)d(W))t m Tغ + øß (1a) kkk k where W is an n-dimensional Brownian motion, j: + fi + is a well-behaved ¡ ¡ deterministic function satisfying j(0)=0, l is an n-dimensional deterministic volatility k function loading each scalar Brownian motion, V(t) is a scalar positive process to be specified, and m is an n-dimensional numeraire-specific drift that ensures lack of arbitrage k across bonds. An expression for m can be found in Andersen and Andreasen (2000); for k the special case where the bond P(t, T ) is used as the numeraire asset m ()t =0 and the k+1 k kth forward rate is a martingale. Notice that (1a) is non-standard as we have allowed the local variance ||(l)||t 2 to be multiplied by a scalar factor V(t). Setting V()t ” 1 (or some k 2 deterministic function) recovers the extended Libor market model of Andersen and Andreasen (2000); letting V be random introduces the desired stochastic movements of rate volatilities. We note that these movements of the volatility surface are "parallel"; additional factors could, in principle, be added to make the fluctuations of volatilities more complicated, but the resulting abundance of hard-to-estimate process parameters would make the model difficult to populate, hard for traders to comprehend, and probably not much more flexible in terms of the types of volatility smiles and smirks that could be generated. The process for V is here assumed to be a classical mean-reverting process: dVt(V)(t)d=(t)V-k(tq)d(eZyt + ) ( ) , (1b) where k,,q e >0 , y : + fi + is a smooth function with y(0)=0, and Z is a scalar ¡ ¡ Brownian motion independent of W. V has the role of a scale factor, so typically (but not necessarily) V(0)=q =1. To ensure proper scale behavior, it is most natural to set y(x)=xp, p >0, but for now we do not explicitly impose this choice. A few comments about the chosen framework: • The volatility functions l ,1k,2=,..., N1 - are assumed exogenous and it is up to the k model builder to decide whether a parametric approach or a non-parametric approach should be used in estimating this function. • The assumptions of y(0) =0 and j(0)= 0 ensure that variances and forward rates cannot go below zero. We shall occasionally violate the latter condition to gain tractability. • The V-process generates a near-symmetric smile which is superimposed onto the base smile produced by the function j . The smile generated by the V-process is "stationary" in the sense that the bottom of the smile will move along with fluctuations in the forward rate. • If j()x/ x is a monotonically downward sloping function in x (e.g. j()x,x0=< p <1p ), the base smile is a downward-sloping skew. If j()x/ x is non- monotonic (e.g. j()x,0=1,+x<wp<xp1 q>q ), the base smile can be a true U- shaped smile, but will be non-stationary: when rates move, the bottom of the smile will not move with the forward. • The smile effects of j typically persist for much longer maturities than those of the stochastic volatility process (1b), a consequence of the fact that j introduces dependency in log-increments of forward rate movements. Working with a general skew function j thus allows the model builder some flexibility in shaping the long-term behavior of the smile (but see the point above for a caveat). • While Z is assumed independent1 of W, as long as j()x/ xc„onst . our model nevertheless is able to generate a range of effective correlations between the local forward rate volatility (defined as sl(t,,)F||(Vj)VtF(=)/ F 2 ) and the kkkk k forward itself. 1 The need to assume independence is primarily technical: without it, any change of probability measure introduces terms depending on forward rates in the process for V, thereby destroying the tractability of the model. 3 As shown in Andersen and Brotherton-Ratcliffe (2001), (1a-b) leads to swap par rate dynamics that for many practical choices of j can be closely approximated by ( ) dRtV(tR)(t)t»(d)W()( ) tj l (2) se,,s,esese , where W is a vector Brownian motion under the measure (“swap measure”) induced by s,e using the annuity factor A as numeraire, and where l (t) as usual can be approximated s,e s,e as a linear combination of the l ’s, ks=s+,1,...,e- 1 . Applications and tests of swap rate k approximations such as (2) can be found in numerous sources; see for instance Andersen and Andreasen (2000), Glasserman and Kou (1999), and Rebonato (2001), to name a few. 3. Cap and Swaption Pricing by Asymptotic Expansions. To enable fast model calibration, we now seek efficient means of computing prices of European caps and swaptions. Specifically, consider now a caplet C paying at time T k k+1 the amount CT(F)(T ()=Xd ) - +, as well as a European payer swaption S paying kkk +k1 k s,e at time T the amount STA()()(T(=) RT) K - + (X and K are thus the caplet strike s se,,ssesse s , and the fixed swap coupon, respectively). We have CPkk(k0k)(=T0k,Ed)(FT) X+1 k+1( - )+, SAE(0R)T(=0)( K) As,e ( - )+, se,,sese s , where Ek+1 and EAs,e denote expectations with respect to the probability measures induced by the numeraires P(,t T ) and A (t), respectively. F is a martingale under the former k+1 s,e k measure; R is a martingale under the latter. Notice that with the processes (1a) and (2) s,e having identical form, caplet and swaption price formulas will be completely equivalent (short of numeraire scaling factors), and in this section we only consider the former. Specifically, we here propose a asymptotic expansion for the general caplet pricing problem; in the next section we present an exact transform-based formula for special choices of the function j and y . To develop an asymptotic expansion, consider first the special case where V is non-random, e.g. V ” 1. Even for this simple case, caplet prices cannot be computed explicitly under our process assumption. However, for a large class of j 's, accurate asymptotic expansions can be constructed. For instance, a small-time expansion around the special case j(x) =x results in the following expression: Caplet Expansion for V” 1. Assume V()t =1 for all t. Defining cT=tdt-1(cid:242) Tk||(l)|| 2 and writing k 0 k CP(0T)g(=0F,)(0);c d ( ) for some function gg=F(c; ), we have k+1 k k ln(FX/)F( c,– W)1 2 g(F;c)()(=FF- d)XF d , d = 2 , (3) + - – W (F,c) 4 where F is the cumulative Gaussian distribution function, and W=W+(WFc,)F()c(T)(FcTO)+T1/21/23/23/25/2 , 0 kk 1 k F F I I ln(F/ X) W (F) G G FX J1/2J W (F)= z ; W (F) =- Fz 0 I lnHW (F)H K K . 0 Fj(u)- 1du 1 H Fj(u)- 1duK2 0 j(F)j(X) X X Proof and tests of this result can be found in Andersen and Brotherton-Ratcliffe (2001) who also list the result for the limit F fi X . While being constructed to be accurate for small k maturities, the expansion turns out to often retain its precision for long maturities, even for options with strikes far from at-the-money. This is, for instance, the case for CEV-like specifications j(x)=xp, as well as truncated exponentials, j()x(1= xae+) - bx , a,b> 0. We now relax the assumption that V ” 1 and instead assume that V follows (1b). A small-time expansion for the resulting caplet pricing equations turns out to be cumbersome; instead we adapt the small-e expansions in Hull and White (1987) and Lewis (2000) to our purposes. The result is: Caplet Expansion for Stochastic V. Define YF=ln((0X)/ ) and cV()T||(=)|+| - 1(cid:242)- Tk[tVledt] 2 (q q - kt) . Then k k 0 k ( ) CP(0T)(=0,)(T)(T0);gF(c(0)V) - * with g given as in Eq. (3), and k+1 k+1 k k cV*2(4c)2V(42Y4=4)Y(e)()( +) aOebe6a+ebebe+++ e + L- e2Y2 (4) 0011 2 for coefficients a ,a , b , b , b , and L an arbitrary positive number. The 0 1 0 1 2 corresponding Black-Scholes implied volatility is given by s =W+W cV**(3c)/(V2)( +T)O T2 . impk 0 k 1 The coefficients a , a , b , b , b depend on the parameters of the process of V, 0 1 0 1 2 including the chosen skew function y . Their computation is essentially a matter of tedious algebra, but the resulting expressions are lengthy and for space reasons we here must refer to Andersen and Brotherton-Ratcliffe (2001) for details (article is available on the Internet). Suffice to say that all coefficients can be computed almost instantaneously on a computer, making the resulting expression extremely efficient and fast enough even for production systems pricing and risk-managing many thousands of caplets. The parameter L is, essentially, a defense mechanism protecting the expansion from possible degeneration for options very deeply out-of-the money (where |Y | can grow large). In most cases, the exact value of L is of little importance. 5 We note that for many typical parameters and option contracts, it often suffices to only include terms in (4) to order O(e2) (that is, set b =b =b =0), but the higher- 0 1 2 order terms often become important when, say, the mean reversion speed (k) is low and the volatility-of-variance (e) is high. Figure 1 demonstrates typical performance of the O(e2) and O(e4) expansions (“O4” and “O2” in the figure). Many further tests can be found in Andersen and Brotherton-Ratcliffe (2001). 4. Exact Transform Solution. In cases where very high volatility of variance and/or very low mean reversions are required, the expansions in Section 3 might have to be continued to inconveniently high order. In such situations it is sometimes safer and easier to instead rely on exact transform- based pricing expressions. To develop such expressions, however, we must simplify our setup somewhat. In particular, we specialize to an affine setting with y(x)= x and j()x(m1=x+m )L- , where L‡ 0 and 0<m£ 1. As we shall show, this setup is tractable, although the tractability comes with the cost of permitting interest rates can become negative with positive probability. We also notice that the choice of y(x)= x allows the process (2) to reach 0 with positive probability if, as will often be the case in applications, 2kq <e2. As 0 is non-absorbing, this has few practical ramifications, however. As demonstrated in Section 3, caplet and swaption formulas have the same form in our setup; for variety we will here develop the latter. Dropping subscripts (i.e. R (t) s,e becomes R(t), and so forth), we write SA(0)(=0E-=)(R-)T()K'(E0)TA” (K )+f A(0)(0)A(x )+ A s m s m where Km'(=1K+m L)- and dx(t) ==+Vt(t)m-()ld()W,(t0m)(R0m)(1L x) . x(t) Using the results of Heston (1993) and Lewis (2000) we find the swaption price formula by numerically solving the inverse Fourier integral fH(0)(=0)x(w0d, -)wK'(cid:242) +¥ e()-lniw(+0)12/ 'x K (5) 2p ¥- w2+1 4 where i= - 1 and Ht(, w) =e at(bt;Vw)(;t+) (w) with a,b given as the solutions to the Riccati ODEs da =- kqba, T( ;w=) 0, dt s db =+1- mt2b2l2b2(b),2Tw(ke; )= 01 w . dt 2 2 s 6 For constant l these ODEs can be solved in closed form (see e.g. Lewis (2000) or Lipton (2001)), a result which can also be used iteratively for piecewise constant l. In the general case, a simple numerical Runge-Kutta scheme can be applied at little additional computational expense. We point out that the transform solution can easily accommodate non-zero correlation r between V and R by modifying the term kb in the second ODE to ( ) b k +(iwe-l1r) , but we shall rarely need this as non-zero correlation will hamper the 2 construction of practical models for the joint dynamics of the full yield curve (see footnote 1). We finally note that when computing the integral (5) numerically, performance can often be improved by using the case e =0 as a “control variate” through the split fH(0e)d(=K0--)x(wg0,) K'(cid:242) +¥ e-()-lniw(+0)w12/ 'x 'K ( - (w)/2+2 14 v2 ) 2p ¥- w2+ 1 BS 4 where gx =F+-F()1( )x(0)x( ); BS - K' + xvv=m–=lu+n(x0)/ K-'Vedu1 ,()((220) (cid:242) Ts2) lq ( q - ku) . – v 2 0 With this trick, the combined scheme of Runge-Kutta and numerical Fourier inversion is fast and accurate, although some attention to step sizes is, as always, needed for long-dated options with strikes far from the at-the-money point. Typically, accurate option prices can be obtained in around 0.03 seconds per option on a regular PC (about three times faster for the constant parameter case). 5. A low-dimensional Markov model with stochastic volatility. In Section 2, we specified our stochastic interest rate model as an extension of the multi-dimensional Libor market model. In many cases, it is useful to work with approximating low-dimensional Markov models that allow for option pricing in finite difference grids. We shall here describe such a model which can be made approximately consistent with the cap and swaption expressions derived earlier. As a starting point, let f(t,T)lPn(t,=)T/-¶ T ¶ be the usual continuously compounded forward rate at time t for deposit over the interva l [,TTdT+ ] . As shown by Heath, Jarrow, and Morton (HJM) (1992) any arbitrage free model with a single Brownian motion driving the yield curve level can be written as df (t,T)(t,T)(t=s,s)d(sd)tdWØ (cid:242)tTs + ø Œº œß t { } where W is a scalar Brownian motion under the risk-neutral measure and s(,t T) is a T‡ t collection of general stochastic processes adapted to the Brownian motion. The model is { } fully specified by the initial forward curve and a given volatility structure s(,t T) . The T‡ t 7 general one-factor HJM model requires as full continuum of all forward rates as state variables, with tree or lattice approximations of the model generally non-recombining and largely impractical. Numerically tractable versions of the HJM model were treated in a range of papers in the early 90s, see for example Babbs (1993), Cheyette (1992), Jamshidian (1991), and Ritchken and Sankarasubrahmaniam (1992). In particular, it was noticed that when the forward rate volatility is of the separable form s(,t)T(g)T=(h) t where g is a deterministic function and h is some adapted stochastic process, a finite state variable Markov representation of the yield curve is possible. Specifically, defining n()tg'(t=g)/(-t) and h()tt(tg=,)(st )h(t) = we get the bond reconstitution formula P(,t)T,(e,G=)tTP(0,T)e-dGu(t,T)x(tG)(- t,T)12(y)t 2 =(cid:242) T - (cid:242)tun(sd)s (6) P(0,t) t where, in the risk-neutral measure, dxt(t)(x=)ty-((td)n(t)t(d)W(),t+(0+)x 0 ) h = , (7a) ( ) dyt(t)t()y=2t-d(ht )y(2),(0)n 0 = . (7b) As long as h()tt,()=x,t(hy)(t ), (6) and (7a-b) show that the resulting model is Markovian in only two state variables, x and y. The first state variable x can be interpreted as a yield curve factor, with the instantaneous short rate given by rt(f)(tt”,)(0,f)( )=tx t + . The second state variable y has the interpretation of a convexity correction term ensuring that the model is arbitrage-free. It is worth noting that the model only collapses to a single state variable model in the case when h is deterministic. In this case the model corresponds to the general Gaussian model as presented in, for example, Jamshidian (1991). In traditional implementations of the model above it is customary to let the volatility of the yield factor be a function of the short rate, for example h(t)r((cid:181)t) . p In this specification, however, the volatility skew induced by the coefficient p will tend to flatten for longer underlying tenors. To avoid this we here follow Andreasen (2000) and let h be a function of a longer tenor forward swap rate, i.e. h()tt(,=(Y))ht where Y is some forward swap rate, the identity of which can change over time. Which swap rate is chosen at 8 any point in time will depend on the instrument that we wish to price. For example when pricing a standard Bermudan swaption2 with final maturity T we let e Yt(R)()=,[,t [t T ˛ T , ke,k k -1 whereas when pricing a Bermudan callable cap we could set Yt(R)(tt)=, T[, T[ ˛ . kk,1k+ k 1- In any case, through (6), we would obviously have h(t)t,(=)x,t(hy)(t ). As a general rule, we notice that calibration and specification of low-dimensional Markovian models typically needs to be instrument-specific, unlike the Libor market models which are “large” enough to make possible calibration to the swaption and cap markets as a whole. To incorporate stochastic volatility to our Markov model, we add to (7a-b) the third state variable V as given by (1b). The actual specification of our model is then hl()tV()t=(j)( )tY t ( ), where l is a deterministic function. As before, we will assume that dZd(cid:215) W = 0. For calibration of the model above to the swaption market closed-form or simple approximations for swaption prices are convenient. We note that under our assumptions the forward swap rate R (t) evolves according to s,e ¶ R (t) dRttd(W)()=( )t s,e h se,se ¶ x , where W is a scalar Brownian motion, and s,e (cid:229) e d P(,t T)(G, t) T ¶ R (t) P(,t T)(G,)(,tTP)(,t- T) G t T i- 1 i i s,e = eessi s +R (t) = +1 . ¶ xAtA t ()( ) s,e se,se , To price swaptions we approximate this derivative as being approximately deterministic3, and arrive at ( ) dRtV(t)()»()(tR)(td)Wl* t j , (8) se,,sese , (cid:230) ¶ R (t) j(Y(t)) (cid:246) l*(t)(=) l t (cid:231)(cid:231)Ł ¶s,ex j(R (0))(cid:247)(cid:247)ł . s,e xt(y)(t=) 0 = 2 One might question whether it is sufficient to use only a single driving Brownian motion for the yield curve when pricing Bermudan swaptions. Andersen and Andreasen (2001) conclude that the answer is generally yes, provided that calibration of the mean reversion parameter n is done carefully. 3 To understand why this is reasonable, we notice that ifR were the continuously compounded yield s,e on a zero-coupon bond, it would be exactly linear in x. 9 We find that the approximation above provides sufficient accuracy for most applications, even for long-dated options out to, say, 30 years. Importantly, (8) allows us to use the results in Sections 3 and 4 to efficiently price swaptions and caps. Due to the structure of the { } model, once the mean-reversion of rate n has been set4, the model volatilities l(t) can be bootstrap-calibrated to match a strip of swaption prices. On a standard PC, bootstrap calibration to, say, the swaption strip {1··29,228,,2·91} can be done in less than 5 K seconds when using transform inversion and in a few tenths of second using expansions (for comparison, the calibration can be done in less than one tenth of a second when volatilities are deterministic). Figure 2 demonstrates the quality of the fit for the case of j()x(m1=x+m )L- , m=0.2. 6. Numerical Methods. Numerical implementation of the general Libor market model (1a-b) presented in Section must virtually always be done through Monte Carlo simulation. Due to the zero correlation between the forward rate processes and the stochastic volatility process, it is possible to split the Monte Carlo generation of interest rate paths into two pieces: 1) draw a path of the variance process V through time; 2) draw a path of forward rates assuming that l ()tV( )t is deterministic. Schemes for step 2) are well-known (see e.g. Andersen and k Andreasen (2000)), so we here focus on simulating (1b) on some discrete time-grid {t} . i i=0,1,.. A direct Euler or log-Euler discretization of (1b) is prone to instability unless either the time step or the mean-reversion parameter k are small. Noticing the exact result E(Vt(V)|tV(t)( ) )e=+q ( - q) - k(ti+-1 ti) ii+1 i it is generally better to resort to a moment matching scheme where, for instance, we can approximate Vt(V)|(t) as a log-normal variable: i+1 i ( ( ) ) Vˆt(V)i(t+e1)=+q eˆ- i q - k(ti+-1 ti) -G+G12 ()tt(ii2)zi % , (cid:230)(cid:231) 1ey2 (Vˆt() )2k - 11( -e -2(k ti+-1 )ti )(cid:246)(cid:247) G=()tln2 1+ (cid:231) 2 ( ( i ) )2 (cid:247) , (cid:231) q + Vˆt( )- qe - k(ti+-1 ti) (cid:247) Ł i ł where the z ’s are a sequence of i.i.d. standard Gaussian draws. %i Monte Carlo schemes for the one-factor Markovian model presented in Section 5 are similar to those of the full LM model. Importantly, however, the limited number of state variables (a total of three) allows us to write down a low-dimensional PDE for option prices 4 To prevent the model from being too non-stationary, it is often best to let the rate mean reversion be constant. For instance, we could best-fit n to reproduce the auto-correlation of the swap rates as computed approximately in a globally fitted LM yield curve model; see Andreasen (2000) for details. 10