Zero Bound, Option-Implied PDFs, and Term Structure Models

Finance and Economics Discussion Series Divisions of Research & Statistics and Monetary Affairs Federal Reserve Board, Washington, D.C. Zero Bound, Option-Implied PDFs, and Term Structure Models Don H. Kim 2008-31 NOTE: Staff working papers in the Finance and Economics Discussion Series (FEDS) are preliminary materials circulated to stimulate discussion and critical comment. The analysis and conclusions set forth are those of the authors and do not indicate concurrence by other members of the research staff or the Board of Governors. References in publications to the Finance and Economics Discussion Series (other than acknowledgement) should be cleared with the author(s) to protect the tentative character of these papers.

ZERO BOUND, OPTION-IMPLIED PDFS, AND TERM STRUCTURE MODELS ∗ DON H. KIM Abstract. This paper points out that several known ways of modeling non-negative nominal interest rates lead to different implications for the risk-neutral distribution of the short rate that can be checked with options data. In particular, Black’s boundary models (“interest rates as options”) imply a probability density function (pdf) that contains a Dirac delta function and a cumulative distribution function (cdf) that is nonzeroatthezeroboundary[P(K) KbehaviorforputoptionpriceforstrikeK near ∝ 0], while models like the CIR and positive-definite quadratic-Gaussian (QG) models haveazerocdfattheboundary[P(K) Kα, α>1]. Eurodollarfuturesoptionsdata ∝ are found to favor Black’s boundary models: the CIR/QG models, even multifactor versions,havedifficulty capturingoptionpricesaccuratelynotonlyinlowinterestrate environments but also in higher interest rate environments, and data in early 2008 provide an almost tangible signature of the Dirac delta function in Black’s boundary pdf models. Options data also contradict the prediction of well-known models whose cdf is zero at the zero boundary, namely that the risk-neutral pdf is always positively skewed. 1. introduction It is well known, at least since Breeden and Litzenberger (1978), that options at a broad range of strikes can provide information about the whole risk-neutral distribution of the underlying security prices, not just mean and variance. Policy makers and market participants have utilized this fact to deduce the risk-neutral distribution of short-term interest rates from eurodollar futures options or federal funds futures options, and many studies have explored techniques that are useful in this regard and discussed their applicationstovarioussettings.1 However, thisliteraturesofarhasnotmademuchconnection with the vast literature on dynamic term structure models, i.e., relatively little work has been done to investigate what the option-implied risk-neutral distribution tells about Date: June 12, 2008. ∗ Division of Monetary Affairs, Federal Reserve Board. email: Don.H.Kim@frb.gov. I thank Alain Chaboud, Jim Clouse, Darrell Duffie, Benson Durham, Young-Ho Eom, Tony Rodrigues, Jonathan Wright, and the seminar participants at Yonsei University for helpful comments, and Dion Chiu and Scott Konzem for help with eurodollar futures options data. The views expressed in this paper do not necessarily reflect those of the Federal Reserve Board. 1See, for example, the papers in BIS (1999). 1

2 D. H.KIM the underlying stochastic process of interest rates and whether known term structure models can capture the risk neutral distribution implicit in option prices.2 The present paper is a contribution toward filling this gap. We shall be focusing in particular on the the zero boundary behavior of the short rate and its ramifications for term structure models. Perhaps the best known model that recognizes the presence of the zero bound is the CIR model, written down by Cox, Ingersoll, and Ross in their seminal paper (1985). This type of model has the feature that the cumulative distribution function (cdf) of the short rate vanishes as the short rate approaches zero. A different treatment of zero-boundary behavior was proposed by Fischer Black in his last paper (1995), in which the nominal short rate is modeled as r = max[x ,0], where x is a “shadow-rate” process that can go below zero. As we shall t t t see, this type of model implies a short rate cdf that is nonzero at the zero boundary. Proper modeling of the behavior near the zero boundary (e.g., whether it is better described as CIR-like, Black-like, or something else) is a basic and historically important problem in term structure modeling, but it has significance beyond that, as different boundary behaviors could signify different underlying macroeconomic mechanism of the interest rate determination. In the case of Japan, studies by Gorovoi and Linetsky (2004) and Ueno, Baba, and Sakurai (2006) using yields data give strong support for Black-type models. However, some might argue that this conclusion does not necessarily carry over to the postwar US economy.3 Indeed, non-negative interest rate modeling with US data has so far been predominantly in terms of models whose short rate cdf is zero at the zero boundary, e.g., the BGM model (Brace, Gutarek, and Musiela (1997)), the multifactor CIR model, the quadratic-Gaussian (QG) model, and the non-affine model of Ahn and Gao (1999).4 As the US short rate in the past 50 years has not been low enough to see a significant effect of the zero bound on the yield curve,5 the distinction between the two types of models is difficult to discern from yields data. A key point of the present paper is that data on out-of-money options can help, as the models may imply strong enough differences in risk-neutral distributions, especially near the zero bound. Recently (early 2008), shortterm interest rates have come down to a fairly low level while uncertainty about interest rates has remained relatively high, thus it is particularly interesting to explore what the 2Ho˝rdahl(2000)and Ho˝rdahland Vestin (2005)explore the risk-neutralandphysicalmeasure probability density functions of certain term structure models, but they do not apply the models to actual options data. 3In the US, interest rates seemed to display a Black-like behavior during the Great Depression of 1930s,whichmayhavebeenoneofthemotivationsforBlack’spaper(1995,seep1376). Theanalysisof the yield curve during that period is complicated by certain institutional effects; see Cecchetti (1988). 4Some results in the literature could be construed as an indication of reflection-like behavior. For example, Goldstein and Kierstead (1997) note that according to Ait-Sahalia (1996) the US short rate in the 1973-95 period had (physical-measure) dynamics that was similar to a Brownian motion with two reflecting barriers. 5Bomfim (2003) notes that even when the federal funds target rate was at the historical low of 1% in 2003, the effect of the zero bound was not very tangible in the yield curve.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 3 options data from this episode tell about term structure models. Beyond fundamental theoretical interest, modeling of risk-neutral distribution in such environment is an important practical problem for market participants as well. This paper makes largely two contributions. The first is to explore the implications of well-known and tractable term structure models that respect the non-negativity condition for nominal interest rates (e.g., the multifactor positive-definite CIR/QG models) for the risk-neutral distribution of the short rate. In particular, this paper derives the asymptotic form of their risk-neutral probability density function (pdf) at vanishing interest rates (relevant for the discussion of the behavior near the zero boundary), which does not seem to have been explored in the literature. This paper also develops a method for fast computation of option prices for the risk-neutral distribution implied by the multifactor CIR model and (positive-definite) QG model, and explores how well these models perform in matching the observed out-of-money eurodollar futures option prices. In addition, this paper points out that known models which have zero cdf at the zero boundary imply a fairly strong restriction on the shape of the distribution, namely that the distribution is always positively skewed, and examines whether this prediction is born out in the options data. The second contribution is to examine the implications of Black-type boundary models for option prices and option-implied pdfs. After deriving the risk-neutral pdf of the simplest Black’s boundary model (“Black-Vasicek model”), this paper develops flexible parametric models of risk-neutral pdf that are consistent with Black’s boundary behavior and lead to convenient calculation of option prices, and explore their practical performance. The main findings are as follows. A flexible parametric pdf model with Black’s boundary behavior (normal-mixture shadow rate model) is found to perform quite well in matching the option prices across various strikes in a variety of interest rate environment from 1998 to 2008. The model produces implied prices of the 1-year-maturity 0.5%-strike eurodollar futures option in February and March of 2008 that agree fairly well with actualsettlement prices. Ontheother hand, the traditionallognormal-mixture pdfmodeland2-factorand3-factor(positive-definite) CIR/QGmodelssubstantially underprice this “low-strike” put option. The CIR/QG models have difficulty matching the option prices accurately not only in low interest rate environments but also in higher interest rate environments, often generating large and biased pricing errors. In addition, options data contradict these models’ prediction that the skewness of the pdf is always positive. The remainder of this paper is organized as follows. Section 2 sets up the stage for the later parts of the paper, discussing empirically important zero boundary behaviors and reviewing the risk-neutral pdf extraction techniques. It also introduces a mathematical object called the Dirac δ-function, which seldom appears in finance literature but is neededforthediscussionofthepdfswithBlack’sboundarybehavior. Section3examines the pdfsandoptionprices inone-factor models, comparing Black’s boundary modelwith CIR/QG models. Section 4 examines richer pdf models and develops useful techniques

