Keywords: Zero bound, term structure of interest rates, options, probability density functions
Abstract:
This paper points out that several known ways of modeling nonnegative nominal interest rates lead to different implications for the riskneutral 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 nonzero at the zero boundary, while models like the CIR and positivedefinite quadraticGaussian (QG) models have a zero cdf at the boundary. Eurodollar futures options data are found to favor Black's boundary models: the CIR/QG models, even multifactor versions, have difficulty capturing option prices accurately not only in low interest rate 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 wellknown models whose cdf is zero at the zero boundary, namely that the riskneutral pdf is always positively skewed.
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 riskneutral distribution of the underlying security prices, not just mean and variance. Policy makers and market participants have utilized this fact to deduce the riskneutral distribution of shortterm 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 applications to various settings.^{1}However, this literature so far has not made much connection with the vast literature on dynamic term structure models, i.e., relatively little work has been done to investigate what the optionimplied riskneutral distribution tells about 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 zeroboundary behavior was proposed by Fischer Black in his last paper (1995), in which the nominal short rate is modeled as , where is a "shadowrate" process that can go below zero. As we shall 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 CIRlike, Blacklike, 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 Blacktype models. However, some might argue that this conclusion does not necessarily carry over to the postwar US economy.^{3} Indeed, nonnegative 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 quadraticGaussian (QG) model, and the nonaffine 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 outofmoney options can help, as the models may imply strong enough differences in riskneutral 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 options data from this episode tell about term structure models. Beyond fundamental theoretical interest, modeling of riskneutral 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 wellknown and tractable term structure models that respect the nonnegativity condition for nominal interest rates (e.g., the multifactor positivedefinite CIR/QG models) for the riskneutral distribution of the short rate. In particular, this paper derives the asymptotic form of their riskneutral 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 riskneutral distribution implied by the multifactor CIR model and (positivedefinite) QG model, and explores how well these models perform in matching the observed outofmoney 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 Blacktype boundary models for option prices and optionimplied pdfs. After deriving the riskneutral pdf of the simplest Black's boundary model ("BlackVasicek model"), this paper develops flexible parametric models of riskneutral 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 (normalmixture 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 1yearmaturity 0.5%strike eurodollar futures option in February and March of 2008 that agree fairly well with actual settlement prices. On the other hand, the traditional lognormalmixture pdf model and 2factor and 3factor (positivedefinite) CIR/QG models substantially underprice this "lowstrike" 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 riskneutral pdf extraction techniques. It also introduces a mathematical object called the Dirac function, which seldom appears in finance literature but is needed for the discussion of the pdfs with Black's boundary behavior. Section 3 examines the pdfs and option prices in onefactor models, comparing Black's boundary model with CIR/QG models. Section 4 examines richer pdf models and develops useful techniques for them; the application of these models to actual options data and the empirical results are discussed in Section 5. Section 6 concludes.
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 socalled CEV model^{6}
(1) 
(2) 
In the accessible boundary case ( , or and ), 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 stay at zero 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 nonnegative 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 , the BGM model, and Ahn and Gao (1999)'s nonaffine model. The latter include the CIR model with and (positivedefinite) quadraticGaussian 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 "zerocdf models", for brevity.
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 ; this means that the probability density function of the conditional distribution, , is infinite at . Mathematically,
(3) 
(4) 
Using the function, the pdf of Black's boundary models can be expressed as^{10}
The cumulative distribution function of the short rate, , takes the form
(7) 
Consider a eurodollar futures option with strike and maturity .^{11}Let denote the futures rate at the maturity
of the option, and let
be the riskneutral pdf of . The put option price
and the call option price
can be expressed in terms of as
Taking successive derivatives of both sides of eq. (8) gives a simple formula for the pdf,
There are many possibilities for parametrizing . The general wisdom is that a flexible parametric form with 4 to 6 parameters is sufficient to match observed option prices at various strikes.^{14} A parametrization substantially richer than this can often lead 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}
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 as , the form (11) can have difficulty when the actual distribution contains a function at . As noted earlier (eq. (5)), the delta function 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), but in such a twocomponent mixture model this may incur a substantial cost in the fit of other aspects of the pdf. Section 4.2 develops parametric forms for 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 distributions in the riskneutral measure (which determines pricing), and not the physical (real world) measure.^{16}Therefore, henceforth we shall often drop the adjective "riskneutral" for brevity.
Let us now derive the pdf and option prices for the simplest Black's boundary model, which has the shadow rate described by the onefactor Vasicek process:
The conditional distribution of this model, i.e., the transition density , can be easily obtained if one thinks in terms of the shadow rate process rather than the process. Note that if , we have ; hence in this case, the distribution of equals that of . On the other hand, if , we have . Thus all scenarios in which is nonpositive maps to a single point . Therefore, we have
(15) 
The transition density function for the Vasicek process is wellknown: Since the stochastic differential equation (13) has the solution
Thus, suppressing time indices and simplifying notations, we have the conditional distribution of the short rate given by
The prices of maturity put/call options with strike are straightforward
to evaluate. Substituting eq. (19) into eq. (8), we have
(21)  
(22) 
Let us now consider two wellknown and tractable positive interest rate models: the CIR model and the QG model. The CIR model is given by
Consider now the positivedefinite 1factor QG model
(27)  
The pdf of the distribution (CIR/QG models) is given by
(28) 
To get a feel for the qualitative difference between Black's boundary models and zerocdf (CIR/QG) models, it is instructive compare the behavior of the pdf in the limit.
For the CIR/QG model, expanding the expression (29) in powers of , using the wellknown series expansion of the modified Bessel function of the first kind
(30) 
(32) 
From eq. (31) it can be seen that depending on , we can expect qualitatively different behavior of :
(33) 
Note also that since for the QG model, we have for the QG model. As the zero boundary is accessible in the QG model (since in eq. (26) is unrestricted, it can reach zero), it shows a qualitatively similar small behavior as the CIR model with violated Feller condition. It should be noted, however, that an accessible zero boundary does not always imply , as we shall see with the case of the 3factor QG model (Sec. 4.1). Note also that even in the case (onefactor CIR model with violated Feller condition and the QG model), the cdf is zero at , since in the limit.
The small behavior of the pdf in eq. (31) implies that the put option price for small strike is given by
These asymptotic behaviors of CIR/QG models can be compared with the those of the BlackVasicek model:
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 BlackVasicek model with parameters in eq. (19) given by , which imply . It also shows, in thin dashed line, the pdf of the model with (corresponding to the CIR model satisfying the Feller condition), whose and were determined such that the model matches the mean and variance of the BlackVasicek model;^{22} requiring the means and variances of the distributions be similar amounts to having comparable atthemoney option prices. Figure 2b shows the pdf of the model with (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 , 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 distribution with (i.e., Figure 2a).
Despite the rough similarity in the shape of the pdf of the BlackVasicek model and the model, the modelimplied put prices are drastically different, with the BlackVasicek model producing far larger numbers, as shown in Figure 2c. This is because the singularity in the model is a much weaker singularity than the function singularity in the BlackVasicek 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 BlackVasicek model shows a linear dependence on , as predicted by eq. (35).
Although the behavior of the pdfs and put option prices derived above for the onefactor BlackVasicek 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 's than zerocdf models (with comparable mean and variance) 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 nonnegligible. 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 model (inaccessible boundary) and the model (accessible boundary) are quite similar, beyond the very small region.
The onefactor 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 zerocdf models and Black's boundary models.
A nonnegative, multifactor version of the CIR model can be constructed by adding up independent CIR factors:
(36)  
Because 's are independent, the results about the 1factor CIR model carries over easily; it is straightforward to show that the distribution of is given by
Let us now consider the general factor positivedefinite QG model,
It is straightforward to show that the conditional distribution of the short rate in the positivedefinite QG model also takes the form (37),^{25} with
(41)  
(42)  
(43) 
(44) 
Models like the multifactor CIR model have been criticized for the restrictive nature of the factor correlation structure. However, the positivedefinite 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 emphasizing that not only the multifactor CIR model but also the positivedefinite QG model has a short rate distribution given by a positive linear combination of independent noncentral random variables, and that this holds no matter how flexible the model is and no matter how many factors it has.
The pdf and cdf of the distribution for the multifactor CIR/QG model (eq. (37)) do not have simple closedform expressions in general, but the characteristic function has a simple expression
Let us now consider the small behavior of . Appendix A shows that
(49) 
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 (
). The mean, variance, and the third central moment of this distribution are easily derived from the cumulants of the noncentral variables (see, e.g., Johnson, Kotz, Balakrishnan (1994, p447)):
More generally speaking, one cannot create a negatively skewed distribution by forming a positive linear combination of independent, positiveskewed random variables. Therefore, a model like
(51) 
This section develops flexible parametric forms for the riskneutral pdf that are consistent with Black's boundary behavior. The earlier discussion with the BlackVasicek model points to the way: write down a flexible form for the shadow rate distribution that has some weight below zero, and take
The distribution can be modeled with similar degree of parsimony as the pdfs that are in use in the extant literature. A natural candidate for is the "mixture of normals" form, i.e.,
These forms, which are richer than that of the BlackVasicek 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 , and F for the mixtureofnormals shadow rate model (eq. (53)) and the GramCharlier/Edgeworth shadow rate model (eq. (54)) are provided in Appendix C.
To examine how well zerocdf 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:


The QG3 model covers all possible 3factor positivedefinite QG models, but it also includes a special case of the 3factor CIR model ( ). The CIR2 model covers all possible cases of the 2factor CIR models, but it also nests the 2factor QG model. This model also nests the LongstaffSchwartz (1992) model, which was regarded by at least Hrdahl (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 BMN 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 ( for QG3) are linked to the elementary parameters of the term structure model and the state variables .^{28}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 BMN 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., and . Empirically, Duffie and Singleton (1997) obtained a negative for the multifactor CIR model, but in that case it is no longer a nonnegative interest rate model. Ahn, Dittmar, and Gallant (2002) obtained a positive in their estimation of the QG model; in this case, there is no probability weight in the region . 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.
For the empirical investigation, we shall confine attention to the eurodollar futures option with maturity of about 1 year. The 1year 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). The prefactor in eq. (8) was computed using the 1year LIBOR. The CME eurodollar futures options are American options. I have converted their price into corresponding European price by subtracting the early excercise premium computed based on BaroniAdesi 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 sawtooth behavior of the time series of option maturities (oscillating between approximately 7/8 and 9/8 year) for the monthly sample.
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 , and the call strike rate that corresponds to the futures rate plus , where the proxy for the standard deviation of the pdf is computed as
(57) 
Table 2 provides some statistics regarding the liquidity of options proxied by open interest, for strikes 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 standard deviation) is fairly active; e.g., the median open interest in the put option is more than a third of that of the option.
Lastly, as a caveat, note that since the here is taken to be the 3month LIBOR, the "zero bound" (more precisely, the lower bound) is not exactly at zero, as the 3month LIBOR contains a premium for credit/liquidity risk, besides some differential between the 3month rate and the instantaneous rate. This means that the proper location of the delta function for the Black's boundary model pdf ( ) should be a small but nonzero number. Furthermore, the delta function itself would be smudged out a bit, in other words,
It is difficult to get a precise sense of for the 3month LIBOR, but the recent Japanese experience can provide some guidance: even at the height of Bank of Japan's Quantitative Easing Program, the 3month 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 is at 0, but I have also computed the pdf for the BMN model with 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., ), but at a very low strike (e.g., 25 basis points), there may be a more visible difference between vs .
To get a sense of the overall fit across strikes, it is useful to look at the rootmeansquare errors, i.e.,
(59) 
Figure 4 shows the time series of RMSEs based on the four models (BMN, MLN, CIR2, QG3) for the 19982008 period. As can be seen in the top panel, the BMN model (Black's boundary model) and the MLN model (conventional benchmark) performed similarly well, the MLN model producing a somewhat smaller RMSE in 20032004 and larger RMSE in 20072008. Note that during most of the data period, with the exceptions around 2003 and 2008, the distributions were sufficiently away from the zero bound (as can be seen in Figure 3), thus the BMN model was simply the "mixture of normals" model (i.e., ); the results indicate that there is little practical difference between the normalmixture model and the lognormalmixture 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 yaxis 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 (strikespecific) pricing errors, and . Table 3 displays summary statistics for the pricing errors for put options with strikes and units below F and call options with strikes and 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 75percentile value of the pricing error for the put option at the strike of F is negative; this negative bias means that the QG3 and CIR2 models tend to systematically underprice the fartheroutofmoney put options. On the other hand, they tend to overprice the nearthemoney call options.
In sum, the pdfs based on zerocdf term structure models (CIR2 and QG3 models) perform notably worse than the MLN and BMN models despite having one more parameter, suggesting that the kind of restrictions on the shape of the riskneutral distribution implied by these models are not well supported by data.
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 , while a positively skewed distribution has .^{33} Since the cumulative distribution function at the median is 1/2 by definition, one can check whether the skewness is always positive by examining whether the condition
(61) 
An example of negatively skewed pdf (April 30, 2007) is shown in Figure 6. The BMN 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 zerocdf models like the QG3 model cannot produce a negatively skewed distribution, the estimated pdf from this model is a nearly symmetric distribution instead. Not 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).
In early 2008, the turmoil in financial markets that had begun in mid2007 intensified, posing considerable systemic risk and raising 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 1yearahead 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 nearthemoney options. The 1yearmaturity 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 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 for the 0.5%strike option from the BMN, 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 BMN 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 the BMN model indicates a substantial weight of the function piece ( in eq. (52)), shown in Figure 7b. Although not shown, the pricing errors for tell a similar story. Figure 7c shows the ratio for the BMN, 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 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 BMN 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 modelimplied 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 from Section 4.1 that the pdf of the 3factor positivedefinite QG model has the behavior , ), thus producing much lower prices of lowstrike options than the BMN model 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 BMN 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. For example, put options with strikes below 1% were dropped in the computation (10), but this led to very similar results as the original computations. In other words, the good fit of the lowstrike put option in the BMN model is not due to their inclusion in the computation.
Another robustness check is with respect to the choice of the lowerboundary; as discussed in Section 5.2, there are reasons to consider a nonzero . The computation of the BMN model pdf with again matched the price of the 0.5%strike put option reasonably well. Figure 7c shows, with inverted triangle symbols, the 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 20023 is another interesting period from which the nearzero behavior of the riskneutral pdf can be studied, but this paper forgoes a treatment comparable to that given to the 2008 episode, for space considerations.
The main conclusions of this paper can be summarized as follows. (1) Black's boundary model can capture the lowstrike put option prices in February and March of 2008, which would have seemed puzzlingly high from the viewpoint of the popular zerocdf 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's boundary behavior provides one possible rational explanation for them. (2) Wellknown nonnegative interest rate models like the multifactor (positivedefinite) CIR/QG models not only have difficulty capturing the lowstrike 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 2factor and 3factor 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 wellknown 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 riskneutral pdfs from the Japanese interest rate options market (euroyen futures options).
(b) The evidence against zerocdf 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 positivedefinite 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 counterintuitive to deliberately ignore a clearly reasonable constraint about nominal interests (positivedefiniteness). 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
The empirical design of this paper had the following features: (i) the focus on the riskneutral measure, and (ii) daybyday implementations using only options data (and the current futures rate) for that day. Though this design was sufficient for detecting the problems with zerocdf models (such as the multifactor positivedefinite 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 timeseries as well as the crosssectional aspects of both options data and yield curve data, to achieve a more detailed understanding of the nearzero boundary processes.
This appendix derives the small asymptotic expansion of the probability density function for the distribution .
For this, we will use the Laplace transform of the probability density function
(64) 
(65) 
The Laplace transform of the pdf, i.e. , is simply related to the characteristic function as . Therefore, from eq. (45), we have
Note that the exponent in the exponential in eq. (66) can be written
(67) 
(68) 
Therefore,
(69)  
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)  
(71) 
(72)  
(73) 
(74)  
It thus remains to derive expressions for cdfs and
. This can be done using the saddlepoint formula of Lugannani and Rice (1980) for cdfs,^{37} which uses the cumulant generating function
:
Note that is related to the characteristic function as . Thus, from eq. (45), we have
(76) 
The cumulative distribution function in the new measure is given by the same expression as that of , except that the cumulant generating function in eq. (75) is replaced by the cumulant generating function for the new measure, , which is given by
(77) 
(78) 
(79) 
I have checked the LugannaniRice formula for the cumulant generating functions against the Fourierinversion 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 , i.e., both in the tail and center regions of the distribution. One caveat is that the formula can be numerically unstable when happens to be the mean of the distribution (i.e., , ). However, such situation occurs rarely in our optimization problem (eq. (10)); if it does occur, one can find a nearby parameter vector that does not have the problem but has a very similar pdf content, or use the conventional (slower) formula (47) with characteristic functions corresponding to .
This appendix provides formulae that are useful for the implementation of flexible pdf models with Black's boundary behavior.
It is straightforward to show that the option prices for the normalmixture shadow rate distribution (eq. (53)) and lower boundary
are given by
(80)  
(81)  
(82) 
(83) 
The option prices for the GramCharlier/Edgeworth shadow rate distribution (eq. (54)) and lower boundary
are given by
(84)  
(85) 
(86)  
(87)  
(88) 
Using the property of Hermite polynomials that , it is straightforward to show
(89) 
(90)  
(91) 
(92)  
(93) 
(94) 
Ahn, D. H., R. F. Dittmar, and A. R. Gallant (2002), "Quadratic term structure models: theory and evidence", Review of Financial Studies, 15, pp24388.
Ahn, D. H., and B. Gao (1999), "A parametric nonlinear model of term structure dynamics," Review of Financial Studies, 12, pp72162.
AitSahalia, Y. (1996), "Testing continuoustime models of the spot interest rate," Review of Financial Studies, 9, pp385426.
Bank for International Settlements (1999), Estimating and Interpreting Probability Density Functions (Proceedings of the workshop held at the BIS on 14 June 1999).
BaroniAdesi, G. and R. Whaley (1987), "Efficient analytic approximation of American option values," Journal of Finance, 42, pp30120.
Beaglehole, D. and M. Tenney (1992), "Corrections and additions to 'A nonlinear equilibrium model of the term structure of interest rates' " Journal of Financial Economics, 32, pp34553.
Black, F. (1976), "The pricing of commodity contracts", Journal of Financial Economics, 3, pp16779.
Black, F. (1995), "Interest rates as options", Journal of Finance, 50, pp137176.
Bomfim, A. (2003), "Interest rate as options: Assessing the markets' view of the liquidity trap", FEDS working paper No 200345.
Brace, A., D. Gatarek, and M. Musiela (1997), "The market model of interest rate dynamics", Mathematical Finance 7, 12755.
Breeden, D. T. and R. H. Litzenberger (1978), "Prices of statecontingent claims implicit in option prices", Journal of Business, 51, 62151.
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, pp111141.
Chan, K. C., G. A. Karolyi, F. A. Longstaff, and A. B. Sanders, "An empirical comparison of alternative models of the shortterm interest rate", Journal of Finance, 47, pp120927.
Cox, J.C., J. Ingersoll, and S. A. Ross (1985), "A theory of the term structure of interest rates", Econometrica, 53, pp385407.
Daniels, H. E. (1987), "Tail probability approximations", International Statistical Review, 55, pp3748.
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, pp12871321.
Feldhtter, P., and D. Lando (2007), "Decomposing the swap spread", working paper.
Feller, W. (1952), "The parabolic differential equations and the associated semigroups of transformations", Annals of Mathematics, 55, pp468519.
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, pp4978.
Hrdahl, P. (2000), "Estimating the implied distribution of the future shortterm interest rate using the LongstaffSchwartz model", working paper.
Hrdahl, P. and D. Vestin (2005), "Interpreting implied riskneutral densities: the role of risk premia", Review of Finance, 9, pp97137.
Jackwerth, J. C. (2004), OptionImplied RiskNeutral 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. and L. Wu (2003), "Design and estimation of quadratic term structure models", European Finance Review, 7, pp4773.
Longstaff, F. A. (1992), "Multiple equilibria and term structure models", Journal of Financial Economics, 32, pp33344.
Longstaff, F. A. and E. S. Schwartz (1992), "Interest rate volatility and term structure: a twofactor general equilibrium model", Journal of Finance, 47, pp125982.
Lugannani, R. and S. Rice (1980), "Saddlepoint approximations for the distribution of the sum of independent random variables", Advances in Applied Probability, 12, pp47590.
McManus, D. J. (1999), "The information content of interest rate futures options", working paper.
Rogers, L. C. G. and O. Zane (1999), "Saddlepoint approximations to option prices", Annals of Applied Probability, 9, pp493503.
Shimko, D. (1993), "Bounds of probability", RISK, 6(April), 3337.
Stuart, A. and K. Ord (1994), Kendall's Advanced Theory of Statistics, Volume 1, Halsted Press.
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.