., thefunctional forms for akt,r) or cp(t,T). Once the volatility structures arespecified, the drift term can be obtained through the above no-arbitragecondition.Heath, Jarrow, and Morton have proposed the following nearproportional model
where W ,(t) and W2(t)are independent one-dimensional Q -Brownianmotions, q and 0 2 are deterministic yield volatilities (i.e., independent ofany stochastic variables), M is a large constant positive number, and themin(.) function returns the minimum of its arguments.This model is similar to a lognormal model in the sense that it insures positive interest rate (positive forward rates in this case). The useof minHf(t,T),M) is to prevent At,r) from diverging with positive probability. As long as it is much greater than the initial forward rate, theexact value of M i s not critical and can be chosen to be, e.g., 10 standarddeviations away from the initial forward rate at the longest expiration ofthe derivatives being modeled. Jarrow used M = 1,000,000 in one of hisexamples (Jarrow, 1996).Typically, 01 and o2are obtained by calibrating the model to pricesor implied volatilities of selected caplets/floorlets and European swaptions.Generally, q (assuming it is associated with the most dominant factor) is roughly constant, and reflects the change in the level of the entireyield curve (forward rate curve in this case), and 0 2 (assuming it is associated with the second dominant factor) is roughly a function of T-t, andreflects the change in the slope of the yield curve. A third factor can alsobe added to reflect the change in the curvature of the yield curve.
The Heath-Jarrow-MortonFramework
222
Chapter 6
The general forms of the HJM models are non-MarkovianY8'andthus, in general, cannot be solved by the computationally efficient recombining lattice and PDE grid techniques. The conventional bushy treetechnique and the conventional Monte Carlo technique that can be usedto solve the general forms of the HJM models suffer from various technicaVnumerica1 difficulties, such as the computation time limitation. Thereare significant efforts in searching for efficient numerical algorithms forsolving the general forms of the HJM models.Attacking this problem from a different angle, some researchers haveidentified special volatility structures, under which the HJM models areMarkovian with respect to a finite number of state variables (see, e.g.,Ritchken and Sankarasubramanian, 1995, and Inui and Kijima, 1998),and thus can be solved by the computationally efficient recombining lattice technique. We review these approaches in the next two sections.One of the important state variables for these approaches is the shortrate. So let us take a look at the short rate process first. Since the shortrate is defined asr(t)=f(t,t)
(6.15)
Y
thus,r ( t )= f (0, t ) -
6
(
( T p ( S , t )0; S , t ) ds
+
6
of(s,t)dW, .
(6.16)
Differentiating both sides of the above equation yields the followingshort rate proced2
XI
To be more precise, the general forms of the HJM models are non-Markovian withrespect to the underlying Brownian motion generated states, but can be Markovian withrespect to all the states (i.e., all the forward rates).The differentiation of the short rate r(t) is quite different from that of the forward ratef(t,7), as in the latter case the second argument Tremains as constant.
Sec. 6.1 The Heath-Jarrow-Morton Model
223
d r ( t ) =-af ( 0 4 dt-op(t,t)o; ( t , t ) d t + o r (t,t)dW,at
Recall that o p ( t , t ) = O and
aoP( S , t ) = -or (s,t) . Then, theat
above short rate process can be rewritten as
af (atO ~ % f t + o(t,t)dW,,
d r ( t ) =-
1
(6.18)
which, except for some special cases, is non-Markovian due to the history or path dependence introduced by the integrals in [ ]. For an nfactor process, o,(t,T), okt,T),W, are all n-dimensional vectors.It is interesting to point out some special cases or special volatilitystructures under which the HJM short rate process is Markovian. Morespecifically, it is easy to prove that the HJM model can be reduced to theHo-Lee and the Hull-White model (the extended Vasicek model) with
of( t , T ) = o o and of(t,T)=ooexp(-a(T-t)), respectively, whereq,and a are positive constants.Exercise 6.1. Prove the above statement.Ritchken and Sankarasubramanian (1995) and Inui and Kijima(1998) developed approaches for reducing the HJM models to Markovian processes with respect to a finite number of state variables, which
224
The Heath-Jarrow-Morton Framework
Chapter 6
we shall discuss in the following sections. Ritchken and Sun (2003) alsoextended the same methodology to the modeling of risky interest rates orcredit spreads.
6.2
The Ritchken-Sankarasubramanian Model
Ritchken and Sankarasubramanian (1995) have identified a class of volatility structures, under which the 1-factor HJM models are Markovianwith respect to two state variables, the short rate r(t) and the cumulativevariance Kt). The Ritchken-Sankarasubramanian (RS) model can besolved by the computationally efficient recombining lattice technique(Li, Ritchken, and Sankarasubramanian, 1995).83The basic idea is to identify a particular volatility structures aAt,T) sothat a minimum number of Markovian state variables (including the shortrate itself) can be used to capture the history or path dependency.Recall that the HJM short rate process is given by
(6.19)
and
r ( t ) =f (0,t)- l ( T p ( s , I ) o > ( s , t ) d s +c a f ( s , t ) d W s .
(6.20)
Now, our goal is to reduce to a simpler form some of the history orpath dependent terms in the HJM short rate process, i.e.,
*'
In the original LRS paper, a recombining binomial tree was used for the numericalimplementation. However, a recombining binomial tree when solving the RS modeloften times does not converge. A recombining trinomial tree, on the other hand, canprovide much better accuracy and convergence.
Sec. 6.2
The Ritchken-Sankarasubramanian Model
225
Noting the resemblance of these terms with the short rate r(t) itself,we restrict the volatility structures oAt,T) to the following special functional form(6.22)or, equivalently,(6.23)where K ( t ) is a scalar. In other words,
of( t , T ) = O f ( t J ) k ( t , T ) ,
(6.24)
k(t,T)=eXp(-~-K(s)ds),
(6.25)
where
Thus, these history or path dependent terms reduce to a simple functionof r(t) as follows (by plugging Eq. (6.20) in the last step)
= --K ( t )( r (t ) -
f (0,t ))
Consequently, by plugging the above equation to Eq. (6.19), theHJM short rate process reduces to the following RS model
(6.27)
where &t) is the cumulative variance given by(6.28)
226
Chapter 6
The Heath-Jarrow-MortonFramework
Differentiating the above equation yields
which is a locally deterministic or predictable process due to the lack ofany incremental Brownian motion term.Thus, in general, the 1-factor RS model is Markovian with respect totwo state variables, r(t) and Kt). If of( t , t ) is state independent, i.e.,does not depend on ~ ( tand) Kf),and is not stochastic by itself, then f i t )is deterministic. Consequently, in this case, the 1-factor RS model is a 1factor normal model and is Markovian with respect to one state variable,
40
*
As always, if we can find a closed-form solution for the zero-couponbond price P(t,T) or the forward rate At,r), it will significantly improvethe computation efficiency for derivatives pricing, as the payoffs of thederivatives generally depend on the future zero-coupon bond prices andforward rates, or the future yield curve. If we can further find closedform solutions or approximations for options, such as caps/floors andEuropean swaptions, then we can also speed up the volatility calibration.Such approximations, though maybe crude, are typically available (see,e g , Section 7.6.1).As shown in Appendix 6.5.1, the closed-form solution for the zerocoupon bond price P(t,T) or the forward rateAt,T) is given by
and
>
where
(6.31)
Sec. 6.2 The Ritchken-SankarasubramanianModel
P(t,T)=r k ( t , u ) d u .
227
(6.32)
Another way of intuitively understanding the RS volatility structurepresented in Eq. (6.24) is to rewrite it as follows taking into account Eq.(6.25), i.e.,0
tt
( ’ ) k(O,T),o tT( ’ ) = k(0,t)
(6.33)
or
where o(t) is an n-dimensional vector and g(T) is a scalar. Essentially,the time and the term dependency of the volatility are separable.So far, on the surface, it is still an n-factor model. We can show thatit actually can be reduced to a 1-factor model. More generally, any HJMmodel with the above separable volatility structure can be reduced to a 1 factor model (under proper technical conditions).
Exercise 6.2. Prove the above statement. What if o(t) is a scalar andg( T ) is an n-dimensional vector?
of( t , t ) can actually take quite general functional forms, i.e.,Of( t J )= O ( r ( t ) , @ ( t ) J* )
(6.35)
For practical applications, it is often chosen to be a CEV function form,i.e.,Of( t , t ) = O ( t ) r a ( t ) ,
(6.36)
where 0 I a I 1 and a(t ) is state independent, i.e., does not depend on
r(t) and K t ) , and is not stochastic by itself.If a=O,then the RS model reduces to the Hull-White model. Ifa = 0 . 5 , then RS model is similar to the CIR square root model. Ifa = 1, then RS model is similar to the BDT model.Alternatively, we also use the following similar shifted lognormal
228
The Heath-JarrowMorton Framework
Chapter 6
volatility structure,Of( t J )= O ( t ) (a r ( t ) + (1- a )f ( 0 , t ) )'
(6.37)
Exercise 6.3. Prove that if a=O,then the RS model reduces to theHull-White model.As we can see that the RS model allows rich volatility structures.Generally, a is chosen (or calibrated) to match or best fit certain volatility skews. However, it does preclude certain intuitively appealing volatility structures, such as those depending on the level of the forwardrates, which are allowed by the general HJM model even with one factor.In other words, in the RS model, the volatilities of the forward rates aredriven by the level of the short rate.The RS model can indeed produce the commonly observed humpedvolatility term structures of the forward rate f ( t , T ), which are increasing functions of the time-to-maturity when the time-to-maturity is short(e.g., less than a few years) and decreasing functions of the time-tomaturity when the time-to-maturity is long. According to Eq. (6.24) and(6.25), this can be achieved with a slightly positive ~ ( t in) the shorttime horizon and a slightly negative ~ ( t in) the longer time horizon.According to Eq. (6.27), ~ ( t is) also the mean reversion parameter forthe short rate. More specifically, a positive K( t ) indicates mean diverging and a negative K( t ) mean reversion.Typically, ~ ( tand) ~ ( t are) obtained by calibrating the model toprices or implied volatilities of selected caplets/floorlets and Europeanswaptions.
6.3
The Inui-Kijima Model
The h i - K i j i m a (Inui and Kijima, 1998) model is an extension of the RSmodel to multi-factors. It is an n-factor HJM model that is Markovianwith respect to 2n state variables.The IK (Inui-Kijima) model can be derived in a very similar ap-
Sec. 6.3
The Inui-Kijima Model
229
proach as the RS model. Recall that the HJM short rate process is givenbyd r ( t )=- af(07t)dt+af(t,t)dW;at
(6.38)
andr ( t ) = f ( O , t ) - [o,(s,t)a>(s,t)ds+
[of(s,t)dW, ,
(6.39)
where aP(s7t),o,(s,t), and W, are all n-dimensional vectors.Let ob ( s , t ) , q(s,t), and Wi(s)denote the ifh component of op(s,t),oAs,t), and W,, andxi(t)--ga;,(s,t)oi(s7t)ds+
g q ( s , t ) d W ,( t )
(6.40)
(I 2 i 2 n )Then, the short rate in Eq. (6.39) can be rewritten asr(t)=J.(o't)+~xi(t).
(6.41)
i=l
Differentiating Eq. (6.40) yields the dynamics for xi ( t ) , i.e.,dXi ( t ) = oj(t,t)dW,( t )
(6.42)
Like before, our goal now is to reduce to simpler forms the two ofthe history or path dependent terms in the process of xi ( t ) , i.e.,
The Heath-Jarrow-MortonFramework
230
Chapter 6
Noting the resemblance of these terms with xi(t) itself, we restrict thevolatility structures q ( t , g to the following special functional form
(6.44)or, equivalentlyaoj
(t , T ) = - K ~( T ) o i( t , T ) ,
dT
(6.45)
where K~ ( t ) is a scalar. In other words,oj ( t ,T ) = oj (t,t)kj ( t , T ),
(6.46)
where
Thus, these history or path dependent terms reduce to a simple functionof as follows
=-Ki(t)Xi(t)dt
Consequently, by plugging the above equation into Eq. (6.42), weget the following IK modela!xi (t ) = (-Ki (t )Xi (t ) + @i (t ))dt + Oi (t ,t ) d p (t ) ,
(6.49)
where 4;(t)is the cumulative variance given by@i ( t )=
10:
(s,t)ds .
Differentiating the above equation yields
(6.50)
Set. 6.3
The Inui-Kijima Model
231
d4j ( t >= (0;(t,t>-2Kj (t>4j(1,)dt
9
(6.51)
which is a locally deterministic or predictable process due to the lack ofthe dWi(t) (incremental Brownian motion) terms.Thus, in general, an n-factor IK model is an n-factor HJM model thatis Markovian with respect to 2n state variables, xi(t) and #i(t). If cri ( t ,f)is state independent, i.e., does not depend on xi(t) and &(t), and is notstochastic by itself, then A(t) is deterministic. Consequently, in this case,the n-factor IK model is an n-factor normal model and is Markovian withrespect to n state variables, xi(t).As always, if we can find a closed-form solution for the zero-couponbond price P(t,T) or the forward rateAt,T), it will significantly improvethe computation efficiency for derivatives pricing, as the payoffs of thederivatives generally depend on the future zero-coupon bond prices andforward rates, or the future yield curve. If we can further find closedform solutions or approximations for options, such as caps/floors andEuropean swaptions, then we can also speed up the volatility calibration.Such approximations, though maybe crude, are typically available (see,e.g., Section 7.6.1).As shown in Appendix 6.5.2 and similar to the case of the RS model,the closed-form solution for the zero-coupon bond price P(t,T) or theforward ratef(t,r) is given by
i=l
and
(6.52)
232
The Heath-Jarrow-Morton Framework
Chapter 6
where
1T
Pj ( t , ~=) kj (t,u>du
(6.54)
*
We can gain more intuitions, if we rewrite Eq. (6.52) as followsf ( t , T )= f ( 0 , T )
where we have assumed that &(t) is the ithmost important driving factor(or corresponding to the thmost important principal components) of theyield curve. Similar to the RS model, k, ( t , T ) is closely related to thevolatility structure (or hump) of the forward rate f ( t , T ) and the meanreversion (or diverging) parameter of the short rate. The ratiosk; ( t ,T )( i = 2,3, ...,n ) largely determine the functional form of thek, ( t ,T )higher order principal components of the yield curve, such as the slopeand the curvature.Similar to the RS model, oj( t , t ) can actually take quite generalfunctional forms. It has been suggested to be a CEV function form, i.e.,o;(t,t)=oi(t)ra'(t),
(6.56)
ajI1 .g4 Generally, a; arewhere oj( t ) is state independent and 0 Ichosen (or calibrated) to match or best fit certain volatility skews.For better numerical computation efficiency, we can identify furtherreduced functional form for the volatility structure oj( t , t ). In other
words, if we restrict oi( t , t ) to only depend on xi(t) and &t), but notother state variables, i.e.,
84 There may be technical difficulties if
1 < ai or ai <
t in that the solution may be
unbounded or non-unique. However, these technical difficulties can be easily circumvented by setting an upper limit to the volatility and the drift, which we shall discuss lateron.
Sec. 6.3
The h i - K i j i m a Model
0;
233
(6.57)
( t ,t ) = oj ( t ,xj 4j)
then, the process of each pair of state variables xi(t) and 4i(t) are completely independent of any other pairs of state variables. Thus, in thiscase, a multi-factor IK model can be solved by solving independent 1factor models and performing a convolution. This significantly simplifies the numerical implementations, particularly, recombining lattice implementations.Some feasible functional forms of oi(f,f) in this category are, in theCEV form,(6.58)
o;(f,t)=oi(t)(Y;(t))aiY
or, in the shifted lognormal form,o j ( t , t ) = o j ( t ) ( a j y j(t)+(1-aj)s(O,t))
3
(6.59)
where(6.60)andn
x u i ( t ) = 1*
(6.61)
i=l
In other words,
(6.62)= t Y i(t)i=l
Another advantage of the above further reduced volatility structure isthat if all the yi(t)'s are positive, then the short rate r(t) is also positive.With proper parameterization, the short rate r(t) can also reach any givenpositive number however small it is.
234
6.4
The Heath-Jarrow-Morton Framework
Chapter 6
Overview of Numerical Implementations of the RS andthe IK Model
Solving an n-factor model that is Markovian with respect to 2n statevariables using recombining lattices or trees is not a trivial task. In thissection, we shall present recombining tree techniques that are suitable forsolving both the RS model and the IK model.Li, Ritchken, and Sankarasubramanian (LRS, 1995) solved RSmodel (1-factor Markovian with respect to 2 state variables) by the computationally efficient recombining tree technique. The basic idea is tobuild a tree for the process of the short rate r(t) and approximate theprocess of the cumulative variance &t). The rationale behind this idea isthat the cumulative variance &t) is a locally deterministic or a predictable process.In the original LRS paper, recombining binomial tree was used forthe numerical implementation.In general, however, recombining trinomial trees can provide muchbetter accuracy and convergence. The basic reason is that a recombiningtrinomial tree has more degrees of freedom than a recombining binomialtree. Consequently, in the case of modeling a stochastic process withstate-independent (but non-stochastic) parameters (i.e., drift and volatility), a recombining binomial tree can match the local mean, variance, andall odd order moments, whereas a recombining trinomial tree can matchthe local kurtosis (the fourth order moment), in addition to matching thelocal mean, variance, and all odd order moments.In the case of modeling a diffusion process with stochastic and statedependent parameters, a recombining trinomial tree can match the localmean and variance, but a recombining binomial tree may not. A recombining binomial tree (or any other types of trees) not matching the localmean and variance in general does not converge to the desired continuous limit. In the case of modeling a diffusion process with stochastic andstate-dependent parameters, a recombining trinomial tree generally cannot match the local skew and kurtosis. In this case, the trinomial treeparameters with best approximation to the local skew and kurtosis shouldbe used for better convergence.
6.4.1
Recombining Trinomial Tree Technique
Here we shall first give a brief introduction to the recombining trinomialtree technique.
Sec. 6.4
Overview of Numerical Implementations of the RS a n d the IK model
235
In a trinomial tree, there are three branches emanating from each andevery tree node.85 In a recombining trinomial tree, most of the branchesrecombine or share the same end points with other branches.86 Figure6.1 shows schematic diagrams of a recombining trinomial tree. Asshown in Figure 6.l(a), from any given tree node X(t), there are threebranches, up (u), middle (m),and down branch (4, leading to the nextstep at positions of X,(t+dt), X,(t+dt), and XAt+dt). pu,pm,and P d denote the branching or transition probabilities.(a)
<
(b)
Xu(t+At)
X(t)
X,,,(t+At)
X&+At)m=l
m=2
m==4
Figure 6. I . Schematic diagrams of a 1-step trinomial tree (a), which serves as thebuilding block for multi-step recombining trinomial trees, a 2-step recombiningtrinomial tree (b), and a 4-step recombining trinomial tree (c). As shown in (a),from any given tree node X(t), there are three branches, up (u), middle (m), anddown branch (d), leading to the next step at positions of Xu(t+At), Xm(t+At), andXd(t+At). pu, pm, and pd denote the branching or transition probabilities.
As a summary, in order for a recombining trinomial tree to accurately model a stochastic process driven by a one-dimensional Brownianmotion, the following conditions need to be satisfied simultaneously.1. The local or conditional mean and variance on the recombiningtrinomial tree need to match those of the stochastic process beingmodeled. For better convergence, one can have the local or conditional skew and kurtosis match or be close to those of the stoIf the number of branches emanating from each and every tree node is two, then it is abinomial tree. The number of branches can be more than three and can vary from node tonode, in which case, it is a multi-nomial tree.86 The recombination of the tree branches is for reducing the number of tree nodes and,thus, better computation efficiency. It is not necessary, nor possible, to have all thebranches recombine. For instance, the outermost tree branches may not recombine withany other branches.
236
The Heath-Jarrow-Morton Framework
Chapter 6
chastic process being modeled.2. All the branching or transition probabilities need to be nonnegative, as negative branching or transition probabilities indicate the existence of arbitrage opportunities.3. Most of the branches recombine or share the same end pointswith other branches.It is trivial to have each of the above conditions satisfied individually. Nonetheless, it can be challenging to have the above conditionssatisfied simultaneously when the parameters are state-dependent (suchas in the case of interest rate modeling).The third condition essentially is for reducing the number of treenodes and, thus, for better computation efficiency. We can see thisthrough the total number of tree nodes. For the recombining trinomialtrees shown in Figure 6.1, the number of tree nodes at the mth tree step isgiven by 2m + 1, or, more generally,(ng - I ) m + l ,
(6.63)
where nB is the number of branches emanating from each node of thetree. The total number of tree nodes for m steps is thus given by
On the other hand, if the tree is completely non-recombining orbushy, then the total number of tree nodes for m steps is thus given by(6.65)Let's consider an example of a trinomial tree with 30 steps. If thetree is a recombining trinomial tree as shown in Figure 6.1, then the totalnumber of tree nodes is 960 as given by Eq. (6.64). If the tree is bushy,then the total number of tree nodes is 330z 2.06 x lOI4 as given by Eq.(6.65). The huge number of nodes in a bushy tree presents significantchallenges to the current computation technologie~.~~87 Suppose that one has a computer that is capable of computing 10 million tree nodes per
second, then computing all the 2.06 x 1 OI4 nodes would take about 0.65 year! There aresome techniques, such as the nonexploding bushy tree (NBT) techniques, that can help tocircumvent such difficulties. The nonexploding bushy tree (NBT) techniques shall bediscussed later on.
Sec. 6.4
Overview of Numerical Implementations of the RS and the IK model
231
A general way of insuring the recombining of the tree is to establisha grid with pre-specified boundaries and to have the end points of thetree branches land only on the grid points, but not necessarily on the gridpoints adjacent to that of the center branch. This is illustrated in Figure
Figure 6.2. Examples of branching of the tree branches of a recombining trinomial tree on a grid.
Various branching configurations illustrated in Figure 6.2 are neededfor handling, e.g., state and time dependent drift (such as mean reversion) and volatility, which are common in the interest rate modeling.The number of tree nodes in a tree involving these general branchingconfigurations can be more or less than what is given in Eq. (6.64).The exact number of the tree branches depends on the conditions formatching the local mean and variance (and approximately matching thelocal kurtosis), and non-negative transition probabilities. These conditions can be explained with respect to the following general Brownianmotion driven stochastic process,
dx ( t ) = p ( t )at + o(t )aw ( t ),
(6.66)
where p ( t ) and o ( t ) are the drift and the volatility, respectively, bothof which can be time- and state-dependent, W ( t ) is a standard onedimensional Brownian motion. In discrete time, the above equation reduces to
AX(t ) = x(t+ at)- x(t ) = p ( t)at + +)AFT(
t ),
(6.67)
Based on the above equation, we can translate the condition formatching the local mean and variance into the following equations, i.e.,
238
The Heath-Jarrow-Morton Framework
Chapter 6
[
IVur, [ A X ( t ) ]= Et ( A X ( t )- Et [ A X ( t ) ] ) 2 ]= ~ ( tAt2) ,~ (6.68)
where E,[.] and Vur,[.]are the expected value operator and the varianceoperator, both conditioned on the information available up to time t andunder the same measure under which the standard one-dimensionalBrownian motion W ( t ) is defined. Note that in deriving the aboveequations, the following properties of the one-dimensional incrementalBrownian motion were used, i.e.,El [ A W ( t ) ] = 0
Vur,[ A W ( f ) ]= At2 '
(6.69)
Using Eq. (6.68) and the symbols in Figure 6.l(a), we can furthertranslate the condition for matching the local mean and variance into thefollowing equations, i.e.,
where
AT, ( t )= x, ( t + ~ t-(x(t))+p ( t ) ~ t )( a = d,rn,u)
(6.71)
The condition for non-negative transition probabilities can be simply expressed asg8Pd 2 0, Pm2 0, Pu 2 0 .
(6.72)
With proper tree construction, solutions satisfying both Eq. (6.70)and Eq. (6.72) exist. These solutions are generally not unique, and thosethat best match the local kurtosis can be selected.88 The condition for non-negative transition probabilities also ensures that the probabilities are no greater than one, as all the probabilities add up to one.
Sec. 6.4
6.4.2
Overview of Numerical Implementations of the RS and the IK model
239
Adaptive Recombining Trinomial Tree Technique
Here we shall introduce an adaptive technique for the tree constructionsuch that the solutions to equations of a recombining tree presented in theprevious section can be easily obtained.The basic idea is that the grid for the tree is constructed adaptively ordynamically (as opposed to being pre-specified). In other words, the gridfor the tree is formed by constructing the middle tree branches exactlyfollowing the (local) drift and the outermost tree branches flowing the(local) volatility at all time steps. More specifically, the grid for the treeis constructed according to the following equation,
1
S J t ) =0
(at all grid points)
= k, Io(t)lAt
(at max(X(t)))
ad(t)=