4 D. H.KIM for them; the applicationof these models to actualoptionsdata and theempirical results are discussed in Section 5. Section 6 concludes. 2. Preliminary discussion 2.1. “Zero-cdf model” vs Black’s model. Nominal interest rates cannot go below zero. The best known model that respects this condition is the CIR model. This model is nested by the so-called CEV model6 (1) dr = κ(θ r )dt+σrαdB , t − t t t which admits two kinds of boundary behavior: when α > 1/2, the short rate r does t not reach zero, i.e., zero is an inaccessible boundary (also called unattainable boundary). When α < 1/2, the short rate r can reach zero; in this case, zero is an accessible t boundary. The CIR model (α = 1/2) can display both behaviors depending on the parameters. More specifically, if the Feller condition (2) kθ/σ2 > 1/2 is satisfied, the zero boundary is inaccessible, and if not, it is accessible. In the accessible boundary case (α < 1/2, or α = 1/2 and kθ/σ2 < 1/2), zero is a regular boundary (according to the terminology of Feller (1952)),7 and additional conditions at the boundary can be imposed. For example, requiring that r stay at zero t forever after it hits the zero boundary amounts to having an absorbing boundary. As the permanent stay of the short rate at zero would not be a promising description of the actual economy, we shall not consider the absorbing boundary scenario further in this paper, and instead focus on the regular “unrestricted boundary behavior”, after Longstaff (1992).8 Most of the non-negative interest rate models applied to the US data have been either inaccessible boundary models or regular unrestricted boundary models. The former include the CIR model with kθ/σ2 > 1/2, the BGM model, and Ahn and Gao (1999)’s non-affine model. The latter include the CIR model with kθ/σ2 < 1/2 and (positivedefinite) quadratic-Gaussian model (Beaglehole and Tenney (1992)). All these models (i.e., both the inaccessible boundary models and regular unrestricted boundary models) share the feature that the cdf of the conditional distribution of the short rate is zero at the zero boundary. Therefore, in this paper we shall refer to these models collectively as “zero-cdf models”, for brevity. 6This model (in the physical measure) was empirically investigated by Chan et al (1992), who obtained an α estimate of 1.5. 7See, also Karlin and Taylor (1981, Chap 15). 8Longstaff (1992) calls the bond pricing for the CIR model without further conditions “the unrestricted equilibrium”. One could consider yet other types of boundary condition, such as those leading to “sticky boundary” or “elastic boundary”, but these yield much less tractable models, and I am not aware of any studies in the finance literature that have explicitly investigated these possibilities.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 5 Let us now discuss Black’s boundary models. As the short rate in these models can stay at the zero boundary for an extended period, they have the feature that a finite probability mass is concentrated at the point r = 0; this means that the probability density function of the conditional distribution, f(r), is infinite at r = 0. Mathematically, ǫ (3) Prob(r = 0) = lim f(r)dr = w > 0, f(r = 0) = . ǫ→0+ Z−ǫ ⇒ ∞ A natural way to describe such a distribution is to use the Dirac δ-function, which is defined as9 ∞ 0 if x = 0 (4) δ(x) = 6 δ(x)dx = 1. if x = 0, (cid:26)∞ Z−∞ ∞ It has the property h(x)δ(x a)dx = h(a), where h( ) is an arbitrary function, and −∞ − · it can be thought of as the pdf of a normal distribution (or any other distribution) with R a vanishing variance, i.e., 1 (x a)2 (5) δ(x a) = lim exp − . − ǫ→0 √2πǫ − 2ǫ2 (cid:18) (cid:19) Using the δ-function, the pdf of Black’s boundary models can be expressed as10 ˜ (6) f(r) = wδ(r)+f(r) ˜ where f(r) is an improper probability density function with support at [0, ] (i.e., ∞ ˜ ∞ ˜ ∞ f(r)dr = f(r)dr = 1 w < 1). −∞ 0 − r The cumulative distribution function of the short rate, F(r) = f(s)ds, takes the R R −∞ form R r ˜ (7) F(r) = w + f(s)ds Z0 (for r > 0). Note that the presence of the δ-function in the pdf translates to a nonzero cdf at the zero boundary (more precisely, lim F(r) = 0). r→0+ 6 2.2. Option-implied risk-neutral pdfs. Consider a eurodollar futures option with strike K and maturity τ.11 Let r˜ denote the futures rate at the maturity of the option, 9See, for example, Dirac (1958,p58). One place the δ-function appears in the finance literature is in thespecificationofinitialconditionsfortheKolmogorov(Fokker-Planck)equations;see,e.g.,Goldstein and Kierstead (1997). 10Note, incidentally, that the pdf of the absorbing boundary models also contains a delta function. 11In the case of eurodollar futures options, the futures price is quoted as 100 [futures rate]. − Therefore,aputoptiononfutures contract (futures price)isacalloptiononfutures rate. Inthispaper, “put option” will always mean put option on the futures rate, and “strike” is in terms of futures rate (not futures price).

6 D. H.KIM and let f (r˜) be the risk-neutral pdf of r˜. The put option price P (K) and the call τ τ option price C (K) can be expressed in terms of f as τ τ K (8) P (K) = e−rfτ (K r˜)f (r˜)dr˜, τ τ − Z−∞ ∞ C (K) = e−rfτ (r˜ K)f (r˜)dr˜, τ τ − ZK where the prefactor e−rfτ denotes the price of a zero-coupon bond maturing τ periods hence.12 In the case of the eurodollar futures options, the futures rate at the option maturity is the 3-month LIBOR. For simplicity, we shall treat this rate as the short rate of interest, and drop the “tilde” symbol from r˜. Taking successive derivatives of both sides of eq. (8) gives a simple formula for the pdf, (9) f(K) = erfτP′′(K) = erfτC′′(K), where the double prime notation denotes second derivatives; the maturity index τ has been dropped for notational simplicity. In practice, options are available at discrete values of K (say, Kp,Kp,... for put options and Kc,Kc,... for call options), and ob- 1 2 1 2 served prices could contain some errors. Taking numerical second derivatives (with an interpolated function) to obtain f(r) typically leads to unrobust results.13 Therefore, a common practice is to assume a flexible functional form for f(r;γ), where γ collectively denotes the parameters of f, and solve the least squares problem (10) γˆ = argminζ(γ), np nc ζ(γ) = (Pob(Kp) P(Kp))2 + (Cob(Kc) C(Kc))2 +(Fob F)2 i − i i − i − i=1 i=1 X X where the superscript ob denotes the observed values. Including the futures pricing error helps to pin down the mean of the distribution; note that the model-implied current ∞ futures value F is E(r) = rf(r;γ)dr. This least squares problem amounts to min- −∞ imizing the pricing errors in absolute terms. Because options far out of money tend R to have very low prices, this metric in effect substantially discounts the information in such options, as compared to near-the-money options. One could also consider minimizing the pricing errors in relative terms (i.e., replacing (Pob P)2 in eq. (10) by − (P Pob)2/Pob2, etc.), but this is not advisable, as the farther out-of-money options − 12More precisely, we have P (K) = E[exp( τ r ds)max[0,K r˜]], where E denotes τ − 0 s −τ the risk-neutral expectation. Making the approximation E[exp( r ds)max[0,K r˜]] τ R − 0 s − ≈ E[exp( r ds)]E[max[0,K r˜]] (which is common in this literature) gives eq. (8). − 0 s − R 13To overcome this problem, Shimko (1993) proposed to fit a smooth functional form for the Black- R Scholes implied volatility and then take derivatives of the Black-Scholesformula, but this method does not tell much about the tail parts of the pdf.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 7 prices are expected to be less reliable, due to their lower liquidity and more pronounced market microstructural effects (minimum tick sizes, etc.). There are many possibilities for parametrizing f(r;γ). The general wisdom is that a flexible parametric form with 4 to 6 parameters is sufficient to match observed option pricesatvariousstrikes.14 Aparametrizationsubstantially richerthanthiscanoftenlead to strange implications at minimal improvements in the fit of option prices (overfitting). Perhaps the most popular parametric form is the “mixture of lognormals”:15 1 (logr µ )2 1 (logr µ )2 1 2 (11) f(r) = β exp − +(1 β) exp − . √2πσ r − 2σ2 − √2πσ r − 2σ2 1 (cid:18) 1 (cid:19) 2 (cid:18) 2 (cid:19) Note that if β = 1, eq. (11) becomes just the lognormaldistribution, which is an implicit assumption behind the so-called Black implied-volatility formula (Black (1976)), which is often used to quote cap/floor and swaption prices. Parametric forms like f(r) in eq. (11) cannot be derived from known classes of term structure models. However, certain parametric forms can be derived from a specific class of term structure models. In Section 4.1, we shall discuss such an example, where an implicit parametric form for f(r) is derived from the multifactor CIR and quadratic-Gaussian models. It is also worth noting that despite the typically good performance, forms like (11) do have a limitation that is particularly relevant to our problem (pdf modeling in the proximity of the zero bound). Because the lognormal distribution has the feature that f(r) 0 as r 0, the form (11) can have difficulty when the actual distribution → → contains a δ-function at r = 0. As noted earlier (eq. (5)), the delta function δ(r) can be approximated by a lognormal distribution with mean and variance close to zero, which can be achieved by having a small σ and a large and negative µ in eq. (11), 1 1 but in such a two-component mixture model this may incur a substantial cost in the fit of other aspects of the pdf. Section 4.2 develops parametric forms for f(r) that can describe Black’s boundary behavior more easily and naturally. To have a manageable scope, this paper will be focusing only on the processes and distributionsintherisk-neutralmeasure(whichdeterminespricing),andnotthephysical (real world) measure.16 Therefore, henceforth we shall often drop the adjective “riskneutral” for brevity. 3. One-factor models 3.1. Black-Vasicek model. Let us now derive the pdf and option prices for the simplest Black’s boundary model, which has the shadow rate described by the one-factor 14See, e.g., Jackwerth (2004)’s survey and the papers in BIS (1999). 15McManus(1999,includedinBIS(1999))hascomparedwithperformanceofthelognormal-mixture pdf with other parametric forms in the context of eurodollar futures options. 16Nonetheless, the risk-neutraland physical processes take the same (or similar) form in many term structure model specifications, thus many of the qualitative findings in the present paper also apply to the physical measure as well.

