A Bayesian regression spline approach to estimation of the term structure of interest rates.

ABSTRACT

This paper presents a case study of interest rate term structure estimation of U.S. Treasury bonds and A T& T corporate bonds from April 1994 to December 1995. We first adopt a Bayesian regression spline model to estimate the term structure of risk-free Treasury bonds where the number and location of the spline knots are adaptively selected using the reversible jump Markov chain Monte Carlo algorithm. We then develop a hierarchical Bayesian approach to estimate the corporate term structure, "borrowing strength" from the estimated Treasury term structure. This is necessary due to only a few corporate bonds are available on any given trading day. An appealing feature of this Bayesian hierarchical framework is that the knowledge of positive credit spread can be naturally incorporated into the model with informative priors. The small sample size of the corporate bonds poses little difficulty in inference.

Keywords: Coupon bonds; Credit spreads; Forward rates; Reversible-jump Markov chain Monte Carlo; Yield curve

1. INTRODUCTION

The term structure of interest rates describes the relationship between the interest rates of a bond and the bond's maturity. It plays a vital role in both macroeconomics and finance. For example, the term structure of government bonds contains information about current macroeconomic conditions and market participants' expectation of future economic conditions. Many central banks use the term structure as an indicator of their monetary policy. Other important applications of the term structure of interest rates are described in Li and Yu (2005). A credit derivative is a financial product that transfers credit risk (risk due to a party's possible inability to meet its obligations) from one party to another. Given the explosive growth in the credit derivatives market, the term structure of corporate bonds is also becoming increasingly important. Unlike risk-free government bonds, corporate bonds bear credit risk in that a company may default on its debt and are less liquid. Therefore, the interest rate of corporate bonds is usually higher than that of government bonds and the difference is called spread (risk premium). Due to increased trading in instruments with credit risk, credit derivatives were created to partially or fully offset the credit risk of a business deal. In fact, the size of the credit derivatives market in 2001 reached around 835.5 billion dollars. Some of the commonly traded credit derivatives include default swaps, credit spread options, credit linked notes, etc. Bielecki and Rutkowski (2002) give a detailed account of credit derivatives. Pricing models for both corporate debt and credit derivatives require the corporate term structure as input (see Jarrow and Turnbull, 1995; Duffie and Singleton, 1999). Moreover, these term structures can be used to assess credit quality for certain risk management procedures. Value at risk (VAR) computations, bond portfolio management, corporate loan considerations, and FDIC insurance premium calculations (see FDIC, 2000) all require credit quality assessment. Thus, obtaining accurate estimates of both Treasury and corporate term structure is essential.

The term structure of interest rates can be determined by any of the discount function D(0,T), the yield curve y(0,T), or the forward rate f(0,T), where T is time to maturity from today (time 0). To simplify the notation, we consider current time to be fixed at 0 and omit 0 from this point on. The relationships among these three functions are:

D(T) = exp{-Ty(T)} = exp {-[[integra].sup.T.sub.0](s)ds}. (1)

Availability of any one of the three determines the rest. The discount function D(T) represents today's price of a zero-coupon bond (paying no interest or principal until maturity) that pays one dollar at maturity time T. The yield curve y(T) = [[integral].sup.T.sub.0]f(s)ds/T gives the average of the forward rates between today (time 0) and the maturity date T. The forward rate f(s) gives the marginal return, which is also the rate one can lock in today for future borrowing or lending at time s. To estimate the term structure, a spline may be used to model the discount function, the yield curve, or the forward rate. However, estimating the term structure of the forward interest rates f(s) appears to be the best approach, as argued in Li and Yu (2005).

The forward interest rates in term structure are not observable but implied by the price of the bonds. The relationship between the forward interest rate f(s) (see equation 1) and the date s is the term structure we are interested in. The "empirical forward rate" (see Jarrow et al., 2004) can be found through f = [DELTA]{-log(D)}/[DELTA]s, where [DELTA] is the differencing operator. It provides a rough estimate of the unobservable forward interest rate. Figure 1 and Figure 2 plot the price and the empirical forward rates versus years to maturity s for the US Treasury STRIPS (Separate Trading of Registered Interest and Principal of Securities) on December 31, 1995. The objective is to estimate the forward rate based on the bond price.

[FIGURES 1-2 OMITTED]

Many methods have been developed for term structure estimation. McCulloch (1971, 1975), Shea (1984) estimate the discount function with regression spline methods. Vasicek and Fong (1982) estimate the discount function with exponential splines. Chambers et al. (1984) estimate the yield curve with polynomials. Linton et al. (2001) present kernel methods in term structure estimation and establish asymptotic properties. Adams and Van Deventer (1994), Fisher et al. (1995) estimate the Treasury term structure using a smoothing (penalized) spline model. Jarrow et al. (2004) and Li and Yu (2005) use penalized splines (see Ruppert et al., 2003) to estimate both Treasury and individual corporate term structures, where the forward rate curve is estimated by minimizing the model mean squared error and putting a penalty on the spline coefficients.

This paper presents a case study of term structure estimation of US Treasury and corporate bonds using a Bayesian regression spline (BRS) approach. Smoothing is achieved through adaptively selecting the number and location of the knots using the reversible jump Markov chain Monte Carlo (MCMC) method. The term structure of US Treasury STRIPS is estimated first. A hierarchical Bayesian model is then developed to estimate the spread between the corporate term structure and the estimated Treasury term structure. This "borrowing strength" from the Treasury term structure is necessary due to the fact that the sample size of corporate bonds available on a given day is 4 or 5 on average while there are numerous Treasury STRIPS securities on any trading day. Treasury securities are considered default-free and less risky than corporate securities given the taxing power of the federal government. Thus, this credit spread should be positive and the Bayesian model we propose accounts for this naturally with an informative prior. On the other hand, within the classical framework, forcing the credit spread to be positive is equivalent to adding n constraints to a nonlinear minimization problem, where n is the number of corporate bonds available. It is possible to solve this constrained nonlinear optimization problem. However, the small sample size of corporate bonds may make it difficult to justify the classical statistical inference procedures as they rely on the assumption of large samples. The Bayesian approach does not assume large samples and statistical inference can be performed based on small samples.

The remainder of the paper is organized as follows. The data are described in Section 2. A Bayesian approach to estimating the Treasury term structure with regression splines is presented in Section 3. We then illustrate the Bayesian estimation of the corporate term structure with informative priors in Section 4. Section 5 presents a case study of term structure estimation of US Treasury STRIPS and AT&T bonds traded from April 1994 to December 1995. Some concluding remarks are provided in Section 6.

2. DATA

The University of Houston Fixed Income database contains over 28,000 instruments, including publicly traded non-convertible debt with principal value over one million dollars. Warga (1995) gives a detailed description of the database. The US Treasury STRIPS data contain current dates, issue dates, maturity dates, and the market prices. The bonds included in the Lehman Brothers Bond Indices are reported with month-end flat prices, accrued interest, coupon, yields, current date, issuance date, maturity date, S&P ratings, Moody's and option-like features. The market price of a corporate bond equals the quoted flat price plus the accrued interest. For the case study, we will apply the proposed Bayesian estimation procedures to US Treasury STRIPS and AT&T Bonds from April 1994 to December 1995 in the data base. Table 1 shows the market prices of all five AT&T bonds available on December 31, 1995.

There are usually a few individual corporate bonds available in a given month. Jarrow et al. (2004) observed that on average there were only 4.3 bonds available per month during the period of April 1994 to December 1995. It is difficult to obtain a meaningful estimate of the corporate term structure based on these few observations. On the other hand, the average number of US Treasury STRIPS per month during the same period was 117, ranging from 115 to 120. "Borrowing strength" from other sources such as US Treasury STRIPS becomes necessary in estimating the term structure of the corporate bonds.

For estimation, we need to convert the time-to-maturity and the coupon payment times to the same unit scale. These can be easily accomplished with the Matlab functions days365() and cfdates(). Days365() counts the number of days between dates based on 365-day year and cfdates0 gives cash flow dates for a fixed-income security. Other formats, such as days360(), can be used if dates are based on 360 days a year.

3. BAYESIAN REGRESSION SPLINE ESTIMATION OF TREASURY TERM STRUCTURE

3.1 Model

The Treasury term structure model is discussed in this section. Let [P.sub.i], i = 1, ..., n, be the market price of the ith bond on date 0 (today). Each bond pays fixed coupons and principal [C.sub.i]([t.sub.i,j]) due on dates [t.sub.i,j], where j = 1, ..., [z.sub.i] and [z.sub.i] is the total number of coupon and principal payments for the ith bond. Thus, [t.sub.i,t] is the first coupon payment date and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the maturity date of bond i. The model price for the coupon bond, [[mu].sub.i], is related to the forward rate f through the discount function:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The forward rate curve f is modeled as splines f(s) = [delta]'B(s), where B(s)is a vector of spline basis functions and [delta] is the coefficient vector. We use the d-th degree power basis:

f(s) = [[delta].sub.0] + [[delta].sub.1]s + ... + [[delta].sub.d][s.sub.d] + [K.summation over (k=1)][[delta].sub.d+k][(s - [t.sub.k]).sup.d.sub.+]. (3)

Here delta = [[[delta].sub.0], [[delta].sub.1], ..., [[delta].sub.d+K]]' and B(s) = [1, s, ..., [s.sup.d], [(s - [t.sub.1]).sup.d.sub.+], ..., [(s - [t.sub.K]).sup.d.sub.+]]',

where [(s - [t.sub.k]).sub.+] = max(0, s - [t.sub.k]) and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are K spline knots. Power basis has the advantage of being simple and allows easy modeling of the credit spread for corporate term structure in the later portion of the paper. Moreover, setting some polynomial coefficients to zero allows convenient modeling of sub-models. Other basis, such as B-spline (de Boor, 1978), can also be used. However, truncated power basis can be easily transformed during computation to B-spline basis for numerical stability, as in many regression packages (Ruppert et al., 2003, P. 71).

3.2 Reversible Jump MCMC

Denison et al. (1998) apply the reversible jump MCMC (Green, 1995) method to select the number and location of the knots in the linear regression model. After knot locations are selected, the coefficient estimates are computed using the least squares method. The model is linear in coefficients. The term structure model is nonlinear since the spline coefficients appear in the exponential term in the discount function as shown in equation (2). We adopt the approach in Denison et al. (2001) and compute the spline coefficients using a nonlinear regression algorithm once the number and location of the knots are chosen. Thus, the method here is semi-Bayesian in nature. A full Bayesian approach will require imposing priors on spline coefficients. However, conjugate priors are not available for this nonlinear model and computation becomes too burdensome.

The model relating the observed market price [P.sub.i] and the model price [[mu].sub.i], defined in (2), for the ith bond is:

[P.sub.i] = [[mu].sub.i] + [[epsilon].sub.i], (4)

with the zero-mean normal errors [[epsilon].sub.i], i = 1, ..., n. In matrix notations, the model is

P = [mu] + e,

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and n is the total number of Treasury bonds available. Some monotonic transformation h(.), such as the identity or log function, may be applied to the market price and the model price. Often, in economics literature, it is more appropriate to consider the bond price following the lognormal distribution. Then the log transformation log(P) = log([mu])+ [epsilon], or P = [mu]exp([epsilon]) may be used. In addition, for zero-coupon bonds such as the US Treasury STRIPS data, there is no coupon payment and only the principal is due on the maturity date. The model price reduces to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (5)

where L is the principal and [T.sub.i] is time until maturity for the ith bond. An advantage of the log transformation h([micro]) is linear in the spline coefficients in estimation. We found the log transformation gives virtually the same fit as (4). Knot selection is more crucial. Hence, for presentation simplicity, we will keep the notation without transformation.

To simplify the notation, we integrate the forward rate:

[[integral].sup.t.sub.0]f(s)ds = [[integral].sup.t.sub.0][delta]'B(s)ds = [delta]'[[integral].sup.t.sub.0]B(s)ds = [delta]'[B.sup.t](t), (6)

where

[B.sup.i](t) = [[integral].sup.t.sub.0]B(s)ds = [t [t.sup.2]/2 ... [t.sup.d+1]/(d + 1) [(t - [t.sub.1]).sup.d+1.sub.+]/(d + 1) ... [(t - [t.sub.k].sup.d+1.sub.+]/(d + 1)]'. (7)

The regression splines minimize the sum of squared deviations between [P.sub.i] and [[mu].sub.i], where [P.sub.i] is the observed bond price and [[mu].sub.i]([delta]) is the model price:

[(1/n)[n.summation over (i=1)][[P.sub.i] - [[mu].sub.i]([delta])].sup.2]. (8)

The parameters to be estimated are K, [t.sub.1], ..., [t.sub.K], and [[sigma].sup.2] The coefficients [[delta].sub.0], [[delta].sub.1], ..., [[delta].sub.d], [[delta].sub.d+1], ..., [[delta].sub.d+K] are determined from the data (knot locations) using the Gauss-Newton method. We use [theta] to denote the parameter vector: [theta] = [K, [t.sub.1], ..., [t.sub.K], [[sigma].sup.2]). Given the normal errors, the likelihood is

[L.sub.k]([theta]\P) = [[1/([square root of 2[pi][sigma]]).sup.n]exp{- [n.summation over (i=1)][{[P.sub.i] - [[mu].sub.i]}.sup.2]/[2[sigma].sup.2]}.

We follow the general framework by Green (1995) and Denison et al. (1998): including a birth step (add a knot), a death step (delete a knot), and a movement step (move a knot) in the reversible jump MCMC algorithm. We model the forward rate curve with a spline as in (3).

Let K be the number of knots in the model with knots located at: ([t.sub.1], [t.sub.2], ..., [t.sub.K]). Both K and knot locations ([t.sub.1], [t.sub.2], ..., [t.sub.K]) are random variables that need to be estimated. Denison et al. (1998) use the Poisson prior on K. However, Hansen and Kooperberg (2002) convincingly argue that the discrete uniform distribution is a better choice and we adopt it as the prior:

p(K) = (K!([k.sub.max] - K)!/[k.sub.max] !)/([k.sub.max] !)/([k.sub.max] + 1), K [less than or equal to] [k.sub.max],

where [k.sub.max] is the maximum number of knots we choose to be in the model. We did not find [k.sub.max] to play a significant role in the final fit in the data we used in the case study. Candidate location of the knots is at data points. The birth and death probabilities are chosen to be

[b.sub.K] = 0.4 min {1, p(K + 1)/ p(K)}, [d.sub.K+1] = 0.4 min {1,p(K)/ p(K + 1)},

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The move probability is [m.sub.K] = 1 - [b.sub.K] - [d.sub.K]. Any of the three birth, death, and move steps gives a candidate model.

The acceptance probability for the birth step is the product of the likelihood ratio, the prior ratio, the proposal ratio, and the Jacobian. The likelihood functions for both the current model and the candidate model are evaluated using the parameter estimates of [[delta].sub.0], [[delta].sub.1], ..., [[delta].sub.d], [[delta].sub.d+1], ..., [[delta].sub.d+K], and [[sigma].sup.2]. The polynomial and spline coefficients estimates are computed using the Gauss-Newton method. The variance can be treated as an auxiliary parameter and its estimate is generated using a Gibbs sampling step. Using the usual conjugate prior distribution on [[sigma].sup.2], the inverted Gamma distribution IG (a, b), the conditional posterior distribution of [[sigma].sup.2] given the other parameters and data is again an inverted Gamma distribution:

IG(n/2 + a, [((P - [mu])' (P - [mu])/2 + 1/b).sup.-1]). (10)

Here a is chosen to be a small positive number and b is chosen to be a large positive number. Therefore, the prior distribution on [[sigma].sup.2] is the usual diffuse prior and has little impact on the posterior distribution.

Once all parameter estimates are computed, the likelihood ratio can be formed:

N ([[mu].sub.c], [[??].sup.2.sub.c] [I.sub.n]) / N ([mu], [[??].sup.2. [I.sub.n]), (11)

where c denotes the candidate model. This ratio is the Bayes factor, approximately. However, DiMatteo et al. (2001) point out the likelihood ratio in the normal (least squares) model always favors the model with more parameters. DiMatteo et al. (2001) argue that the likelihood ratio based on the least squares methods in Denison et al. (1998) should include a penalty factor as in BIC (Schwarz, 1978). We also replace (11) with BIC to penalize the dimensionality of the model. According to Schwarz (1978), -2log (Bayes Factor) can be approximated by BIC:

BIC = -2log [N ([[mu].sub.c], [[??].sup.2.sub.c] [I.sub.n]) / N ([mu], [[??].sup.2] [I.sub.n])] - ([Param.sub.c] - Param) log n,

where the likelihood is evaluated at the maximum likelihood estimates, [Param.sub.c] is the number of parameters in the candidate model, and Param is the number of parameters in the current model. Thus, the Bayes factor is approximately exp(-BIC/2). This relationship is used to approximate (11) and thus to penalize more complex models.

The prior ratio is

(K + 1)/([k.sub.max] - K). (12)

The proposal ratio is

[b.sub.k]/([k.sub.max] - K). (13)

The Jacobian equals 1 in that the model space is discrete. Thus, the acceptance probability for the birth step is

[alpha] = min {1, exp(- BIC/2)x([d.sub.k+1]/[b.sub.k])}. (14)

Similarly, the acceptance probability for the death step is

[alpha] = min {1, exp(- BIC/2)x([b.sub.k-1]/[d.sub.k])}. (15)

Each collection of knot locations represents a possible model and our task is to choose the correct model. One approach to find the correct model is to choose the mode, i.e., the knot locations that give the maximum likelihood for the model. An alternative approach is to adopt the view of Bayesian model averaging and find the weighted average of the estimated values for f based on a large number of simulations from the predictive distributions. It is well-known that given the predictive squared-error loss function for the model response, the Bayes estimator is the expected response under the predictive distribution. This is equivalent to finding the average of all models considered. Thus, given N reversible jump MCMC samples, the posterior mean of f can be approximated by

(1/N) [N.summation over (i=1)] E{ f | B, P, [[theta].sup.(i)}. (16)

The algorithm is outlined in the next section.

3.3 Algorithm

1. Initially, select [k.sub.0] knot locations on data points.

2. Compute the coefficient estimates using the Gauss-Newton method.

3. Generate a uniform (0,1) random number u:

(i) If u [less than or equal to] [b.sub.k], perform the birth step. Add a knot by randomly picking a knot from the candidates. Perform step 2 and compute the acceptance probability (14). Accept the birth step with the acceptance probability.

(ii) If [b.sub.k] < u [less than or equal to] [b.sub.k] + [d.sub.k], perform the death step. Randomly pick a knot from the current model and delete it. Perform step 2 and compute the acceptance probability (15). Accept the deletion with the acceptance probability.

(iii) If u > [b.sub.k] + [d.sub.k], perform the move step. Relocate a knot and compute the acceptance probability (the likelihood ratio), which is (11). The move is accepted with the acceptance probability.

4. Repeat the steps 2 and 3 above until the mean-squared error MSE = (1/n) [n.summation over (i=1)] [([P.sub.i] - [[mu].sub.i]).sup.2], converges.

5. Find the forward rate curve either using the mode [??] or using the model averaging estimate (16).

4. HIERARCHICAL BAYESIAN ESTIMATION OF THE CORPORATE TERM STRUCTURE WITH INFORMATIVE PRIORS

We have estimated the term structure of Treasury bonds using Bayesian regression splines. In this section, we illustrate Bayesian estimation of the term structure of individual corporate bonds with informative prior distributions on the spread parameters. This term structure is needed in pricing and hedging corporate debt, among other applications.

Since only a few bonds of a corporation are traded on any given day (see Section 2), it is difficult to estimate the term structure based on these few observations. Corporate bonds are assigned credit grades such as AAA, BBB, etc. by rating agencies. In some situations, one may estimate the term structure of many bonds with the same credit grade to alleviate the small sample size problem. However, in many situations, it is essential to know the term structure of the individual corporate bond. To estimate the individual corporate term structure, one approach is to "borrow strength" from the estimated Treasury term structure from Section 3. We "borrow strength" by adding a credit spread to the estimated Treasury term structure [f.sub.Tr] (s) = [??]' [B.sup.1] ([t.sub.i,j]) where [??] is a vector of coefficient estimates from Section 3:

[f.sub.C] (s) = [f.sub.Tr] (s) + spread. (17)

This credit spread represents the excess return from a corporate bond over the return from an equivalent Treasury bond. Thus, the resulting term structure is in fact the term structure of credit spreads. There is one Treasury term structure curve on a given day with a large sample size but many corporate bonds with very small sample sizes each. "Borrowing strength" from the Treasury term structure to estimate the corporate term structures is both desired and natural. The Bayesian approach facilitates this "borrowing strength."

Given the small sample size of individual corporate bonds, we believe constant, linear, and quadratic spreads are sufficient. For a constant spread, [f.sub.C] (s) = [[??].sub.Tr] (s) + [[beta].sub.0] = [[delta]'.sub.C] B (s) where

[[delta].sub.C] = [[??].sub.0] + [[beta].sub.0], [[??].sub.1], [[??].sub.2], ..., [[??].sub.d+K]]' and [[??].sub.0], [[??].sub.1], [[??].sub.2], ..., [[??].sub.d+K] are the mode that gives the maximum likelihood in the Treasury term structure estimation. The linear and quadratic spreads can be similarly expressed.

We first consider the estimation of a constant spread. Since the Treasury securities are backed by the taxing power of the US government, they should be less risky than corporate debt. As a result, the forward rates for corporate debt should be higher than those for Treasury securities, as a premium is required to compensate for the risk. Consequently, we incorporate the knowledge of positive credit spread into our Bayesian model as an informative prior. We assume that the spread is between 0 and a positive constant g prior to estimation, where the prior for a constant spread [[beta].sub.0] is assumed to follow a uniform distribution (0, g):

[[[beta].sub.0] = (1/g) [I.sub.(0,g)] ([[beta].sub.0] [varies] [I.sub.(0,g)] ([[beta].sub.0]. (18)

The posterior distribution for the constant spread is

[[[beta].sub.0] | P, [[sigma].sup.2]] [varies] [[beta].sub.0] [P | [[beta].sub.0], [[sigma].sup.2]] [varies] [I.sub.(0,g)] ([[beta].sub.0] exp (-[n.summation over (i=1)] [([P.sub.i] - [[mu].sub.i].sup.2] / 2[[sigma].sup.2]), (19)

given the error (and thus the likelihood) is normally distributed. This is not of any known distribution and the Metropolis algorithm is applied. Adopting the usual conjugate inverted Gamma distribution IG (a,b) as the prior for the error variance, the corresponding posterior distribution is again an inverted Gamma distribution:

[[[sigma].sup.2] | P, [[beta].sub.0], a, b] ~ IG(n/2 + a, [([n.summation over (i=1)] [([P.sub.i] - [[mu].sub.i]).sup.2] / 2 + 1/b).sup.-1]). (20)

To estimate the linear spread, we adopt a flat prior and add a constraint: 0 < [[beta].sub.0] + [[beta].sub.1] s < g. The prior distribution can be written as: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The constraint forces the credit spread between AT&T bonds and Treasury STRIPS to be positive. For the quadratic spread, the constraint is 0 < [[beta].sub.0] + [[beta].sub.2] [s.sup.2] < g and the prior becomes: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. For both the linear and quadratic spreads, the constraint is: 0 < [[beta].sub.0] + [[beta].sub.1] s + [[beta].sub.2] [s.sup.2] < g and the prior becomes: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. One needs to choose g to complete the prior specification.

This prior information can come from one's subjective knowledge. For example, a trader who has previously traded AT&T bonds before has some idea regarding the spread between STRIPS and AT&T bonds. The prior information can also come from previous empirical studies on the AT&T bonds credit spread. We suppose g to be 0.2 as we believe it is reasonable to assume that the credit spread is between 0 and 0.2 prior to estimation.

We note that the priors we have adopted are equivalent to constraints in a nonlinear optimization problem within a classical framework. Adding these constraints to such a problem may make the asymptotic justification difficult as the average sample size of the individual corporate bonds on a given day is only 4 to 5. On the other hand, the Bayesian approach proposed here allows one to perform valid statistical inference by obtaining the posterior interval, a 100(1 - [alpha])% central interval, from MCMC sample draws.

Choosing the type of spread (constant, linear, or quadratic) is the same as choosing the correct model. The posterior odds can be adopted. For example, to select between a constant spread or a linear spread, the posterior odds ratio

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (21)

can be computed to choose the model. Here [[sigma].sup.2.sub.Const] and [[sigma].sup.2.sub.Lin] are error variances from the constant and linear spread models respectively. The priors on the model parameters are flat and only the likelihood ratio needs to be computed. However, the posterior odds ratio is only an approximation from the MCMC sample draws.

5. CASE STUDY

We present a case study of the end-of-the-month U.S. Treasury STRIPS and AT&T bond prices for each of the 21 months from April 1994 to December 1995. The Treasury term structure is estimated first by Bayesian regression splines and then the AT&T corporate term structure is estimated by adding a positive spread to the former in a hierarchical Bayesian setup.

We focus on one STRIPS data set on December 31, 1995 first. One of the commonly used term structure estimation approaches for large samples such as Treasury bonds in practice is the smoothing spline approach by FNZ (Fisher et al., 1995). The number of knots in FNZ is one-third of the sample size. Thus, FNZ is not strictly a smoothing spline method but more like penalized splines (Ruppert et al., 2003; Jarrow et al., 2004) which do not include as many knots. FNZ smooth the forward rate curve by using a modified GCV (Generalized Cross-Validation) to select the smoothing parameter [lambda]:

GCV([lambda]) = [n.sup.-1] [[summation].sup.n.sub.i=1] [{[P.sub.i] - [[mu].sub.i] ([delta])}.sup.2] / [{1 - [n.sup.-1] [theta]trA ([lambda])}.sup.2], (22)

where trA([lambda]), the trace of the "smoother" A([lambda]), is the effective degrees of freedom and [theta] is a tuning parameter. In standard GCV, [theta] = 1. However, FNZ found that standard GCV with [theta] = 1 did not smooth the forward rate curve enough and [theta] = 2 was a better choice. Indeed, the simulation study in FNZ demonstrates the effectiveness of their smoothing spline method with this modified GCV. However, there is no formal theoretical justification for the nonstandard [theta] = 2 in GCV and FNZ offered no explanation.

JRY (Jarrow et al., 2004), using the classical penalized spline approach, carefully examined different methods of selecting the smoothing parameter, in particular, GCV v.s. EBBS of Ruppert (1997). They found that EBBS gives more smoothing while GCV with [theta] = 1 undersmoothes. GCV with larger value of [theta] mimics EBBS. It is the forward rate f in (2), not the discount function D that needs to be estimated. The standard GCV with [theta] = 1 smoothes the discount function D adequately, but not the forward rate f. The modified GCV with [theta] = 2 in FNZ smoothes the forward rate better and is closer to the EBBS method Jarrow et al. (2004) proposed. EBBS selects the smoothing parameter [lambda] by minimizing the mean square error of the forward rate as a function of [lambda]. BRS smoothes the forward rate curve by adaptively selecting knot locations.

The fitted curves using BRS, FNZ ([theta] = 2 in GCV), and EBBS are displayed in Figure 2. The BRS fit is based on the model averaging estimate (16). On average, 17 knots were selected by the reversible jump MCMC. The empirical forward rates are not observed data but approximations. The true forward rates are not observable. However, the empirical forward rates approximate the true forward rates had them been observable. We found the mean squared error (MSE) or the mean absolute deviation

(MAD = (1/n) [n.summation over (i=1)] [absolute value of [P.sub.i] - [[mu].sub.i])from BRS and FNZ to be virtually the same. The fitted forward rate curves

from both methods are similar but not as smooth as that from EBBS. However, how smooth the forward rate curve should be is still an open question and no one knows the answer. From a limited simulation study (not reported here), BRS performs mostly well or even better than EBBS for i.i.d, errors. However, it is still an open question whether BRS or EBBS can recover the true forward rate curve better and whether the correlated and possibly nonstationary residuals are the reason causing the standard GCV to undersmooth.

After estimating this Treasury term structure of interest rates, we add a spread (see Section 4) to it to estimate the AT&T term structure using all five AT&T bonds traded on December 31, 1995. 1,000 simulation draws from the Metropolis method are used to find the spread estimates and construct the posterior intervals, with another 1,000 draws as the burn-in phase. Our experience is that 1,000 burn-in draws are enough before the Markov chains mix adequately.

First, a constant spread is added to the estimated Treasury term structure and is found to be 0.004 with a 95% posterior interval (0.0036, 0.0043). The proposal distribution is normal and only the samples that satisfy the constraint [[beta].sub.0] > 0 are used. A linear spread, [[beta].sub.0] + [[beta].sub.1] s, is then considered. The two parameter estimates are found to be 0.004 and 3.25x[10.sup.-6]. The 95% posterior intervals for both parameters are (0.0036, 0.0043) and (-1.7x[10.sup.-6], 8.6x[10.sup.-6]). The posterior odds ratio for a constant spread compared to a linear spread is found to be 1.0084 and the simpler constant spread model is favored. At maturity time 0, the risk of default by AT&T should be negligible and the credit spread should be due to liquidity risk, not credit risk. The name "credit spread" is a misnomer since it can be due to other factors in addition to credit differences (see Marshall, 2000). AT&T bonds are not as liquid as Treasury bonds and they may not be sold immediately if necessity arises. A premium in the interest rate is needed to compensate for this liquidity risk and find an immediate buyer. The intercept [[beta].sub.0] of the credit spread can be interpreted as the liquidity risk and it does exist based on the posterior interval and odds ratio. However, the posterior interval and odds ratio merely mean that a constant spread is more likely and they do not mean the spread is indeed constant. A sample size of only 5 may not be sufficient to detect a difference from a constant spread. A constant spread model can serve as a good first approximation for AT&T and many other corporate bonds with a minimal number of bonds available. For individual corporate bonds with larger sample size, non-constant spreads may very well be found.

We also consider a quadratic spread, [[beta].sub.0] + [[beta].sub.1] [s.sup.2], and find the posterior odds ratio to be 16.32. The quadratic spread model is not favored. The 95% posterior intervals for both parameters are (-1.2x[10.sup.-3], 6.6x[10.sup.-5]) and (-1.4x[10.sup.-7], 4.9x[10.sup.-5]). The case of including both the linear and the quadratic spread is not considered for this AT&T data set. It is not desirable to estimate three parameters with only five data points. However, for other individual corporate bonds with a larger sample size, including both the linear and the quadratic spreads may be considered.

We have only focused on the term structure of the December 1995 data. However, modeling the evolution of the term structures of both Treasury and individual corporate bonds is an important problem in pricing interest rate derivatives, among other applications in finance (Jarrow, 2002). To study this evolution, the method described in Section 3 is applied 21 times independently to the 21 monthly U.S. Treasury STRIPS data sets from April 1994 to December 1995. Figure 3 displays the 21 fitted forward rate curves. The STRIPS data are sparse after 25 years to maturity and only the forward rate curves up to 25 years to maturity are plotted. We do not recommend using the fitted values after 25 years to maturity. By graphing the 21 fitted curves in one figure, we show both the evolution of the end-of-month forward rates from April 30, 1994 to December 31, 1995 and the rates based on years to maturity from 0 to 25 years. If maturity is fixed, the forward rate is a function of time.

[FIGURE 3 OMITTED]

Given the 21 estimated STRIPS forward rate curves, we apply (21 times independently) the Bayesian method developed in Section 4 to the 21 monthly AT&T bond data sets from April 1994 to December 1995. A constant spread is favored according to the posterior interval and odds ratios. Figure 4 plots the fitted forward rates for U.S. STRIPS and AT&T bonds over the 21-month period of April 1994 to December 1995. The upper sheet contains 21 AT&T forward rate curves with constant spreads and the bottom sheet contains the 21 STRIPS forward rate curves. The range of years to maturity for these 21 monthly data sets is from 0.9 to 12 years and is used in the figure. Extrapolation beyond this range is not desirable. The figure shows how the forward rates evolve from month to month as well as over time to maturity. On the other hand, if the maturity is fixed, the forward rate as a function of time is observed and the curves are rough. This is not surprising in that interest rates move randomly and abruptly. Thus, smoothing both maturity and time with a bivariate function is not desired.

[FIGURE 4 OMITTED]

6. CONCLUSION

We have presented a case study of Treasury and corporate term structure estimation with Bayesian regression splines. We estimate the Treasury term structure, adaptively choosing location of the knots. Due to the small sample size of corporate bonds, credit spreads are added to the estimated Treasury term structure to obtain the corporate term structure. The fact that the positive credit spreads between the default-free Treasury securities and the risky corporate bonds is too important to ignore. The Bayesian approach naturally takes this into account as a prior distribution and also gives valid statistical inference. On the other hand, achieving these objectives within the classical framework may be difficult to justify.

REFERENCES:

Adams, K. J. and Van Deventer, D. R., "Fitting yield curves and forward rate curves with maximum smoothness", Journal of Fixed Income, 1994, 52-56.

Bielecki, T. and Rutkowski, M., Credit Risk: Modelling, and Hedging, Springer, New York, 2002.

Chambers, D. R., Carleton, W. T., and Waldman, D. W., "A new approach to estimation of the term structure of interest rates", Journal of Financial and Quantitative Analysis, Vol. XIX (3), 233-252.

de Boor, C., 2001, A Practical Guide to Splines. New York: Springer.

Denison, D.G.T., Mallick, B. K., and Smith, A. F. M., "Automatic Bayesian curve fitting", Journal of the Royal Statistical Society, B, Vol. LX, 330-350.

DiMatteo, I., Genovese, C. R., and Kass, R. E., "Bayesian curve fitting with free-knot splines", Biometrika, Vol. LXXXVIII (4), 1055-1071.

Duffie, D. and Singleton, K., "Modeling Term Structures of Defaultable Bonds", Review Financial Studies, Vol. XII (4), 1999, 197-226.

FDIC (Federal Deposit Insurance Corporation), "Options Paper", 2000, August issue.

Fisher, M., Nychka, D. and Zervos, D., "Fitting the term structure of interest rates with smoothing splines", FEDS 95-1, 1995, Federal Reserve Board, Washington DC. (available at www.markfisher.net).

Green, P. J., "Reversible jump Markov chain Monte Carlo computation and Bayesian model determination", Biometrika, Vol. LXXXII (4), 1995, 711-732.

Hansen, M. H. and Kooperberg, C., "Spline Adaptation in Extended Linear Models", Statistical Science, Vol. XVII (1), 2002, 2-51.

Jarrow R., Modeling Fixed-Income Securities and Interest Rate Options, 2nd ed., Stanford University Press, Palo Alto, 2002.

Jarrow, R., Ruppert, D., and Yu, Y., "Estimating the term structure of corporate debt with a semiparametric penalized spline model", Journal of the American Statistical Association Vol. XCIX (465), 2004, 57-66.

Jarrow, R. and Turnbull, S., "Pricing derivatives on financial securities subject to credit risk", Journal of Finance, Vol. L, 1995, 53-85.

Li, M. and Yu, Y., "Estimating the Interest Rate Term Structures of Treasury and Corporate Debt with Bayesian Penalized Splines", Journal of Data Science, 2005, to appear.

Linton, O., Mammen, E., Nielsen, J. P., and Tanggaard, C., "Yield curve estimation by kernel smoothing methods", Journal of Econometrics, Vol. CV (1), 2001, 185-223.

Marshall, J. F., Dictionary of Financial End, Wiley, New York, 2000.

McCulloch, J. H., "Measuring the term structure of interest rates", Journal of Business, Vol. XIV (1), 1971, 19-31.

McCulloch, J. H., "The tax-adjusted yield curve", Journal of Finance, Vol. XXX (3), 1975, 811-830.

Nelson, C. R. and Siegel, A. F., "Parsimonious modeling of yield curves", Journal of Business Vol. LX (4), 1987, 473-489.

Ruppert, D., Wand, M. P., and Carroll, R. J., Semiparametric Regression, Cambridge University Press, Cambridge, 2003.

Schwarz G., "Estimating the Dimension of a Model", Annals of Statistics, Vol. VI (2), 1978, 461-464.

Shea, G., "Pitfalls in smoothing interest rate term structure data: equilibrium models and spline approximations", Journal of Financial and Quantitative Analysis, Vol. XIX (3), 1984, 253-269.

Vasicek, O. A., and Fong, H. G., "Term structure modeling using exponential splines", Journal of Finance, Vol. XXXVII (2), 1982, 339-356.

Wahba, G., Spline Models for Observational Data, SIAM, Philadelphia, 1990.

Warga, A., Fixed Income Securities Database. Fixed Income Research Program, University of Houston, 1995. (Request information at [e-mail]).

AUTHOR PROFILES:

Dr. Min Li earned his Ph.D. at the University of Cincinnati in 2002. Currently he is an assistant professor of quantitative methods at California State University, Sacramento.

Dr. Yan Yu earned his Ph.D. at Cornell University in 2000. Currently he is an assistant professor of quantitative analysis and operations management at the University of Cincinnati.

TABLE 1. AT&T BONDS ON DECEMBER 31, 1995.
COUPON LISTED IS THE SEMI-ANNUAL DOLLAR COUPON PAYMENT

Current Date     Issue Date    Maturity Date    Market Price    Coupon

19951231          19920114       20020115         109.458       7.125
19951231          19940324       20040401         106.284       6.750
19951231          19940601       20060601         111.436       7.500
19951231          19950228       20070301         115.509       7.750
19951231          19950512       20050515         107.659       7.000