Zero Coupon Yield Curve Estimation with the Package Termstrc

Since zero-coupon rates are rarely directly observable, they have to be estimated from market data. In this paper we review several widely-used parametric term structure estimation methods. We propose a weighted constrained optimization procedure with analytical gradients and a globally optimal start parameter search algorithm. Moreover, we introduce the R package termstrc, which offers a wide range of functions for term structure estimation based on static and dynamic coupon bond and yield data sets. It provides extensive summary statistics and plots to compare the results of the different estimation methods. We illustrate the application of the package through practical examples using market data from European government bonds and yields.


Introduction
The term structure of interest rates defines the relationship between the yield of a fixed income investment and the time to maturity of its cash flows. Accordingly, the zero-coupon yield curve provides the relationship for investments with only one payment at maturity. It serves as the basis for the valuation of other fixed income instruments and as an input for various models, e.g., for risk management, monetary policy, derivatives pricing. Although zero-coupon prices can be directly used to construct the term structure, the lack of market liquidity and the limited available maturity spectrum necessitates the estimation based on observed coupon bond prices.
In this paper we briefly review important term structure estimation methods. We focus on the cubic splines approach of McCulloch (1971McCulloch ( , 1975 and the Nelson and Siegel (1987) method with extensions by Svensson (1994), Diebold and Li (2006) and De Pooter (2007). According to a survey of Bank for International Settlements (2005) they are widely applied at central banks. We give a detailed discussion of the model estimation, and propose a nonconvex weighted global optimization procedure with linear constraints, analytical gradients and an appropriate start parameter search algorithm. Moreover, we introduce the package termstrc which is written in the R system for statistical computing (R Development Core Team 2010). It is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/ package=termstrc and from the R-Forge development platform (Theußl and Zeileis 2009) at http://R-Forge.R-project.org/projects/termstrc/. The package provides a wide range of functions for estimating and analyzing of the term structure of interest rates. All methods, except the cubic splines approach, can be applied statically and iteratively on dynamic coupon bond price and zero yield data sets.
In the following two subsections we review the relevant literature and the currently available term structure estimation software. In Section 2 we introduce the notation and discuss different models and estimation procedures. The structure of the package termstrc (i.e., the available functions, S3 classes and methods) is presented in Section 3. Several examples for the different models based on dynamic/static bond/yield data sets are presented in Section 4. Section 5 concludes the paper.

Literature review
The first term structure model goes back to Vasicek (1977, Vasicek in the following), in which he assumed that the short rate followed a univariate continuous Markovian diffusion process and all riskless zero-coupon bond prices were functions of the short rate, the current time and the maturity date of the bond. Whereas Vasicek used an Uhlenbeck and Ornstein (1930) process for the short rate, Cox, Ingersoll, and Ross (1985, CIR in the following) proposed a one-dimensional square root diffusion process, which ensured positive short rates. CIR derived their model by applying an equilibrium from utility-maximizing identical individuals with a log-utility function. Nawalkha, Beliaeva, and Soto (2007) classified term structure models, which assume a time-homogeneous short rate process and require and explicit market price of risk specification as fundamental models. In addition to the Vasicek and CIR model, the multifactor affine term structure models (ATSMs) of Duffie and Kan (1996) and Dai and Singleton (2000) also belong to the class of fundamental models. To estimate these models a time series of a cross section of fixed income security prices is required. The estimation of multifactor ATSMs can be challenging, especially with nonlinear market price of risk specifications. Moreover, the model implied prices may or may not converge to the observed market prices. However, for the purpose of pricing fixed income derivatives, models which exactly fit the term structure are required. Based on Vasicek and CIR, models were proposed, which assume an endogenous term structure. They price the initially observed zero-bond prices without errors by allowing time inhomogeneity in the parameters of the stochastic differential equation for the short rate. Examples of so-called endogenous or no-arbitrage models are Ho and Lee (1986), Hull and White (1990) and Black and Karasinski (1991).
Another class of models, which has not been deduced from equilibrium conditions or/and no-arbitrage conditions, takes a more empirical approach. In general these models assume a parametric form of the spot rate, forward rate or discount function. The unknown parameters are then estimated by minimizing the error between theoretical and observed prices of a cross section of coupon bonds at a certain point in time. The method of Fama and Bliss (1987) iteratively extracts the forward rates by extending the discount function at each step of the calculation, i.e., the forward rates are computed to price bonds with increasing maturity given the discount function fitted to the previous bonds with shorter maturity. The obtained spot or forward rates are referred to as unsmoothed Fama/Bliss yields. The discount function is piecewise linear with the number of parameters equal to the number of included securities. The procedure works well only if all cash flows have the same maturity intervals (see, e.g., Hagan and West 2006). McCulloch (1971McCulloch ( , 1975 proposed using splines to fit the discount function of the segmented term structure. Several different types of splines have been suggested as well as the use of penalty functions, e.g., the penalized spline model of Fisher, Nychka, and Zervos (1995). Nelson and Siegel (1987, NS in the following) proposed a four parametric parsimonious function from the family of general exponential-polynomial functions for the forward curve, which was later extended by Svensson (1994, SV in the following). Both specifications make it possible to reproduce a wide range of possible term structure shapes. The estimation of the parameters typically involves a nonconvex optimization procedure and is based on coupon bond prices. However, it also possible to use interpolated zero or Fama/Bliss yields. Several papers compare the performance of empirical term structure estimation methods, (see, e.g., Bliss 1997;Bolder and Streliski 1999;Ioannides 2003). The models of NS and SV are not formulated in a dynamic framework and the first model is not consistent with arbitrage-free pricing theory (see, Björk and Christensen 1999;Filipovic 1999). The first disadvantage was addressed by Diebold and Li (2006), who iteratively applied a simplified NS model (DL in the following) to a dynamic zero yield data set. Subsequently, they estimated common time series models for the parameter time series in order to forecast future zero-coupon yield curves. Rudebusch (2007, 2009) tried to correct for the second drawback. They developed two ATSMs which belonged to the class of Dai and Singleton (2000), where the spot rate function is similar to the version of NS and SV.

Software review
The package sde developed by Iacus (2009) provides functions for the parametric estimation and numerical methods for the simulation of univariate stochastic processes. The afore mentioned Vasicek and CIR term structure model assumes that the instantaneous spot rate follows an one-dimensional diffusion process. The former model requires the short rate under the real measure to follow an Uhlenbeck and Ornstein (1930) process with constant coefficients, whereas the latter postulates a square-root diffusion process. The package makes it possible to simulate the short rate based on the underlying differential equation and to estimate the parameters. However, it is not possible to estimate the term structure based on coupon bonds or to price fixed income securities with the offered functions. Moreover, only univariate processes are covered. The yuima project is currently developing a more general framework for the simulation and estimation of multivariate stochastic processes. For more details concerning sde we refer to Iacus (2009) and concerning the yuima project to Iacus (2010).
The existing software, which allows the estimation of the mentioned empirical term structure models can be distinguished according to the type of data used. The R package YieldCurve developed by Guirreri (2010) offers functions for the estimation based on zero yields. The available models include NS, SV and DL. The functions require a user-specified range of certain parameters of the spot rate function. All other parameters are subsequently obtained by a linear grid search. Unfortunately, nececessary constraints on the parameters were not considered. The package fBonds of Wuertz (2009) requires forward rates for the estimation and the available methods are NS and SV. In contrast to the YieldCurve package, fBonds obtains globally optimal start parameters from a grid search and afterwards an unconstrained estimation procedure leads to the optimal parameters. Neither optimization procedures allow the user to specify the step size of the grid. The previous two packages consist of a limited number of highly similar functions which do not allow the estimation based on observed prices of coupon bonds and ignore constraints on the parameters. Moreover, they provide no standardized S3 methods. In contrast, the library QuantLib and its R interface RQuantLib developed by Eddelbuettel and Nguyen (2010) offers several methods to estimate the term structure based on coupon bonds (i.e., simple polynomial fitting, various types of splines and the NS model). Additionally, popular one-and two-factor short rate models are covered. For the estimation of the NS model RQuantLib offers a weighted optimization algorithm. The start parameters of the optimization are selected automatically, however the underlying logic is not documented. The estimation is currently limited to fixed-coupon bonds. Moreover, the ability in RQuantLib to analyse and compare the results of the estimation appears to be limited, e.g., it is impossible to obtain the optimal parameters of the NS spot rate function. The commercial software MATLAB provides in the Fixed-Income toolbox functions (The MathWorks, Inc. 2010) which allow the term structure estimation based on observed fixed income prices and yields. The available models include NS, SV and smoothed splines. The estimation is carried out by applying a constrained weighted optimization procedure with box constraints. However, a start parameter search algorithm and linear constraints are not provided.

Zero-coupon yield curve estimation
Before we discuss the problem of model specific zero-coupon yield curve estimation, we introduce the definitions of a few basic terms used in the fixed income literature and the associated notation. Two very fundamental fixed income securities are discount or zero-coupon bonds and fixed-coupon bonds. The first type of bond is a fixed income investment with only one payment at maturity, whereas the second type provides fixed periodical coupons and a redemption payment at maturity. For a coupon bond j we denote the occurring cash flows, i.e., the coupons and the redemption payment, with c ij and the associated maturities with m ij . For a group of k bonds we summarize the cash flows and maturities in the matrices M = {m ij } and C = {c ij } with t rows and k columns, with i = 1, . . . , t and j = 1, . . . , k. The number of rows t is determined by the number of cash flows of the j-th bond with the longest maturity. Dates after the maturity of each bond j are completed with zeros until the maturity date of the bond with the longest maturity. One element m ij of M refers, therefore, to the time to occurrence (in years) of the i-th cash flow of the j-th bond. The maturities of the last cash flows m j of the bonds are collected in the row vector m = {m j } of dimension k .
The spot rate, zero-coupon rate or zero yield s(m ij ) denotes the interest paid on a discount bond for a certain maturity m ij . With continuous compounding and in the absence of arbitrage opportunities the fair price of a discount bond p j paying one Euro at maturity m j is given by The spot curve (or zero-coupon yield curve) shows the spot rates for different maturities. The forward rate f (r, t) is the interest contracted now to be paid for a future investment between the time r and t (r < t), denoted by m r and m t . The forward rate as a function of maturity is the forward curve. Assuming continuous compounding, we observe the following relationship between spot rates and forward rates e s(mr)mr e f (r,t)(mt−mr) = e s(mt)mt .
We can solve this for the forward rate The instantaneous forward rate describes the return for an infinitesimal investment period after r.
Another interpretation is the marginal increase in the total return from a marginal increase in the length of the investment period. Thus, the spot rate can be seen as the average of the instantaneous forward rates where m is the the time to maturity.
At the market the following information is typically available for a bond j: the clean price p c j (given as percentage of the nominal value), the cash flows c ij (coupons and redemption payment) and their maturity dates m ij . When an investor buys a bond j, she obtains the right to receive all its future cash flows. If the purchase occurs between two coupon dates, the seller must be compensated for the fraction of the next coupon, the so-called accrued interest a j . Depending on the market specific day-count convention, e.g., 30/360, Actual/360, a basic form for the accrued interest of bond j is a j = number of days since last coupon payment number of days in current coupon period × c ij .
Therefore, the purchase price also named dirty price p j is the sum of the quoted market price (clean price) p c j and the accrued interest a j . We again summarize the dirty and clean prices and accrued interest of a group of bonds into the k-dimensional row vectors p, p c , a, where p = p c + a = {p j = p c j a j }. By considering the accrued interest, similar to (1), the bond pricing equation for a bond j under continuous compounding is the present value of all cash flows where δ(m ij ) is the discount factor for a certain maturity m ij . Applying the concept of spot rates the discount factors can be expressed as δ(m ij ) = e −s(m ij )m ij . The t × k matrix D = {δ(m ij )} consists of all discount factors for a group of bonds and the bond pricing equation in matrix notation is therefore given by Throughout the paper "·" denotes an element-wise multiplication, () is the transpose of a vector or matrix and ι is a column vector filled with ones. Another possibility to express the bond pricing equation, is to use the internal rate of return of the cash flows. The so-called yield to maturity (YTM) of a bond j is the solution for y j in the following equation We summarize the yield to maturities of a group of bonds in the k-dimensional row vector y = {y j }. The YTM for a discount bond is equal to the spot rate. However, this is not valid for coupon bonds. Plotting just the yield to maturity for coupon bonds with different maturities does not result in a yield curve which can be used to discount cash flows or price any other fixed income security, except the bond from which it was calculated. Therefore, estimating the term structure of interest rates from a set of coupon bonds can not be seen as a simple curve-fitting of the YTMs.
To quantify the sensitivity of a bond's price against changes in the interest rate, one needs to account for the fact that coupons are paid during the lifetime of a coupon bond. A standard measure of risk is the (Macaulay) duration which computes the average maturity of a bond using the present values of its cash flows as weights. The k-dimensional vector d consists of the individual durations and is calculated by where the discount matrix D contains the discount factors calculated with the yield to maturity of each bond, i.e., δ(m ij ) = e −y j m ij (the division is applied elementwise).
The objective of term structure estimation is to extract the spot, forward or discount curve out of a set of coupon bonds. The simplest methods to derive the zero-coupon yield curve are direct methods, e.g., bootstrapping. Such techniques price the bonds without error, which leads to overfitting. The direct methods require bonds with nearly identical cash flow dates and lead to nonsmooth spot, forward or discount curves. Therefore, to overcome the mentioned disadvantages indirect methods have been proposed. They postulate a parametric parsimonious form of the spot, forward or discount function. In the following we denote the vector which contains the parameters of the indirect methods by β. Because bond prices are observed with idiosyncratic errors we define the theoretical bond prices aŝ Thus, the market prices of set of coupon bonds p can be expressed as the sum of the theoretical bond pricesp plus an idiosyncratic error vector . The estimated parameter vectorβ is then obtained by a possible weighted nonconvex optimization procedure. So far we have only mentioned the estimation of the zero-coupon yield curve based on a set of coupon bonds. However, another application is based on zero yields, which may be obtained from a bootstrap procedure. The optimization then aims at minimizing the error between the observed y and the theoretical zero yieldsŷ.
The following two subsections focus on two popular indirect estimation methods, i.e., the method of NS and SV which propose a specific parametric forward and spot rate function and the cubic splines method which models the discount function.

Nelson/Siegel and Svensson method
Spot curve parameterizations Nelson and Siegel (1987) propose a parsimonious function for modelling the instantaneous forward rate as a solution to a second-order differential equation for the case of equal roots with the parameter vector β = (β 0 , β 1 , β 2 , τ 1 ). As in (2), the spot rate is the average of the instantaneous forward rates Instead of one specific maturity m ij the spot rate function can also be written in terms of a maturity vector m, i.e., s(m, β). Consequently, all calculations in (3) have to be applied element by element. The specification of NS can produce a wide range of possible curve shapes, including monotonic, humped, U -shapes or S-shapes. Svensson (1994) adds the term with two new parameters β 3 and τ 2 to increase the flexibility. The parameter vector is then given by β = (β 0 , β 1 , β 2 , τ 1 , β 3 , τ 2 ). This specification of the spot rate allows for a second hump in the curve. Based on the proposed spot rate function of NS and SV several extensions have been developed. In order to simplify the estimation procedure DL suggest to reduce the parameter vector to β = (β 0 , β 1 , β 2 ) by fixing τ 1 = 1 λ to a prespecified value. A potential multicollinearity problem in the SV model arises if the decay parameters τ 1 and τ 2 have similar values. As a result only the sum of β 2 and β 3 can be estimated efficiently. To circumvent the multicollinearity problem De Pooter (2007) proposes similar to Björk and Christensen (1999) to replace the last term in (4) . In the following we refer to the SV model with the replaced term as adjusted SV. Because the named models are nested in the formulation of the (adjusted) SV spot rate function the following parameter interpretations also hold: β 0 is the asymptotic value of the spot rate function lim m ij →∞ s(m ij , β), which can be seen as the long-term interest rate. Due to the assumption of positive interest rates it is required that β 0 > 0. β 1 determines the rate of convergence with which the spot rate function approaches its long-term trend. The slope will be negative if β 1 > 0 and vice versa.
The instantaneous short rate is given by lim m ij →0 s(m ij , β) = β 0 + β 1 , which is constrained to be greater zero.
β 2 determines the size and the form of the hump. β 2 > 0 results in a hump at τ 1 , whereas β 2 < 0 produces a U -shape.
τ 1 specifies the location of the first hump or the U -shape on the curve. β 3 , analogously to β 2 , determines the size and form of the second hump.
τ 2 specifies the position of the second hump.
For the NS (DL), SV (adjusted SV) spot rate function defined in (3) and (4) the elements of the discount matrix D can be calculated as follows:

Weighted least squares and globally optimal parameter estimation
The unknown parameter vector of the spot rate functions β can be obtained by minimizing the error between the observed bond prices p and the theoretical bond pricesp. However, when minimizing the unweighted price errors, bonds with a longer maturity obtain a higher weighting, due to higher degree of price sensitivity, which leads to a less accurate fit at the short end. Therefore, a weighting of the price errors has to be introduced to solve this problem or at least to reduce the degree of heteroscedasticity. As Martellini, Priaulet, and Priaulet (2003) point out, the choice of the weights is a crucial part of the optimization procedure. We use a specification which is based on the inverse of the duration and was proposed by Bliss (1997). The weight ω j for a bond j is given by Alternative weighting schemes also make use of the duration (see, e.g., Vasicek and Fong 1982) or the bid-ask spread (see, e.g., Nawalkha, Soto, and Beliaeva 2005). We summarize the individual weights for k bonds into the k-dimensional row vector ω. The consideration of the weights leads to the first objective function, which minimizes the sum of the weighted squared coupon bond price errors The second minimizes the sum of the squared zero yield errors. The objective function F (β) is given by The entire optimization problem for the SV spot rate function considering all mentioned constraints can be written as The estimation problems using the NS, DL and adjusted SV spot rate function are nested in the above formulation. Therefore, we refrain from reporting upon them. We derive the analytical gradients ∇F (β) for both objective functions and report them in Appendix A.
The optimization problem stated in (6)-(8) is nonconvex and may has multiple local optima, which increases the dependence of the numeric solution on the starting values considerably. Choosing the start parameters arbitrarily and possibly not reaching a global optimum can be avoided by using the following reformulation of the problem based on an idea of Werner and Ferenczi (2006). The parameter vector β is split in a local component β * and a global component τ . Then the original optimization problem in (6)-(8) is equivalent to: where The inner problem (12)-(14) is a (weighted) least squares problem with four variables under linear constraints, which is convex and therefore easy to solve numerically. In contrast, the outer problem (9)-(11) is nonconvex in two variables. Werner and Ferenczi (2006) formulate the optimization problem with the usual parameter constraints given in Svensson (1994), i.e., with only box-constraints for the τ -parameters. Under these constraints the nonconvex outer problem is solved with a grid search method. Werner and Ferenczi (2006) use various global search algorithms (e.g., sparse grids, HCP algorithm) for this purpose. We apply a full grid search in which we are able to reduce the size by at least the half. The reduction is a consequence of the imposed minimum distance constraint ∆ τ between τ 1 and τ 2 . The additional linear constraint for the τ -parameters was proposed by De Pooter (2007) and avoids identification problems that could arise from the similar factor loading structure associated with the parameters β 2 and β 3 . In a static setting, the only consequence of not imposing this additional constraint can be rather extreme parameter estimates, where the total factor contributions virtually cancel each other out. However, this does not influence the accuracy of the fit. De Pooter (2007) notes that the real advantage of the additional constraint are smoother parameter time series when performing a two-step dynamic estimation. Furthermore, the interpretation of τ 1 and τ 2 as the position of the humps is ensured, because they are not able to switch their positions. We can confirm these results and have included the necessary constraints in the code of the optimization algorithm. To sum up, the local optimization problem in six parameters defined in (6)-(8) is solved with the globally optimal starting point from the grid search procedure. An advantage of our method is that the user does not have to specify start paramaters in advance. We apply the above algorithm, with a reduced parameter space (β ∈ R 4 , β * ∈ R 3 , τ ∈ R 1 ), also for the Nelson/Siegel model.

Cubic splines
The cubic splines based term structure estimation method divides the term structure into segments using a series of so-called knot points. Cubic polynomial functions are then used to fit the term structure over these segments. The polynomial functions ensure the continuity and smoothness of the discount function within each interval. McCulloch (1971McCulloch ( , 1975 use the following definition of the discount factors: Where g l (m ij ) (l = 1, . . . , n) defines a set of piecewise cubic functions, the so-called basis functions, which satisfy g l (0) = 0. The functions have to be twice-differentiable at each knot point to ensure a smooth and continuous curve around the points. The unknown parameter vector β can be estimated with ordinary least squares (OLS).

Knot point selection
McCulloch (1975) defines a n-parameter spline with n − 1 knot points q l . We sort the cash flow matrix C and the maturity matrix M such that the k bonds are arranged in ascending order by their maturity dates m. The following specification places an approximately equal number of bonds between adjacent knots. For 1 < l < n − 1 the knot points are defined as and θ = (l−1)k n−2 − h. The first knot point q 1 = 0 and the last one is equal to the maximum maturity, i.e., q n−1 = m k . McCulloch (1971) proposes to set the number of basis functions n to the integer nearest to the square root of the number of observed bonds k, i.e., n = √ k + 0.5 .

Basis functions for cubic splines
For l < n the set of basis functions is defined by For l = 1 we set q l−1 = q l = 0. When l = n, the basis function is given by g l (m ij ) = m ij . For a set of bonds we summarize the basis functions for all cash flows m ij in the t × k matrix G l = {g l ij (m ij )}.

Regression fitting of the discount function
For a set of k coupon bonds the dirty price vector can be expressed as the sum of the discounted cash flows, i.e., the theoretical prices plus an idiosyncratic error where the discount factor matrix is defined as the weighted sum of the l = 1 . . . n basis functions D = 1 + β 1 G 1 + · · · + β n G n .
The difference between the observed price vector p and the sum of the cash flows is summarized by a t-dimensional column vector z, which is given by p − ι C . The unknown parameters are denoted by the t × 1 vector β. The multivariate linear regression equation is therefore given by z = Xβ + . The parameter estimates are obtained by simple OLS, i.e., the n × 1 vectorβ = X X −1 X z.
We can use the resulting parameters to calculate the discount function in (15) for any given maturity m ij between the first and the last knot point. The discount factors can then be converted into the spot rates by the following relationship Confidence intervals for the discount function McCulloch (1975) plots error bands one standard error above and below the best estimate. However, it is also possible to derive a confidence interval for the estimated discount function. Under the assumption of n.i.i.d. disturbances with variance σ 2 , the ordinary least squares coefficient estimator is normally distributed with mean β and variance-covariance matrix σ 2 X X −1 . If the disturbances exhibit heteroscedasticity or/and autocorrelation the estimated covariance matrix σ 2 X X −1 will be inappropriate. The appropriate matrix is then given by Ψ = σ 2 X X −1 (X ΩX)(XX) −1 , where Ω is the weights matrix. Therefore, autocorrelation and heteroscedasticity consistent estimators for Ψ have been proposed. Within the package and for the calculations in this paper we use an estimator developed by Andrews (1991) and provided by the R package sandwich of Zeileis (2004Zeileis ( , 2006, . According to Greene (2002), the confidence interval for a linear combination of OLS coefficients follows a normal distribution. Therefore, the discount function in (15) is normally distributed with mean µ d = 1+g(m ij ) β and variance-covariance matrix σ 2 d = g(m ij ) Ψg(m ij ), where g(m ij ) is a n-dimensional column vector with consists of the elements g l (m ij ) (l = 1, . . . , n). The 1 − α confidence interval for the mean of the discount function can now be constructed as follows where s d is the estimate for σ d , t α/2 the appropriate critical value from the t distribution with k − n degrees of freedom and 1 − α the desired level of confidence.

Software implementation in R
In this section, we describe the structure of the implementation of the package termstrc in the R system for statistical computing. The matrix-oriented notation introduced in Section 2 can efficiently be realized with R base functions. However, a few functions depend on other packages, which we appropriately reference. Table 1 lists the basic structure of the package. It contains of three data classes and two associated estimation methods. Subsequently, we explain in detail the structure of our data sets, the core functions and describe the available standard and generic S3 classes and methods provided by the package.

Data structure
The package termstrc makes it possible to estimate the term structure based on observed data of coupon bonds and yields. Accordingly, we have defined two basic types of data sets and corresponding classes, i.e., the class "couponbonds" for a set of coupon bonds and the class "zeroyields" for zero yields . An object of the class "couponbonds" has to be organized in a list form, which contains sublists for each country or group of bonds. The structure of the sublist is illustrated in Table 2.
The first part consists of the general specifications of the bonds, its market prices and accrued interest. The sublist CASHFLOWS contains all cash flows for the available bonds sorted by their ISIN and their maturity date. The list structure allows us to add additional data easily, e.g., bid and ask prices, volumes. Due to the different day-count-conventions, business holidays and market conventions, no functions for the calculation of accrued interests, maturity dates and cash flows are provided. The data structure provides itself to be convenient when building the cash flow matrices and the maturity matrices defined in Section 2. Moreover, the creation of dynamic data sets can be easily achieved by merging static "couponbonds" data sets. For the resulting object we have invented the class "dyncouponbonds". Objects of the class "zeroyields" are structured in a similar way. However, only sublists with the yields and the associated maturities are required. Due to the simple data structure the package provides no separate class for static zero yield data. Moreover, because of the different data formats of the data providers no general class constructor functions, except for the class "zeroyields", are included. For the classes "couponbonds" and "dyncouponbonds" the package offers a rm_bond() method, which allows us to remove a specified vector of bonds and their associated Data class Method Returned object Method for class "couponbonds" estim_cs() "termstrc_cs" print(), plot(), estim_nss() "termstrc_nss" summary() "dyncouponbonds" estim_nss() "dyntermstrc_nss" print(), plot(), summary(), param() "zeroyields" estim_nss() "dyntermstrc_yields" print(), plot(), summary(), param()  Table 2: Structure of a sublist of an object of the class "couponbonds". data from the data set.

Auxiliary functions
The term structure estimation functions are based on several functions which perform typical fixed income mathematics operations. The majority of them can be applied independently. The most important functions are create_cashflows_matrix() and create_maturities_matrix() for the creation of the cash flows and maturities matrices. Based on the obtained matrices the theoretical prices and yield to maturities of the bonds can be calculated with bond_prices() and bond_yields(). The implemented spot rate and forward rate functions for the different methods allow the calculation of the required interest rate type based on the parameter and maturity vector. The function duration() calculates basic interest rate risk measures, i.e., the Macauly and modified duration.

Core functions
The package offers functions for the estimation of the term structure based on NS, SV, DL, adjusted SV implemented in the method estim_nss() and cubic splines implemented in estim_cs(). The first method is available for objects of the class "couponbonds", "zeroyields" and "dyncouponbonds", whereas the latter one only for "couponbonds" objects. All methods require the appropriate data set as an first input argument. The methods available for coupon bond data sets allow the restriction of the estimation up to a certain maturity range and group of bonds. Further inputs of the estim_nss() method enable the user to provide constraints for the start parameter grid search, to set the options of the used optimizer appropriately and to define the fixed parameter λ = 1 τ 1 of the DL spot rate function. The estim_nss() method for "couponbonds" objects gives the user the option to provide start parameters for the optimization. Without them the grid search algorithm is automatically applied. For the iterative estimation of the term structure the package includes the estim_nss() methods for "dyncouponbonds" and "zeroyields" objects. The first calls at every estimation step the static estim_nss() estimation method for "couponbonds" objects, because an object of the class "dyncouponbonds" consists of subobjects of the class "couponbonds". For the iterative estimation a global optimization at every or the first point in time is possible. For the second alternative the results of the previous estimation are used as start parameters for the current one. All the mentioned functions follow the same internal structure, i.e., data preprocessing, start parameter grid search, optimization and postprocessing of the results. In order to improve the speed of the term structure estimation based on coupon bonds and the associated start parameter grid search we outsourced performancecritical functions, i.e., the objective functions and gradients for the NS and SV methods, into C++ by applying the package Rcpp of Eddelbuettel and Francois (2010).

Exploring the estimation results
The methods for objects of the class "couponbonds" and "dyncouponbonds" return objects of the class "termstrc_nss" or "termstrc_cs" and "dyntermstrc_nss", whereas the method applied on the class "zeroyields" returns an object of the class "dyntermstrc_yields". All the obtained objects from a certain term structure estimation method contain the estimation results, pricing errors, input parameters, pre and postprocessed data, results of the start parameter grid search and method specific information in a list format. The sublists follow the structure of the input data sets, i.e., they are classified according to countries or an other group argument. For all the mentioned result classes we provide appropriate print(), plot() and summary() methods. The summary() method gives goodness of fit measures, i.e., the RMSE and AABSE for the pricing and yield errors and a convergence information from the solver is printed. For static results objects the print() method shows the estimated parameters and method specific information, whereas for dynamic result objects the print() method offers aggregated summary statistics on the estimated parameters. The parameters of the cubic splines approach are estimated with OLS by using the lm() function. Consequently, the print() method for "termstrc_cs" also provides statistics for the OLS estimation by applying the summary() method for "lm" objects or in case of rse = TRUE robust standard errors are used for the calculation of the statistics. The method applies the function coeftest() from the package lmtest of Zeileis and Hothorn (2002). Additionally, we implemented a param() method for objects of the class "dyntermstrc_yields" and "dyntermstrc_nss". The method allows to conveniently extract the estimated parameter time series and an object of the class "dyntermstrc_param" is returned. Using the summary() method leads to an overview about the correlation and possible unit roots of the estimated parameter time series for the levels and the first differences. The unit root tests are performed with functions form the urca package of Pfaff (2008).

Visualization of the results
The package contains several S3 plot() methods which allow us to visualize different interest rate curves, estimation errors, start parameter grid search results and the parameter estimates. The plot() methods for "termstrc_cs" and "termstrc_nss" objects offer the following possibilities: plot zero-coupon, forward, discount or spread curves, return single/multiple plots for the estimated group of bonds, show error plots for pricing/yield errors to identify outliers.
Plots of the zero-coupon yield curve for a single country include additional information. In detail, the yield to maturities are also plotted and for objects of the class "termstrc_cs", the knot points used for the estimation and the 95% confidence interval, based on robust standard errors, of the zero-coupon yield curve are added to the figure. The robust standard errors are calculated by applying functions of the package sandwich of Zeileis (2004). Both plot() methods depend on methods of the classes "spot_curves", "fwr_curves" and "df_curves" which are itself based on the plot() method of the class "ir_curve", which is the most granular class. An object of the class "spot_curves", "fwr_curves" or "df_curves" can therefore consist of several objects of the class "ir_curve". The advantage for the user is that when exploring a "termstrc_nss" or "termstrc_cs" object, she can plot the different curves at every hierarchical level of the object. For example, the command plot(<object>$spot$<COUNTRY>) creates a figure of the spot curve of <COUNTRY>, while plot(<object>$spot) plots the spot curves of all available countries. The plot() method for dynamic estimation results objects of the class "dyntermstrc_nss" and "dyntermstrc_yields" allows us to illustrate the evolution of the spot curves over time. The three-dimensional plots depend on the rgl package of Adler and Murdoch (2010). Objects of the class "dyntermstrc_param" contain the extracted time series of the estimated parameters. By applying the associated plot() method the time series of the levels and first differences of the parameters or the empirical autocorrelation function can be plotted. An inspection of the factor contribution plot enables the user to identify possible identification problems. The method fcontrib() is available for "dyntermstrc_param" objects. Moreover, the package offers a plot() method for "spsearch" objects, which include results on the start parameter grid search. Every static and dynamic results object except "termstrc_cs" includes a "spsearch" object. The command plot(<object>$spsearch$<COUNTRY>) plots the objective function of the start parameter grid search and a contour plot. The inspection of both helps us to understand the optimization problem in terms of multiple local minima.

Examples
In the following section we will estimate zero-coupon yield curves in R from various data sets with the afore mentioned methods. The examples are available as demos in the package.

Nelson/Siegel and Svensson method
After loading the package termstrc, we explore a static data set consisting of price information on coupon bonds from several countries.

R> class(govbonds)
[1] "couponbonds" It includes price data for government bonds of three European countries. The bonds are classified by their International Securities Identifying Number (ISIN), and all the necessary information on the future cash flows is provided. Next, we see the structure described in Table 2  Lets suppose we want to estimate the zero-coupon yield curve for the three included countries with the Nelson and Siegel (1987) method minimizing the duration weighted pricing errors. The sample of bonds is restricted to a maximum maturity of 30 years. We perform the estimation with the estim_nss() method available for the "couponbonds" class. The optional argument tauconstr defines the box constraint 0.2 < τ 1 ≤ 5 and a grid step size of 0.1 years for the start parameter search procedure. If not supplied by the user, economically meaningful constraints based on the available maturities in the data set are chosen.

R> class(ns_res$spsearch$GERMANY)
[1] "spsearch" R> par(mfrow = c(1, 3)) R> plot(ns_res$spsearch$GERMANY, main = "GERMANY") R> plot(ns_res$spsearch$AUSTRIA, main = "AUSTRIA") R> plot(ns_res$spsearch$FRANCE, main = "FRANCE") Next, we see the globally optimal estimated parameters for each country. They have the same order of magnitude, which implies similar shapes for the zero-coupon yield curves. The summary() method gives goodness of fit measures for the pricing and the yield errors, i.e., the root mean squared error (RMSE) and the average absolute error (AABSE). Moreover, it shows the convergence information from the solver to check whether a solution to the nonconvex optimization problem has been found.

R> summary(ns_res)
- Our package offers several options to plot the results, e.g., spot rate, forward rate, discount and spread curves. Figure 2 shows the estimated zero-coupon yield curves. The dashed lines indicate that the curve was extrapolated to match the longest available maturity.

Cubic splines
In this section, we demonstrate how to estimate the term structure of interest rates with the McCulloch (1975) cubic splines approach applied to French governement bonds.
R> cs_res <-estim_cs(govbonds, "FRANCE", matrange = c(0, 30)) R> cs_res FRANCE q q q q q q q q q q qq qq qq qq qq qq q q q q q q q q q q q q q q q q q q q q q q Zero−coupon yield curve 95 % Confidence interval (robust s.e.) Yield−to−maturity Knot points Compared to goodness of fit obtained in Section 4.1 for the French government bonds the cubic splines estimation leads to smaller errors. However, the number of estimated parameters is nearly twice as high and moreover assigning a clear economic interpretation to the parameters appears complicated. Therefore, fitting the zero-coupon yield curve could be preferable to the discount curve. Figure 3 shows the yield to maturities and the estimated zero-coupon yield curve together with the automatically selected knot points.

R> plot(cs_res)
As we can see from the plotted pricing errors in Figure 4, there seems to be a misspricing of two bonds.
R> plot(cs_res,errors="price") They can be removed by applying the rm_bond() method and the estimation is redone. As expected, the goodness of fit improves. q q q q q q q q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q Error (prices)

Rolling estimation procedure
In the last section, we will perform a Nelson/Siegel and Svensson estimation on data sets that store price/yield information for several days, i.e., for objects of the "zeroyields" or the "dyncouponbonds" class.

Zero yield data
The termstrc package includes a data set of German zero yields in CSV (comma-separated values) format. We load it and construct a "zeroyields" object.
For this class an estim_nss() method is available. We estimate all different model parameterizations, and by default a global start parameter search is performed for the first observation. Later time stages use the optimal parameter vector of the previous estimation as a starting value. With this procedure, the estimation tends to stay in a global optimum over time. As mentioned in Section 2.1, choosing an appropriate upper bound for the τ 1 and τ 2 parameter and a minimum distance between them, avoids identification problems and leads to smooth parameter time series. The argument tauconstr has the following structure: (lower bound τ 1 , τ 2 ; upper bound τ 1 , τ 2 ; grid step size; minimum distance ∆τ between τ 1 and τ 2 ). Performing the estimation with a start parameter search in each time step by setting optimtype = "allglobal" leads to only marginally different parameter values (for a small number of observations) in our examples. Therefore, we set optimtype = "firstglobal" as default value, and recommend performing a start parameter search only if an abrupt jump in the parameter values occurs or if the solver encounters convergence problems. We observed that for looser constraints on the τ -parameters, the solution can alternate between two optima. However, the improvement in the fit is negligible, and the smooth time evolution of the parameters is lost. Therefore, we set the upper bound tight enough to avoid alternating optima but ensured that the optimal solution does not actually hit the bound. In particular for the (adjusted) Svensson model, limiting the grid to a sufficient small size is desirable for performance reasons.

R> plot(param(ns_res))
Directly plotting the results object (ns_res) returns a 3D plot of the zero-coupon yield curve evolution created with the rgl package. We choose this over the usual persp() command,

Zero-Coupon Yield Curve Estimation with the Package termstrc
because it is very convenient to freely zoom in and out and rotate the plot. Figure 6 shows the estimated three-dimensional zero-coupon yield curve.

R> plot(ns_res)
Finally, we compare the goodness of fit for the four estimated spot curve parameterizations.

Coupon bond data
The data set GermanBonds contains daily price observations for German coupon bonds for several months.

R> plot(param(sv_res))
Finally, we also estimate the Adjusted Svensson version for the dynamic coupon bonds data set.
R> asv_res <-estim_nss(datadyncouponbonds, "GERMANY", method = "asv", + tauconstr = list(c(0.2, 10, 0.2))) [1] "Searching startparameters for GERMANY" beta0 beta1 beta2 tau1 beta3 tau2 5.907375 -4.544134 -7.227484 0.600100 -4.029217 5.600100 An interesting method for detecting possible overparameterization is plotting the factor contributions. Figure 9 shows this for the first observation in the dynamic data set. The first panel shows the separation in level, slope and curvature of the Diebold/Li model. Comparing this to the Nelson/Siegel method in the second panel could let us conclude that including more parameters would just lead to an identification problem due to the small contribution of the second factor. However, we can see in the third and fourth panel that the six parameter

Nelson/Siegel
Time to maturity Factor contribution

Adj. Svensson
Time to maturity Factor contribution τ 2 )) Figure 9: Factor contribution of several models. models split the individual factor contributions in a different way. This is primarily caused by the added flexibility for a second hump in the curve.

Conclusion
In this paper we discussed the two most widely-used parametric term structure estimation methods. We covered their estimation based on a set of coupon bonds which resulted in a nonconvex optimization problem with linear constraints. Due to the presence of multiple local optima, a global start parameter search algorithm was developed. The presented R extension package termstrc provides functions for the estimation of the zero-coupon yield curve based on static and dynamic coupon bond and zero yield data sets. Detailed examples illustrated the usage and functionality of the package.

A. Analytical gradients of the objective functions
We consider the two different objective functions defined in (5) and report upon the associated elements of the gradient ∇F (β). In the following "·" denotes an element wise multiplication. Divisions and the function exp(·) are also applied to each element of a matrix or vector. If not explicitly stated the elements are valid for the DL, NS, SV and adjusted SV model.
In the case of a term structure estimation based on coupon bond prices the elements of the gradient ∇F (β), with F (β) = p − ι (C · D) 2 ω , are given by  For a term structure estimation based on zero yields the elements of the gradient ∇F (β), with F (β) = (y − s(m, β)) 2 ι, are given by ∂ F ∂β 1 = (−2(y − b)) ι,