8 D. H.KIM Vasicek process: (12) r = max[0,x ], t t (13) dx = κ(θ x )dt+σdB . t t t − This model, which we shall refer to as the “Black-Vasicek model,” has been theoretically and empirically studied by Gorovoi and Linetsky (2004) and Ueno, Baba, and Sakurai (2006) and others, but its risk-neutral pdf has not been discussed in the literature, to my knowledge.17 The conditional distribution of this model, i.e., the transition density f(r ), can t+τ t |F be easily obtained if one thinks in terms of the shadow rate process x rather than the r t t process. Note that if x > 0, we have r = x ; hence in this case, the distribution t+τ t+τ t+τ of r equals that of x . On the other hand, if x 0, we have r = 0. Thus all t+τ t+τ t+τ t+τ ≤ scenarios in which x is non-positive maps to a single point r = 0. Therefore, we t+τ t+τ have (14) f(r x ) = fS(x x )I(x )+w(x )δ(x ), t+τ t t+τ t t+τ t t+τ | | where fS(x x ) denotes the transition density of the x process, w(x ) is the weight t+τ t t t | of fS(x x ) in the nonpositive region, i.e., t+τ t | 0 (15) w(x ) = fS(x x )dx , t t+τ t t+τ | Z−∞ and I(x) is the indicator function (I(x) = 1 for x 0 and I(x) = 0 for x < 0). It ≥ may be of interest to note that the distribution in eq. (14) can be thought of as a special case of the “mixture” distribution, consisting of a truncated normal distribution 1 fS(x x )I(x ) and an infinitely sharp distribution δ(x ) with weights 1 1−w(xt) t+τ | t t+τ t+τ − w(x ) and w(x ), respectively. t t The transition density function fS(x x ) for the Vasicek process is well-known: t+τ t | Since the stochastic differential equation (13) has the solution τ (16) x = e−κτx +(1 e−κτ)θ+ e−κ(τ−s)σdB t+τ t s − Z0 the random variable x = x is conditionally normally distributed18 t+τ (17) x N(µ˜,σ˜2), ∼ with conditional mean and variance given by (18) µ˜ = e−κτx +(1 e−κτ)θ, t − σ2 σ˜2 = (1 e−2κτ). 2κ − 17GorovoiandLinetsky (2004)derive anexpressionfor the optionpricesin the Black-Vasicekmodel using eigenfunction expansions, but they do not discuss pdfs. 18In this paper, the symbol means “is distributed as”. ∼

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 9 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 −2 0 2 4 6 fdp (a) 1 0.8 0.6 0.4 0.2 0 −2 0 2 4 6 fdc (b) Figure 1: Black-Vasicek model pdf (a) and cdf (b). The verticalline at 0 in Figure (a) pictures a delta function. Thus, suppressing time indices and simplifying notations, we have the conditional distribution of the short rate given by 1 r µ˜ (19) f(r) = wδ(r)+ φ − I(r), σ˜ σ˜ (cid:18) (cid:19) wherew = Φ( µ˜/σ˜),andφ( )andΦ( )arethestandardnormalpdfandcdf, respectively − · · (i.e., φ(x) = e−x2/2/√2π, Φ(x) = x dte−t2/2/√2π). Integrating eq. (19) gives the −∞ cdf,19 R (20) F(r) = [w+Φ((r µ˜)/σ˜) Φ( µ˜/σ˜)]I(r) = Φ((r µ˜)/σ˜)I(r). − − − − The pdf and cdf of the short rate in the Black-Vasicek model are illustrated in Figure 1a,b. The dashed line in Figure 1b is the part of the shadow rate distribution that collapses onto the delta function at r = 0, denoted by a vertical line at zero in Figure (a). Note also that the cdf F(r) approaches a finite number as r 0+ (Figure (b)). → The prices of τ-maturity put/call options with strike K are straightforward to evaluate. Substituting eq. (19) into eq. (8), we have (21) P(K) = e−rfτ[wK +(K µ˜)(Φ(k) Φ( µ˜/σ˜))+σ˜(φ(k) φ( µ˜/σ˜))], − − − − − (22) C(K) = e−rfτ[σ˜φ(k) (K µ˜)(1 Φ(k))], − − − where k (K µ˜)/σ˜. ≡ − 19Note that the indicator function I(r) and the Dirac delta function δ(r) are related by the identity dI(r) =δ(r). Thus, taking the derivative of eq. (20), we recover the expression (19). dr

10 D. H.KIM 3.2. Two zero-cdf models. Let us nowconsider two well-known andtractable positive interest rate models: the CIR model and the QG model. The CIR model is given by (23) dr = κ(θ r )dt+σ√r dB . t t t t − It is well known (CIR, 1985) that the random variable r r in this model is condit+τ ≡ tionally distributed (at time t) as (24) r bχ2 , ∼ ν,λ where χ2 denotesthenoncentral chi-squared randomvariablewithν degrees offreedom ν,λ and noncentrality parameter λ, and σ2(1 e−κτ) (25) b = − 4κ 4κθ ν = σ2 4κ λ = e−κτr . σ2(1 e−κτ) t − Consider now the positive-definite 1-factor QG model (26) r = x2, t t dx = κ(θ x )dt+σdB . t t t − Since the random variable x = x is conditionally normally distributed (eq. (17)), it t+τ can be seen immediately that the distribution of random variable r = r is given by t+τ eq. (24), i.e., the same form as the CIR model, with (27) b = σ˜2 ν = 1 λ = µ˜2/σ˜2, where µ˜ and σ˜ were given in eq. (18). The pdf of the r bχ2 distribution (CIR/QG models) is given by ∼ ν,λ 1 r (28) f(r) = g I(r), b b (cid:16) (cid:17) with 1 x (ν−2)/4 (29) g(x) = e−(λ+x)/2 I (√λx), (ν−2)/2 2 λ (cid:16) (cid:17) whereI ( )isthemodifiedBessel functionofthefirstkind ofordera.20 Theoptionprices a · and cdf are not available in simple form for this distribution, but it is straightforward to compute them by one-dimensional numerical integrations with f(r). 20See, for example, Johnson, Kotz, and Balakrishnan (1994).

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 11 3.3. Comparing the model implications. To get a feel for the qualitative difference between Black’s boundary models and zero-cdf (CIR/QG) models, it is instructive compare the behavior of the pdf f(r) in the r 0 limit. → For the CIR/QG model, expanding the expression (29) in powers of r, using the well-known series expansion of the modified Bessel function of the first kind 1 ∞ (z2/4)j (30) I (z) = ( z)a , a 2 Γ(j +1)Γ(a+j +1) j=0 X gives e−λ/2 e−λ/2(λ ν) (31) f(r) = r ν 2 −1 + − r ν 2 +O(r2 ν+1), (2b) ν 2Γ(ν) 2(2b) ν 2 +1Γ(ν +1) 2 2 where the notation O(rn) denotes terms of order n or higher.21 Thus, the cdf is given by e−λ/2 e−λ/2(λ ν) (32) F(r) = r ν 2 + − r ν 2 +1 +O(r ν 2 +2). (2b) ν 2 νΓ(ν) 2(2b) ν 2 +1(ν +1)Γ(ν +1) 2 2 2 2 Fromeq. (31) it canbeseen that depending onν, we canexpect qualitatively different behavior of f(r): , ν < 2 (33) limf(r) = ∞ 0, ν > 2. r→0 (cid:26) Recall that, for the CIR model (eq. (23)), ν is given by eq. (25), therefore, the condition ν > 2isthesameasFellercondition(κθ > 1σ2). Itmakessensethataccessible boundary 2 case (Feller condition violated) shows more probability density in the low r region than the inaccessible boundary case. Note also that since ν = 1 for the QG model, we have lim f(r) = for the QG r→0 ∞ model. As the zero boundary is accessible in the QG model (since x in eq. (26) is t unrestricted, it can reach zero), it shows a qualitatively similar small-r behavior as the CIRmodelwithviolatedFellercondition. Itshouldbenoted, however, thatanaccessible zero boundary does not always imply lim f(r) = , as we shall see with the case r→0 ∞ of the 3-factor QG model (Sec. 4.1). Note also that even in the case lim f(r) = r→0 ∞ (one-factor CIR model with violated Feller condition and the QG model), the cdf is zero at r = 0, since F(r) rν/2 in the r 0 limit. ∝ → The small-r behavior of the pdf in eq. (31) implies that the put option price for small strike K is given by K e−λ/2 4 (34) P(K) = e−rfτ (K r)f(r)dr = e−rfτ K ν 2 +1 +O(K ν 2 +1). Z0 − (2b) ν 2Γ(ν 2 )ν(ν +2) 21This result can be also derived using the Laplace transform technique, as shown in Appendix A.

12 D. H.KIM 0.5 0.4 0.3 0.2 0.1 0 0 2 4 6 8 fdp (a) 0.5 0.4 0.3 0.2 0.1 0 0 2 4 6 8 fdp (b) −3 x 10 2 1.5 1 0.5 0 0 0.05 0.1 0.15 0.2 strike ecirp tup (c) 0.2 0.15 0.1 0.05 0 0 0.5 1 strike ecirp tup (d) Figure 2: (a) The Black-Vasicek pdf (thick solid line) and inacessible boundary CIR/QG model pdf (thin dashed line). (b) Accessible boundary model pdf (thin solid line). (c) Put option prices for the Black-Vasicek model (thick solid line), accessible boundary CIR/QG model (thin solid line), and inaccessible boundary model (thin dashed line). (d) Same as (c). These asymptotic behaviors of CIR/QG models can be compared with the those of the Black-Vasicek model: φ( µ˜/σ˜) (35) f(r) = wδ(r)+ − +O(r1) σ˜ φ( µ˜/σ˜) F(r) = w + − r +O(r2) σ˜ P(K) = e−rfτwK +O(K2). (for r,K > 0), which are easily derived from the formulae in Section 3.1. Figure 2 graphically illustrates the qualitative difference between Black’s boundary models and CIR/QG models. Figure 2a plots, in thick solid line, the pdf of the Black- Vasicek model with parameters in eq. (19) given by µ˜ = 2,σ˜ = 1.2, which imply w = 0.0478. It also shows, in thin dashed line, the pdf of the r bχ2 model with ν = 3 ∼ ν,λ (corresponding to the CIR model satisfying the Feller condition), whose b and λ were

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 13 determined such that the model matches the mean and variance of the Black-Vasicek model;22 requiring the means and variances of the distributions be similar amounts to having comparable at-the-money option prices. Figure 2b shows the pdf of the r bχ2 ∼ ν,λ model with ν = 1 (corresponding to the CIR model violating the Feller condition and the QG model) in thin solid line. As discussed with eq. (31), this pdf has a singularity of the form r−1/2, but it has a local minimum at a point very close to zero,23 and above that it looks similar to the pdf of the r bχ2 distribution with ν = 3 (i.e., Figure 2a). ∼ ν,λ Despite the rough similarity in the shape of the pdf of the Black-Vasicek model and the r bχ2 model, the model-implied put prices are drastically different, with the ∼ ν=1,λ Black-Vasicek model producing far larger numbers, as shown in Figure 2c. This is because the r−1/2 singularity in the r bχ2 model is a much weaker singularity than ∼ ν=1,λ the δ-function singularity in the Black-Vasicek model; to put it simply, Black’s boundary model has a much larger probability weight at or near zero than the CIR/QG model. Note also that in the low interest rate region the Black-Vasicek model shows a linear dependence on K, as predicted by eq. (35). Although the r,K 0 behavior of the pdfs and put option prices derived above for → the one-factor Black-Vasicek and CIR/QG models are useful for getting a firm handle on the zero boundary behavior of these models, these predictions might not be cleanly testable empirically, as options might not exist (trade) at very low strikes, and even if they existed their quality may be in doubt. Furthermore, certain market realities (to be discussed at the end of Sec. 5.2) may blur the clean asymptotic behaviors. Still, even if one moves onto regions where the asymptotic behaviors like (34) and (35) are no longer accurate, the basic intuition would still carry over, and one would expect Black’s boundary model to produce higher put option prices for small K’s than zero-cdf models (with comparable mean andvariance) when the rates are low enough (or the distribution is wide enough) that the weight of the δ-function piece in Black’s boundary model is non-negligible. This is illustrated in Figure 2d (which plots the same objects as Figure 2c but on a larger scale). It is also interesting to note that the put option prices for the r bχ2 model (inaccessible boundary) and the r bχ2 model (accessible ∼ ν=3,λ ∼ ν=1,λ boundary) are quite similar, beyond the very small K region. 4. Richer parametric forms The one-factor models considered in the previous section, though instructive, may be too parsimonious to capture option prices accurately for a broad range of strikes. Let us now therefore consider richer versions of zero-cdf models and Black’s boundary models. 22Since the mean and variance of the r bχ2 distribution are given by b(ν+λ) and 2b2(ν +2λ), ∼ ν,λ respectively, for a given ν (=1,3 here) one can easily solve for b and λ. 23The approximate location of this point can be determined by setting the derivative of eq. (31) to zero and solving for r.

14 D. H.KIM 4.1. Multifactor zero-cdf models. A non-negative, multifactor version of the CIR model can be constructed by adding up independent CIR factors: (36) r = x +x + +x , t 1t 2t nt ··· dx = κ (θ x )dt+σ √x dB it i i it i it it − where theBrownianmotionsB , B ,.. areindependent ofeach other. This “multifactor 1t 2t CIR” model has been used in many settings, e.g., Duffie and Singleton (1997) and Feldhu˝tter and Lando (2007).24 Because x ’s are independent, the results about the 1-factor CIR model carries over i easily; it is straightforward to show that the distribution of r = r is given by t+τ (37) r b χ2 +b χ2 + +b χ2 , (b > 0), ∼ 1 ν1,λ1 2 ν2,λ2 ··· n νn,λn i where χ2 ’s are independent noncentral chi-squared random variables with ν degrees νi,λi i of freedom and noncentrality parameter λ , and i σ2(1 e−κiτ) (38) b = i − i 4κ i 4κ θ i i ν = i σ2 i 4κ e−κiτ i λ = x . i σ2(1 e−κiτ) it i − Let us now consider the general n-factor positive-definite QG model, (39) r = (x α)TΨ(x α), t t t − − dx = (θ x )dt+ΣdB , t t t K − where the superscript T denotes matrix transpose; x is an n-dimensional vector x = t t [x ,...,x ]T; and Σ are n n constant matrices; α and θ are n-dimensional constant 1t nt K × vectors, and B is an n-dimensional vector of standard Brownian motions. In order to t insure the non-negativity of r , we require Ψ to be a symmetric positive-definite matrix. t It is straightforward to show that the conditional distribution of the short rate in the positive-definite QG model also takes the form (37),25 with (40) b = d (> 0) i i ν = 1 i λ = ([ TC−1(µ˜ α)] )2, i i P − 24The multifactor CIR model belongs to Dai and Singleton (2000)’s A (n) affine model classificaj=n tion. Other kinds of multifactor affine models (the A (n) models) are not considered in this paper, j<n as they do not respect the non-negativity condition for nominal interest rates. The concluding section of this paper, however,does make some comments relevant to these models. 25For a derivation, see Ahn, Dittmar, and Gallant (2002, Appendix B).

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 15 where C, µ˜ are given by (41) µ˜ = e−Kτx +(I e−Kτ)θ, t − (42) C = Chol(Ω˜), (i.e., Ω˜ = CCT), τ (43) Ω˜ = e−K(τ−s)ΣΣTe−KT(τ−s)ds, Z0 (the notation eX denoting the matrix exponential), and d and are from the diagonali P ization of the symmetric positive-definite matrix CTΨC: (44) CTΨC = D −1 = D T, P P P P where D is a diagonal matrix with diagonal elements d ,d ,...,d (the eigenvalues of 1 2 n the CTΨC matrix). Note that the dependence of λ on x has been suppressed in the i t notation. Models like the multifactor CIR model have been criticized for the restrictive nature of the factor correlation structure. However, the positive-definite QG models have been widely believed to accommodate a rich factor correlation structure, due to the ability to specify and Σflexibly; see, e.g., Ahn, Dittmar, and Gallant (2002). Therefore, it bears K emphasizing that not only the multifactor CIR model but also the positive-definite QG modelhas a short ratedistribution given by a positive linear combination of independent noncentral χ2 random variables, and that this holds no matter how flexible the model is and no matter how many factors it has. The pdf f(r) and cdf F(r) of the distribution for the multifactor CIR/QG model (eq. (37)) do not have simple closed-form expressions in general, but the characteristic function f ˆ (w) E(eiωr) has a simple expression ≡ n λ b iω (45) f ˆ (w) = (1 2b iω)−νj/2exp j j , j − 1 2b iω j=1 (cid:18) − j (cid:19) Y due to the independence of the noncentral χ2 variables in eq. (37). From this, f(r) and F(r) can be computed by evaluating the well-known Fourier-inversion integrals (see, e.g., Stuart and Ord (1994, Chap. 4)), 1 ∞ (46) f(r) = Re(e−iωrf ˆ (ω))dω π Z0 1 1 ∞ Im(eiωrf ˆ ( ω)) (47) F(r) = + − dω. 2 π ω Z0 The option prices P(K) and C(K) for this model (integrals (8)) are also not available in closed-form, but they can be expressed in terms of certain cumulative distribution functions. For the computation of pdfs from options data (eq. (10)), it is desirable to have a formula that allows fast computation of option prices. Appendix B derives such a formula, using the the celebrated saddle-point formula of Lugannani and Rice (1980).

16 D. H.KIM Let us now consider the small-r behavior of f(r). Appendix A shows that (48) f(r) = e−Σiλi/2 1 r 1 2 Σiνi−1 + Σ i (λ i − ν i )/4b i r2 1Σiνi + . Π i (2b i )νi/2 (cid:18) Γ( 2 1Σ i ν i ) Γ(1 2 Σ i ν i +1) ··· (cid:19) This form has an interesting dependence on the number of factors: the exponent in the leading term increases with number of factors. For the QG model, we have , n = 1, ∞ (49) f(r 0) = finite, n = 2, →  0, n = 3.  This pattern is due to the fact that when there are more factors it is less likely for r  to approach 0 since all of the factors have to become small. Therefore, higher-factor QG models are even more different from Black’s boundary models insofar as the zeroboundary behaviors are concerned. To get further insights into the distribution described by eq. (37), it is useful to examine its cumulants (hence moments), which is straightforward to calculate from the ˆ cumulant generating function (log(f( iω))). The mean, variance, and the third central − moment of this distribution are easily derived from the cumulants of the noncentral χ2 variables (see, e.g., Johnson, Kotz, Balakrishnan (1994, p447)): (50) E(r) = Σn b (ν +λ ) i=1 i i i Var(r) = Σn 2b2(ν +2λ ) i=1 i i i E[(r E(r))3] = Σn 8b3(ν +3λ ). − i=1 i i i It can be seen that because b ,ν ,λ > 0 for all i, the skewness i i i E[(r E(r))3]/Var(r)3/2 is always positive. − More generally speaking, one cannot create a negatively skewed distribution by forming a positive linear combination of independent, positive-skewed random variables. Therefore, a model like (51) r = b x +b x + +b x , (b > 0), t 1 1t 2 2t n nt i ··· where x is Ahn and Gao (1999)’s one-factor process for the short rate, would again it have the feature that the short rate distribution is always positively skewed.26 4.2. Flexible pdfs with Black’s boundary behavior. This section develops flexible parametric forms for the risk-neutral pdf that are consistent with Black’s boundary behavior. The earlier discussion with the Black-Vasicek model points to the way: write 26I have not seriously attempted to prove that the conditional distribution of the short rate in Ahn and Gao (1999)’smodel [r =1/x , dx =κ˜(θ˜ x )dt+σ˜√x dB , (κ˜,θ˜>0, 2κ˜θ˜>σ˜2)] is positively t t t t t t − skewed for all feasible values of the parameters, but I believe so, as the right tail of the distribution is thicker and longer than the left tail. (It can be shown that f(r) r −ae −b/r (a > 7/4, b > 0) for r 0, and f(r) r −a (a>2) for r .) ∝ → ∝ →∞

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 17 down a flexible form for the shadow rate distribution fS that has some weight below zero, and take r (52) f(r) = fS(r)I(r r)+wδ(r r), w = dxfS(x). − − Z−∞ Here we have a slightly more general expression in which the lower boundary r is not necessarily zero. The distribution fS can be modeled with similar degree of parsimony as the pdfs that are in use in the extant literature. A natural candidate for fS is the “mixture of normals” form, i.e., 1 (r µ )2 1 (r µ )2 (53) fS(r) = β exp − 1 +β exp − 2 , 1 √2πσ − 2σ2 2 √2πσ − 2σ2 1 (cid:18) 1 (cid:19) 2 (cid:18) 2 (cid:19) where β = 1 β . Another possibility is the Gram-Charlier/Edgeworth expansion 2 1 − 1 r µ γ r µ γ r µ (54) fS(r) = φ − 1+ 3 H − + 4 H − + , 3 4 σ σ 6 σ 24 σ ··· (cid:18) (cid:19)(cid:18) (cid:18) (cid:19) (cid:18) (cid:19) (cid:19) where H are Hermite polynomials (H (x) = x3 3x, H (x) = x4 6x2+3, etc.). This j 3 4 − − expansion differs from the Gram-Charlier/Edgeworth expansion commonly used in the option-implied pdf literature by the fact that it uses unit normal variable, rather than unit lognormal variable. Note that in both eq. (53) and (54), fS is allowed to have some probability weight below r = 0. These forms, which are richer than that of the Black-Vasicek model (eq. (19)), can be viewed as accommodating a more complicated process for the shadow rate, e.g., a multifactor model with stochastic volatility. The formula for P(K), C(K) and F for the mixture-of-normals shadow rate model (eq. (53)) and the Gram-Charlier/Edgeworth shadow rate model (eq. (54)) are provided in Appendix C. 5. Empirical evidence 5.1. Empirical strategy. To examine how well zero-cdf models perform in matching option prices, I have implemented two versions of pdf models, which we shall refer to as the QG3 model and the CIR2 model: (55) QG3 : r b χ2 +b χ2 +b χ2 , ∼ 1 1,λ1 2 1,λ2 3 1,λ3 (56) CIR2 : r b χ2 +b χ2 , ∼ 1 ν1,λ1 2 ν2,λ2 (b ,ν ,λ > 0). For Black’s boundary model, I have implemented the model in which the i i i shadow rate distribution is that of the mixture of normals, i.e., eq. (53), which we shall refer to as the B-MN model.27 For comparison with a popular benchmark pdf, I have 27I have also explored empirically the Black’s boundary model with Gram-Charlier/Edgeworth shadow rate distribution (eq. (54)), and obtained results that are similar to the B-MN model. For brevity, only the B-MN model results will be reported in this paper.

18 D. H.KIM model pdf parameters description QG3 b ,b ,b ,λ ,λ ,λ 3-factor QG(CIR) models. 1 2 3 1 2 3 CIR2 b ,b ,ν ,ν ,λ ,λ 2-factor CIR(QG) models. 1 2 1 2 1 2 B-MN β ,µ ,µ ,σ ,σ Black’s boundary; mixture of normals. 1 1 2 1 2 MLN β ,µ ,µ ,σ ,σ Mixture of lognormals. 1 1 2 1 2 Table 1: Summary of the models. The pdf parameters for the QG3 and CIR2 models are defined in eqs. (55) and (56). The pdf parameters for the MLN and B-MN models are defined in eqs. (11) and (53), respectively. also computed the pdf based on the mixture of lognormals, i.e., eq. (11), which we shall denote as the MLN model. These four models are summarized in Table 1. The QG3 model covers all possible 3-factor positive-definite QG models, but it also includes a special case of the 3-factor CIR model (κ θ = 1σ2). The CIR2 model covers i i 4 i all possible cases of the 2-factor CIR models, but it also nests the 2-factor QG model. This model also nests the Longstaff-Schwartz (1992) model, which was regarded by at least H˝ordahl (2000) as a promising model for describing the short rate distribution. In terms of the number of parameters, the QG3 and CIR2 models are richer than the B-MN and MLN models by one parameter. It is also worth noting that, because the QG3 and CIR2 models originate from specific term structure models, the pdf parameters b ,ν ,λ (ν = 1 for QG3) are linked to i i i i the elementary parameters of the term structure model and the state variables x .28 it Therefore the pdf parameters can be obtained from an estimated term structure model. However, there are many ways of estimating term structure models,29 which could lead to different estimated elementary parameters and state variables, and in turn, different pdfs. Therefore, in this paper we shall simply treat the pdf parameters for the QG3 model and the CIR2 model as free parameters to be determined in the optimization problem (10), in the same way as the pdf parameters for the B-MN and MLN models are determined. This way can be viewed as giving the QG3 and CIR2 model a maximum chance to fit options data.30 Lastly, note that the estimation of multifactor CIR and QG models in the literature have often included a constant term in the short rate, i.e., r = c+x + +x and t 1t nt ··· 28Recall, from eqs. (38) and (40), that λ in the multifactor CIR/QG models depend on x . i it 29For example, one could use just the current yield curve information (daily recalibration) or use historical yields data to estimate the model in a time-consistent way. The estimation technique itself couldbe ofmany variety(e.g., MLE,GMM, etc.). One couldalso optoropt notto use options data in the estimation. 30This is analogous to the daily recalibration of term structure models that is often done in finance industry. Note that after setting the left hand side of eq. (38) or eq. (40) to the pdf parameters determined this way,one can searchfor the elementary parametersandstate variablesthat best fit the yield curve, but the yield curve fit may be poorer than the unconstrained case.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 19 r = c+(x α)TΨ(x α). Empirically, Duffie and Singleton (1997) obtained a negative t t t − − c for the multifactor CIR model, but in that case it is no longer a non-negative interest rate model. Ahn, Dittmar, and Gallant (2002) obtained a positive c in their estimation of the QG model; in this case, there is no probability weight in the region [0,c]. In any case, the inclusion of the constant term does not change major qualitative behaviors of these models, including the prediction that the pdfs are always positively skewed. 5.2. Data. For the empirical investigation, we shall confine attention to the eurodollar futures option with maturity of about 1 year. The 1-year maturity was chosen so as to have a horizon that is reasonably short but at the same time long enough to have a sizeable width of the distribution. To get a sense of the performance of the models in a variety of contingencies, I have computed the pdfs for the above four models for every last business day of the month from May 1998 to March 2008. Besides these monthly data, I also have computed the pdfs from a second dataset consisting of daily data from Feb 4, 2008 to April 28, 2008. As we shall see, this period is particularly interesting as regards the zero boundary behavior question. The options data used in this paper are from CME (Chicago Mercantile Exchange). Theprefactore−rfτ ineq. (8)wascomputed using the1-yearLIBOR.TheCMEeurodollar futures options are American options. I have converted their price into corresponding European price by subtracting the early excercise premium computed based on Baroni- Adesi and Whaley (1987)’s formula; these corrections tended to be fairly small. The eurodollar futures options markets have fixed maturity dates (identical to maturity dates for the eurodollar futures contracts, 4 days a year); for each date in the sample, I chose the option maturity that is closest to 365.25 days. This creates a sawtoothbehavior of the time series of option maturities (oscillating between approximately 7 and 9 year) for the monthly sample. 8 8 The options are typically available at strikes with multiples of 25 basis points, but sometimes there are options with strikes at odd multiples of 12.5 basis points, which are dropped as they tend to be not as liquid as the even multiples. Also, in order to avoid illiquid options, only the options with strikes at which there is a nonzero open interest are used. Thus the number of call and put options used in the pdf computations varies over time. Figure 3 shows the minimum put strike and the maximum call strike that were available on each day of the monthly sample. Also shown is the put strike rate that corresponds to the current futures rate minus 1.5σ∗, and the call strike rate that corresponds to the futures rate plus 1.5σ∗, where the proxy σ∗ for the standard deviation of the pdf is computed as (57) σ∗ = √τ σBF,

20 D. H.KIM 10 9 8 7 6 5 4 3 2 1 0 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 Figure 3: Available option strikes for 1-year maturity options. Maximum call strike and minimum put strike are shown in thick solid lines. The current futures rate F is shown in thin solid lines. The thin ∗ dashed lines are F 1.5σ . ± where σB is the Black implied volatility (Black (1976)) and τ is the option maturity.31 As can be seen, these numbers (thin dashed line) tend to lie inside the [min(put strike), max(call strike)] band, indicating that typically a fairly wide range of strikes are available, spanning at least 1.5σ∗ to 1.5σ∗ units around the mean of the distribution. I − have used only the out-of-money puts and calls (puts at or below the current futures rate and calls above the current futures rate). Thus at each multiples of 25bp (between the minimum put strike and maximum call strike), there is only one option (either put or call), or none (when there is no open interest). Table 2 provides some statistics regarding the liquidity of options proxied by open interest, for strikes 0.25σ∗,0.5σ∗,σ∗,1.5σ∗ away from the mean.32 Options about 1 2 standard deviation from the current futures rate are the most active, but the whole strike range (at least up to 11 standard deviation) is fairly active; e.g., the median open 2 interest in the 1.5σ∗ put option is more than a third of that of the 0.5σ∗ option. Lastly, as a caveat, note that since the r here is taken to be the 3-month LIBOR, the “zero bound” (more precisely, the lower bound) is not exactly at zero, as the 3-month LIBOR contains a premium for credit/liquidity risk, besides some differential between 31Though simple, this number agrees fairly well with the standard deviation computed from the B-MN model. 32Since options are available at discrete strikes in multiples of 0.25%, the strikes closest to the respective numbers were chosen. If there is no strike within 25 basis points of the number, it is not ∗ counted. This will be the convention also when the pricing errors for strikes 1.5σ away from the ± mean are examined in Sec. 5.3.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 21 put call dist. from F mean 25% 50% 75% mean 25% 50% 75% 0.25σ∗ 48664 16320 32749 58906 28072 10106 18352 36747 0.5σ∗ 49449 15436 37047 57529 32491 9994 20434 37852 σ∗ 45561 4625 25362 99226 27141 9868 19236 34176 1.5σ∗ 31512 4738 13262 37920 19890 3726 10412 26950 Table 2: Mean open interest, and 25-,50-,and 75-percentile open interest. the 3-month rate and the instantaneous rate. This means that the proper location of the delta function for the Black’s boundary model pdf (r) should be a small but nonzero number. Furthermore, the delta function itself would be smudged out a bit, in other words, 1 (x r)2 ˜ (58) δ(x r) = exp − − √2πǫ − 2ǫ2 (cid:18) (cid:19) where ǫ is a small number (but not ǫ 0 as in eq. (5)), would be a more accurate → description. However, as long as ǫ is small enough, the delta function δ(x r) should − ˜ provide a good approximation of δ(x r). − It is difficult to get a precise sense of r for the 3-month LIBOR, but the recent Japanese experience can provide some guidance: even at the height of Bank of Japan’s Quantitative Easing Program, the 3-month yen TIBOR rate (an analogue of the US dollar LIBOR) did not go below 8 basis points. For simplicity, the main computations in this paper assumed that the lower bound r is at 0, but I have also computed the pdf for the B-MN model with r = 0.15% to check the robustness of the results. Intuitively one would expect that the difference between these two settings would be immaterial for option strikes that are substantially larger than the effective lower bound (i.e., K r), ≫ but at a very low strike (e.g., 25 basis points), there may be a more visible difference between r = 0 vs r = 0.15%. 5.3. Option pricing performance. To get a sense of the overall fit across strikes, it is useful to look at the root-mean-square errors, i.e., (59) RMSE = (ζ ˆ /(1+n +n ))1/2, p c ˆ where ζ is the minimized value of ζ in eq. (10). Figure 4 shows the time series of RMSEs based on the four models (B-MN, MLN, CIR2, QG3) for the 1998-2008 period. As can be seen in the top panel, the B-MN model (Black’s boundary model) and the MLN model (conventional benchmark) performed similarly well, the MLN model producing a somewhat smaller RMSE in 2003-2004 and largerRMSEin2007-2008. Notethatduringmostofthedataperiod,withtheexceptions around 2003 and 2008, the distributions were sufficiently away from the zero bound (as

22 D. H.KIM 0.01 B−MN 0.008 MLN 0.006 0.004 0.002 0 98 99 00 01 02 03 04 05 06 07 08 09 0.04 CIR2 QG3 0.03 0.02 0.01 0 98 99 00 01 02 03 04 05 06 07 08 09 Figure 4: Rootmean square optionpricing errors. Top panel: B-MN andMLN models. Bottompanel: QG3 and CIR2 models. can be seen in Figure 3), thus the B-MN model was simply the “mixture of normals” model (i.e., f = fS); the results indicate that there is little practical difference between the normal-mixture model and the lognormal-mixture model most of the times. On the other hand, the QG3 model and CIR2 model led to notably larger RMSEs. (Note the difference in the y-axis scales in top panel and bottom panel). Between the two models, the RMSE numbers were very similar, practically lying on top of each other. It is also useful to examine the individual (strike-specific) pricing errors, Pmodel Pob − and Cmodel Cob. Table 3 displays summary statistics for the pricing errors for put − options with strikes 0.25σ∗ and 1.5σ∗ units below F and call options with strikes 0.25σ∗ and 1.5σ∗ units above F. It can be seen that the errors for the QG3 and CIR2 models are not only large, but they also deviate systematically from zero. For example, even the 75-percentile value of the pricing error for the put option at the strike of F 1.5σ∗ is − negative; thisnegative biasmeans thattheQG3andCIR2models tendtosystematically underprice the farther-out-of-money put options. On the other hand, they tend to overprice the near-the-money call options.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 23 0.25σ∗ 1.5σ∗ mean 25% 50% 75% mean 25% 50% 75% put B-MN -0.0014 -0.0030 -0.0010 0.0009 -0.0003 -0.0019 -0.0001 0.0012 MLN -0.0008 -0.0026 -0.0007 0.0017 -0.0017 -0.0032 -0.0009 0.0013 CIR2 0.0031 -0.0007 0.0031 0.0076 -0.0155 -0.0198 -0.0140 -0.0074 QG3 0.0028 -0.0011 0.0027 0.0077 -0.0157 -0.0198 -0.0141 -0.0083 call B-MN -0.0014 -0.0036 -0.0011 0.0006 0.0011 -0.0008 0.0010 0.0028 MLN -0.0011 -0.0027 -0.0010 0.0007 0.0009 -0.0010 0.0008 0.0027 CIR2 0.0043 0.0009 0.0044 0.0080 0.0005 -0.0056 0.0011 0.0056 QG3 0.0041 0.0008 0.0041 0.0079 0.0011 -0.0050 0.0010 0.0058 Table 3: Statistics (mean, 25-, 50-, and 75-percentile values) about put option pricing errors at strikes ∗ ∗ ∗ F 0.25σ andF 1.5σ (tophalfofthe table),andcalloptionpricingerrorsatstrikesF+0.25σ and − ∗ − F+1.5σ (bottom half of the table). In sum, the pdfs based on zero-cdf term structure models (CIR2 and QG3 models) performnotablyworsethantheMLNandB-MNmodelsdespite having onemoreparameter, suggesting that the kind of restrictions on the shape of the risk-neutral distribution implied by these models are not well supported by data. 5.4. Is the pdf always positively skewed? As pointed outin Section 4.1, well-known zero-cdf term structure models, including the multifactor CIR and positive-definite QG models, imply that the skewness of the risk-neutral distribution of the short rate is always positive. To check this prediction, we can examine the model-implied skewness E[(r E(r))3]/E[(r E(r))2]3/2 from the B-MN and MLN pdfs, shown in Figure 5a. − − The sharp vertical lines that occur in the MLN measure are a small number of occasions in which the MLN model produced very large skewness values that went out of bounds of the figure, even though option pricing errors were small. Ignoring these pathological cases, the B-MN and MLN models tell fairly similar story. In particular, the risk neutral distribution was negative during an extended period of time in recent years (2006-07). Negatively-skewed pdfs occurred also around 1998 and 2001. A test of skewness that does not need an assumption about the functional form of the pdf can be constructed as follows. Note that a negatively skewed distribution has mean < median, while a positively skewed distribution has mean > median.33 Since the cumulative distribution function F(r) at the median is 1/2 by definition, one can check whether the skewness is always positive by examining whether the condition (60) F(F) > 1/2 33This is the idea behind the Pearson measure of skewness = 3 [mean-median]/[standard deviation]. This definition of skewness is not identical to the more common definition [third central moment]/[standarddeviation]3/2. However,thetwodefinitionsarequalitativelyandquantitativelysimilar for unimodal distributions with moderate asymmetry.

24 D. H.KIM holds, where we have used the fact that the mean of the risk-neutral distribution is the current futures rate. Since (61) F(K) = erfτP′(K) = 1+erfτC′(K), we can check the condition (60) by evaluating the first derivatives P′ and C′ at F using the piecewise-cubic Hermite interpolation scheme for P(K) and C(K). Figure 5 shows F(F)’s based on puts and calls. It can be seen that F(F) can move below 1/2; indeed, where this occurs is quite close to where the skewness measures intersect the 0-line in Figure (a). An example of negatively skewed pdf (April 30, 2007) is shown in Figure 6. The B-MN and MLN models produce fairly similar pdf (except for a small “wiggle” in the MLN pdf), and both clearly show a notable negative skew. As discussed earlier, since zero-cdf models like the QG3 model cannot produce a negatively skewed distribution, the estimated pdf from this model is a nearly symmetric distribution instead. Not 4 B−MN (a) 3 MLN 2 1 0 −1 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 0.7 (b) put call 0.6 0.5 0.4 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 Figure 5: (a) Model-implied skewness from the B-MN and MLN pdfs. (b) Alternative test of skewness based on the Pearson measure.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 25 0.8 B−MN 0.7 MLN QG3 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 1 2 3 4 5 6 7 8 Figure 6: Pdfs on April 30, 2007 surprisingly, pricing errors of the QG3 and CIR2 models tend to be particular large during the times of negatively skewed pdfs (as can be seen in Figure 4b). 5.5. Financial market turmoil of 2008. In early 2008, the turmoil in financial marketsthathadbeguninmid-2007intensified, posingconsiderablesystemic riskandraising much concern about the broader economy. This lead to a substantial lowering of the federal funds target rate (and of the market’s expectation of the future target rate). In March, the 1-year-ahead eurodollar futures rate had dropped to levels as low as 1.9%. Amid heightened uncertainty about the path of the interest rate, eurodollar futures options with very low strikes, in particular 0.25% and 0.5%, had begun to trade. This section explores whether the above models, if any, can explain the prices observed for these options. It befits to start with a cautionary remark that the liquidity in these put options was significantly lower than near-the-money options. The 1-year-maturity 0.25%-strike option, in particular, had an open interest of only about 30. However, the market in the 0.5%-strike option was more active, with decent open interest (ranging from 2000 to 8000 in the months of February and March, comparable to the 25 percentile open interest in the put options with strikes σ∗ to 1.5σ∗ away from the current futures rate in Table 2). Therefore, at least the 0.5%-strike option merits a serious look. Figure 7a shows the pricing error Pmodel Pob for the 0.5%-strike option from the − B-MN, MLN, and QG3 models, for every business day from February 4 to April 28, 2008. The MLN and QG3 models (and the CIR2 model, not shown) generated much larger errors than the B-MN model in this period, and these errors were systematically negative, i.e., the QG3 and MLN models tended to underprice the option. These errors were particularly large during the first half of March; this is precisely the period where

26 D. H.KIM the B-MN model indicates a substantial weight of the δ-function piece (w in eq. (52)), 0.02 0.01 (a) 0 −0.01 −0.02 MLN QG3 −0.03 B−MN −0.04 02/01 03/01 04/01 05/01 0.1 0.08 (b) 0.06 0.04 0.02 0 02/01 03/01 04/01 05/01 2 MLN (c) QG3 1.5 B−MN(0%) B−MN*(0.15%) 1 0.5 0 02/01 03/01 04/01 05/01 Figure 7: (a) Pricing errors Pmodel Pob for K = 0.5%. (b) Weight of the Dirac delta function in the B-MN pdf (w in eq. (52)). (c) − The ratio Pmodel/Pob for K = 0.5%. The B-MN ∗ (0.15%) pdf computations were done without using put options with strikes below 1%.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 27 0.5 (a) market 0.4 QG3 MLN 0.3 0.2 0.1 0 0 1 2 3 4 5 6 0.5 (b) market 0.4 B−MN 0.3 0.2 0.1 0 0 1 2 3 4 5 6 Figure 8: Observed put option prices (circle) and call option prices (diamond) and model predictions on March 10, 2008. shown in Figure 7b. Although not shown, the pricing errors for K = 0.25% tell a similar story. Figure 7c shows the ratio Pmodel/Pob for the B-MN, MLN, and QG3 models. Note that the ratio is much smaller than 1 (sometimes orders of magnitude smaller) for the MLN and QG3 models, i.e., these models grossly underpriced the K = 0.5% option. To show a representative day from this period, the fit of put and call option prices on March 10, 2008 is shown in Figure 8. Note that the QG3 model fits the put and call option prices rather poorly; the CIR2 model (not shown) gives similar results. The B-MN model prices put and call options fairly well across strikes. The MLN model also prices put and call option price well, except at the strikes of 0.25% and 0.5%. The model-implied pdfs for this date, shown in Figure 9, help explain this result. The MLN model pdf has a bimodal structure with a pronounced peak in the low interest rate region (centered at about 0.5%) and a broader peak at about 2.5%. This produces a fairly good fit of option prices for most of the strikes, but since the left tails of both peaks have to vanish at 0, it leads to low prices of put options at strikes of 0.5% and 0.25%. The QG3 model also has a pdf that vanishes in the zero limit (recall fromSection 4.1 that the pdf of the 3-factor positive-definite QG model has the behavior f(r) √r, ∝ r 0), thus producing much lower prices of low-strike options than the B-MN model →

28 D. H.KIM 0.4 B−MN 0.35 MLN QG3 0.3 0.25 0.2 0.15 0.1 0.05 0 −1 0 1 2 3 4 5 6 7 Figure 9: Pdfs on March 10, 2008. which has a delta function singularity with weight of about 0.04. It is also interesting to note that the modal implications of the pdfs are quite different, e.g., the QG3 pdf indicates that the most likely value is less than 2%, while the B-MN pdf indicates that the most likely value is larger than 2%. Several alternative computations of the pdf were done to check the robustness of the results. Forexample, putoptionswithstrikesbelow1%weredroppedinthecomputation (10), but this led to very similar results as the original computations. In other words, the good fit of the low-strike put option in the B-MN model is not due to their inclusion in the computation. Another robustness check is with respect to the choice of the lower-boundary; as discussed in Section 5.2, there are reasons to consider a nonzero r. The computation of the B-MN model pdf with r = 0.15% again matched the price of the 0.5%-strike put option reasonably well. Figure 7c shows, with inverted triangle symbols, the Pmodel/Pob ratio for this case (this computation also dropped options with strikes less than 1%); it can be seen that the ratio is still fairly close to 1, and much larger than the values predicted by the QG3 and MLN models.34 Lastly, note that 2002-3 is another interesting period from which the near-zero behavior of the risk-neutral pdf can be studied, but this paper forgoes a treatment comparable to that given to the 2008 episode, for space considerations. 34Tworemarks: (1)NotethatthePmodel/Pob ratiosfortheQG3andMLNmodelsinFigure7cwere computed from the pdfs that assumed r = 0. If these models also used r = 0.15% (via the functional form f(r r)), they would have predicted even smaller ratios. (2) In the case of the 0.25%-strike put − option(notshowninfigures),ther =0.15%lowerboundfortheB-MNmodelledtoaPmodel/Pob ratio that is about half of the r =0 case, but this still represents much larger values than the corresponding numbers from the QG3, CIR2, and MLN models.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 29 6. Concluding remarks The main conclusions of this paper can be summarized as follows. (1) Black’s boundary model can capture the low-strike put option prices in February and March of 2008, which would have seemed puzzlingly high from the viewpoint of the popular zero-cdf term structure models. While some might alternatively argue that the “high” prices of these options are due to low liquidity or irrational pricing, we can say, at least, that the Black’sboundarybehaviorprovides onepossiblerationalexplanationforthem. (2)Wellknown non-negative interest rate models like the multifactor (positive-definite) CIR/QG models not only have difficulty capturing the low-strike option prices in early 2008 but also generate large and biased option pricing errors at other times as well. We have seen this specifically in the context of 2-factor and 3-factor CIR/QG models, but the theoretical arguments presented in Sec 4.1 suggest that having more factors would not change this conclusion. We have also seen that options data contradict the prediction that the pdf is always positively skewed. These constitute strong evidence against well-known term structure models that have zero cdf at the zero boundary. This paper also has broader implications: (a) The encouraging results from the flexible parametric pdf model with Black’s boundary behavior developed here suggests that such a model may be useful for practical applications to other settings, for example, the extraction of risk-neutral pdfs from the Japanese interest rate options market (euroyen futures options). (b) The evidence against zero-cdf models (and evidence for Black’s boundary models) documented here indicates that, in modeling the postwar US economy there is little motivation to search for “special” macroeconomic mechanisms that would generate behaviors in which the short rate always stays above zero or “bounces off” from zero.35 (c) While the tractability of the affine and QG models may provide some incentive to consider the positive-definite versions of these models as approximate models (even if the true boundary behavior were Black’s boundary behavior), the results in this paper cautions that they would be poor approximations, at least in certain contexts. (d) This paper provides some justification for favoring models that can have negative nominal interest rates. It may appear counter-intuitive to deliberately ignore a clearly reasonable constraint about nominal interests (positive-definiteness). Indeed, Ahn, Dittmar, and Gallant (2002) and Leippold and Wu (2003) considered it a virtue that the QG model’s parameters can be restricted so as to guarantee the positivity of nominal interest rates, and to my knowledge, all published papers that studied the QG model empirically have imposed such restriction. In so doing, however, they rule out specifications such as (62) r = c+x2 x2 . t 1t − 2t 35Note that macroeconomic models that optimize the standard cost function E [ βi(αx2 + (1 α)π2 )] (the GDP gap x and inflation π being unbounded random t i t+i − t+i t t variables), with the condition r 0 introduced via a Lagrange multiplier, would lead to a Black-like t P ≥ boundary behavior, as opposed to a CIR-like behavior.

30 D. H.KIM Though it allows negative nominal rates, a specification like this accommodates negatively skewed pdfs (as observed in Sec. 5.4). The key point is that such specification is very similar to Black’s model specification (which guarantees non-negative interest) (63) r = max[0,c+x2 x2 ], t 1t − 2t except when rates are close enough to zero for the zero bound to have substantial effects;36 hence most of the times it corresponds to Black’s boundary model. In situations in which the zero bound matters, eq. (62) would no longer be a good approximation of (63), but in this case the positive-definite versions of the QG model would not be a good approximation of (63) either. The empirical design of this paper had the following features: (i) the focus on the risk-neutral measure, and (ii) day-by-day implementations using only options data (and the current futures rate) for that day. Though this design was sufficient for detecting the problems with zero-cdf models (such as the multifactor positive-definite CIR/QG models) and yielded useful techniques for modeling pdfs for low interest rate scenarios, in the future it would be useful to extend the design to also examine the physicalmeasure processes and distributions, and utilize the time-series as well as the crosssectional aspects of both options data and yield curve data, to achieve a more detailed understanding of the near-zero boundary processes. Appendix A. The r 0 behavior of the pdf for → multifactor CIR/QG models This appendix derives the small-r asymptotic expansion of the probability density function f(r) for the distribution r b χ2 + +b χ2 . ∼ 1 ν1,λ1 ··· n νn,λn For this, we will use the Laplace transform of the probability density function ∞ (64) f¯(ω) = e−rωf(r)dr. Z0 The inverse Laplace transform 1 c+i∞ (65) f(r)= erωf¯(ω)dω 2πi Zc−i∞ gives back the pdf from f¯(ω). The Laplace transform of the pdf, i.e. E(e−ωr), is simply related to the characteristic function E(eiωr) as f¯(ω) = fˆ(iω). Therefore, from eq. (45), we have n e−λjbjω/(1+2bjω) (66) f¯(ω)= . (1+2b j ω)νj/2 j=1 Y The basic idea is to expand this in powers of 1/ω, which will then lead to an inverse transform in powers of r. 36How close is “close enoughto zero”may depend on the contextof the problem; e.g., out-of-money put option pricing would be more sensitive to this issue than at-the-money option pricing.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 31 Note that the exponent in the exponential in eq. (66) can be written λ b ω λ 1 1 (67) j j = j 1 +O(1/ω2) −1+2b ω − 2 − 2b ω j (cid:18) j (cid:19) and that ν 1 (68) (1+2b ω)νj/2 = (2b ω)νj/2 1+ j +O(1/ω2) . j j 4b ω (cid:18) j (cid:19) Therefore, exp( Σ λ /2+(Σ λ /4b )1/ω+O(1/ω2)) (69) f¯(ω) = − j j j j j (Π j (2b j )νj/2)ωΣjνj/2(1+(Σ j ν j /4b j )1/ω+O(1/ω2)) exp( Σ λ /2) λ ν = − j j ω−Σjνj/2+ Σ j − j ω−Σjνj/2−1+O(ω−Σjνj/2−2) . j (Π j (2b j )νj/2) (cid:18) (cid:18) 4b j (cid:19) (cid:19) Because the inverse Laplace transform of ω−q is simply rq−1/Γ(q), we obtain immediately eq. (48). Appendix B. Option prices for risk-neutral distribution described by positive linear combination of independent noncentral χ2 random variables This appendix derives an approximate (but very accurate) formula for option prices for the distribution of short rate implied by CIR/QG models. We want to evaluate the expressions (70) C(K) = e−rfτE((r K) I(r K)) − − (71) P(K) = e−rfτE((K r) I(K r)) − − where E is an expectation with respect to the distribution r b χ2 + +b χ2 . Since ∼ 1 ν1,λ1 ··· n νn,λn E(I(K r)) is the cumulative distribution function F(r), we have − (72) C(K) = e−rfτ(E(r I(r K)) K(1 F(K)) − − − P(K) = e−rfτ(KF(K) E(r I(K r))). − − Thus we need to evaluate E(r I()). For this, it is useful to define a new probability measure, · with Radon-Nikodym derivative r/E(r). (This can be done because r is a positive random variable.) In terms of the new measure (denoted with superscript ), we have † (73) E(r I()) = E(r)E†(I()), · · therefore (74) C(K) = e−rfτ(E(r)(1 F†(K)) K(1 F(K)) − − − P(K) = e−rfτ(KF(K) E(r)F†(K)), − where F†(K) denotes the cdf in the new measure. From eq. (50), E(r) is simply E(r) = Σ b (ν +λ ). i i i i

32 D. H.KIM It thus remains to derive expressions for cdfs F(K) and F†(K). This can be done using the saddle-point formula of Lugannani and Rice (1980) for cdfs,37 which uses the cumulant generating function (ω) = logE(eωr): C 1 1 (75) F(K) = Φ(z) φ(z) , − y − z (cid:18) (cid:19) z = 2(ω∗K (ω∗))sgn(ω∗) −C y =pω∗ ′′(ω∗), C where ω∗ is the solution of the equation p ′(ω) = K. C Note that (ω) is related to the characteristic function as (ω) = logfˆ( iω). Thus, from C C − eq. (45), we have n eλibiω/(1−2biω) n λ b ω ν i i i (76) (ω) = log = log(1 2b ω) . i C i=1 (1 − 2b i ω)νi/2 ! i=1(cid:18) 1 − 2b i ω − 2 − (cid:19) Y X ThecumulativedistributionfunctioninthenewmeasureF†(K)isgiven bythesameexpression as that of F(K), except that the cumulant generating function in eq. (75) is replaced C by the cumulant generating function for the new measure, †, which is given by C (77) †(ω) = logE†(eωr)= log(E(eωrr)/E(r)). C Using the identity ∂ ∂ (78) E(eωrr)= E(eωr)= fˆ( iω), ∂ω ∂ω − it is straightforward to show n n λ b ν b (79) †(ω) = (ω)+log i i + i i log( b (λ +ν )). C C (1 2b ω)2 1 2b ω − i i i " i=1(cid:18) − i − i (cid:19) # i=1 X X I have checked the Lugannani-Rice formula for the cumulant generating functions , † C C against the Fourier-inversion formula for cdf (eq. (47)) at several sets of parameter values, and found that in all cases the agreements were excellent across the entire range of K, i.e., both in the tail and center regions of the distribution. One caveat is that the formula can be numerically unstable when K happens to be the mean of the distribution (i.e., K = ′(ω = 0), K = †′ (ω = 0)). However, suchsituationoccursrarelyinouroptimizationproblem( C eq. (10)); C if it does occur, one can find a nearby parameter vector that does not have the problem but hasavery similar pdfcontent, orusetheconventional (slower) formula(47) withcharacteristic functions corresponding to , †. C C 37Lugannani and Rice (1980) derived their formula via a judicious choice of the path of integration for the inverse Fourier (Laplace) transform representation of the cdf. This formula found much use in the statistics literature (e.g., Daniels (1987)); its application to stock option pricing (lognormal-type distributions) was explored by Rogers and Zane (1999).

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 33 Appendix C. Options and futures prices for flexible pdfs with Black’s boundary behavior Thisappendixprovidesformulaethatareusefulfortheimplementationofflexiblepdfmodels with Black’s boundary behavior. It is straightforward to show that the option prices for the normal-mixture shadow rate distribution fS (eq. (53)) and lower boundary r are given by 2 (80) P(K) = e−rfτ[w(K r)+ β σ (k (Φ(k ) Φ(ρ ))+φ(k ) φ(ρ ))] i i i i i − − i − i i=1 X 2 C(K) = e−rfτ β σ (φ(k ) k (1 Φ(k ))), i i i i i − − i=1 X where (81) k (K µ )/σ , ρ (r µ )/σ , i i i i i ≡ − i ≡ − (82) w = β Φ(ρ )+β Φ(ρ ). 1 2 1 2 ∞ The current futures rate F= rf(r)dr is given by −∞ R 2 (83) F = wr+ β (σ φ(ρ )+µ (1 Φ(ρ ))). i i i i − i i=1 X Theoption prices for theGram-Charlier/Edgeworth shadowrate distribution(eq. (54)) and lower boundary r are given by γ γ (84) P(K) = e−rfτ[w(K r)+σ(H p (k;ρ)+ 3 H p (k;ρ)+ 4 H p (k;ρ) )], − 0 6 3 24 4 ··· γ γ (85) C(K) = e−rfτσ(Hc(k)+ 3 Hc(k)+ 4 Hc(k) ), 0 6 3 24 4 ··· where k p (86) H (k;ρ) (k x)H (x)φ(x)dx j ≡ − j Zρ ∞ (87) Hc(k) (x k)H (x)φ(x)dx j ≡ − j Zk (88) k (K µ)/σ, ρ (r µ)/σ. ≡ − ≡ − Using the property of Hermite polynomials that H (x)φ(x) = ( 1)j dj φ(x), it is straightj − dxj forward to show γ γ 3 4 (89) w = Φ(ρ) H (ρ)φ(ρ) H (ρ)φ(ρ) . 2 3 − 6 − 24 −··· Similarly, and using integration by parts, we have (90) Hc(k) = H (k)φ(k) j j−2 p (91) H (k;ρ) = (k ρ)H (ρ)φ(ρ)+H (k)φ(k) H (ρ)φ(ρ) j − j−1 j−2 − j−2

34 D. H.KIM for j 2. The integrals Hc(k) and H p (k) are given by ≥ 0 0 (92) Hc(k) = φ(k) k(1 Φ(k)), 0 − − p (93) H (k;ρ) = k(Φ(k) Φ(ρ))+(φ(k) φ(ρ)). 0 − − The current futures rate is given by F = wr+µ(1 Φ(ρ))+σφ(ρ) − γ γ 3 4 (94) + φ(ρ)(rH (ρ)+σH (ρ))+ φ(ρ)(rH (ρ)+σH (ρ))+ . 2 1 3 2 6 24 ··· References Ahn, D. -H., R. F. Dittmar, and A. R. Gallant (2002), “Quadratic term structure models: theory and evidence”, Review of Financial Studies, 15, pp243-88. Ahn, D. -H., and B. Gao (1999), “A parametric nonlinear model of term structure dynamics,” Review of Financial Studies, 12, pp721-62. Ait-Sahalia,Y. (1996),“Testing continuous-timemodels ofthe spotinterestrate,” Review of Financial Studies, 9, pp385-426. Bank for International Settlements (1999), Estimating and Interpreting Probability Density Functions (Proceedings of the workshopheld at the BIS on 14 June 1999). Baroni-Adesi,G.andR. Whaley (1987),“EfficientanalyticapproximationofAmericanoptionvalues,” Journal of Finance, 42, pp301-20. Beaglehole,D. andM. Tenney (1992),“Correctionsandadditions to ’A nonlinearequilibrium modelof the term structure of interest rates’ ” Journal of Financial Economics, 32, pp345-53. Black, F. (1976), “The pricing of commodity contracts”,Journal of Financial Economics, 3, pp167-79. Black, F. (1995), “Interest rates as options”, Journal of Finance, 50, pp1371-76. Bomfim, A. (2003),“Interestrate as options: Assessingthe markets’view ofthe liquidity trap”,FEDS working paper No 2003-45. Brace, A., D. Gatarek, and M. Musiela (1997), “The market model of interest rate dynamics”, Mathematical Finance 7, 127-55. Breeden, D. T. and R. H. Litzenberger (1978), “Prices of state-contingent claims implicit in option prices”, Journal of Business, 51, 621-51. Cecchetti, S. G. (1988), “The case of the negative nominal interest rates: new estimates of the term structure of interest rates during the Great Depression”,Journal of Political Economy, 96,pp1111-41. Chan,K.C.,G.A.Karolyi,F.A.Longstaff,andA.B.Sanders,“Anempiricalcomparisonofalternative models of the short-term interest rate”, Journal of Finance, 47, pp1209-27. Cox, J.C., J. Ingersoll, and S. A. Ross (1985), “A theory of the term structure of interest rates”, Econometrica, 53, pp385-407.

ZERO BOUND, PDFS, TERM STRUCTURE MODELS 35 Daniels, H. E. (1987),“Tailprobabilityapproximations”,International Statistical Review, 55,pp37-48. Dirac, P. A. M. (1958), Principles of Quantum Mechanics, (4th ed.) Oxford Univ. Press. Duffie, D. and K. Singleton (1997), “An econometric model of the term structure of interest rate swap yields”, Journal of Finance, 52, pp1287-1321. Feldhu˝tter, P., and D. Lando (2007), “Decomposing the swap spread”, working paper. Feller, W. (1952), “The parabolic differential equations and the associated semi-groups of transformations”, Annals of Mathematics, 55, pp468-519. Goldstein, R. S. and W. P. Kierstead (1997), “On the term structure of interest rates in the presence of reflecting and absorbing boundaries”, working paper. Gorovoi, V. and V. Linetsky (2004), “Black’s model of interest rates as options, eigenfunction expansions, and Japanese interest rates”, Mathematical Finance, 14, pp49-78. Ho˝rdahl, P. (2000), “Estimating the implied distribution of the future short-term interest rate using the Longstaff-Schwartz model”, working paper. Ho˝rdahl,P. and D. Vestin (2005),“Interpretingimplied risk-neutraldensities: the role of risk premia”, Review of Finance, 9, pp97-137. Jackwerth, J. C. (2004), Option-Implied Risk-Neutral Distributions and Risk Aversion, The Research Foundation of AIMR. Johnson, N., S. Kotz, and N. Balakrishnan (1994), Continuous Univariate Distributions, Volume 2, Wiley. Karlin, S. and H. M. Taylor (1981), A Second Course in Stochastic Processes, Academic Press. Leippold,M.andL.Wu(2003),“Designandestimationofquadratictermstructuremodels”,European Finance Review, 7, pp47-73. Longstaff, F. A. (1992), “Multiple equilibria and term structure models”, Journal of Financial Economics, 32, pp333-44. Longstaff, F. A. and E. S. Schwartz (1992), “Interest rate volatility and term structure: a two-factor general equilibrium model”, Journal of Finance, 47, pp1259-82. Lugannani, R. and S. Rice (1980), “Saddlepoint approximations for the distribution of the sum of independent random variables”, Advances in Applied Probability, 12, pp475-90. McManus, D. J. (1999), “The information content of interest rate futures options”, working paper. Rogers,L.C.G.andO.Zane(1999),“Saddlepointapproximationstooptionprices”,Annals of Applied Probability, 9, pp493-503. Shimko, D. (1993), “Bounds of probability”, RISK, 6(April), 33-37. Stuart, A. and K. Ord (1994), Kendall’s Advanced Theory of Statistics, Volume 1, Halsted Press.

36 D. H.KIM Ueno, Y., N. Baba,and Y. Sakurai(2006), “The use of the Black model of interest rates as options for monitoring the JGB market expectations”, Bank of Japan Working Paper Series.