Tsukuba Economics Working Papers
No. 2010-008

Term Structure Models Can Predict Interest Rate Volatility.
But How?
by
Hideyuki Takamizawa
November 2010
This version: April 2011

UNIVERSITY OF TSUKUBA
Department of Economics
1-1-1 Tennodai
Tsukuba, Ibaraki 305-8571
JAPANTerm Structure Models Can Predict Interest Rate Volatility.
But How?

Hideyuki Takamizawa
Graduate School of Humanities and Social Sciences,
University of Tsukuba
Tsukuba Ibaraki 305-8571, JAPANE-mail: takamiza[at]social.tsukuba.ac.jp

April 17, 2011

Abstract
This paper attempts to predict the volatility of interest rates through dynamic
term structure models. For this attempt, the models are improved, based on thethree-factor Gaussian model, to have level-dependent volatilities supported by data.
The empirical results show that the predictive power of the proposed models is higher
than that of the aﬃne models. Compared with time-series models, it is low for the
four-week forecasting horizon but can be comparable for middle to long term rates byextending the horizon up to 32 weeks. The combination of these two diﬀerent typesof forecasts can lead to higher predictive power.

Keywords: Term structure, Stochastic volatility, Realized volatility, Approximation of
conditional moments.

1

1

Introduction

It does not seem unreasonable to think that the current yield curve contains some information on the volatility of changes in interest rates. In making bond portfolios or managing
interest rate risks, investors will take account of conditional second moments of bond returns or yield changes. The resulting shape of the yield curve will then reflect investors’
views toward the volatility. This paper attempts to predict the volatility by utilizing suchinformation through dynamic term structure models. For this attempt, the models are
improved without sacrificing the goodness-of-fit to the yield curve or the predictive power
to the level of interest rates. This improvement will contribute to extending the versatilityof the term structure models.
The idea of relating interest rate volatility to the yield curve is not new. Brown and
Schaefer (1994), Christiansen and Lund (2005), Litterman, Scheinkman, and Weiss (1991),
and Phoa (1997) relate the volatility to the curvature, or convexity, of the yield curve.
Time-series studies using long historical data on U.S. interest rates find a relation between
the volatility and the level of a particular yield, especially the short-term rate, such thathigh volatility is accompanied by high level; see, e.g., Andersen and Lund (1997a), Balland Torous (1999), Chan, Karolyi, Longstaﬀ, and Sanders (1992), Durham (2003), andGallant and Tauchen (1998).
This simple level-volatility relation, however, no longer seems to be a decisive feature
for relatively recent data. Figure 1 shows the time series of interest rates and realized
volatilities (annualized standard deviations) constructed from U.S. dollar LIBOR and swap
rates over 1991–2009: the details of the construction of the realized volatility measure areprovided in Section 2. We notice that the sharp rise in the volatility, observed around
2001–03 and 2008–09, is actually accompanied by the fall in the level of interest rates.
It is, therefore, not surprising that more recent studies using these data are skeptical
about the possibility of extracting volatility information from the yield curve. Andersen
and Benzoni (AB) (2010) test aﬃne spanning conditions that the yield variance, bothex ante and ex post, can be expressed by a linear combination of yields if aﬃne term
structure models are true, and reject these conditions. Collin-Dufresne, Goldstein, andJones (CDGJ) (2009), and Jacobs and Karoui (JK) (2009) report that yield variances
extracted from the cross-section of yields through aﬃne term structure models do not
behave similarly to typical variance measures in time series.

2

Nevertheless, it may be too early to draw the conclusion. Some important issues areleft unaddressed. The first issue is the predictability of the volatility over long forecastinghorizons: AB (2010) considered up to one-month horizons and CDGJ (2009) considered theone-week horizon. Longer horizons are worth considering, however. In predicting the levelof yields such as in the context of the eﬃcient hypothesis of the term structure (EHT),
the horizons longer than a month are often considered; see, e.g., Bekaert and Hodrick(2001), Campbell and Shiller (1991), Duﬀee (2002), and Fama and Bliss (1987). Then,if we match the horizon of the volatility forecast with that of the level forecast typicallyconsidered, the predictability of the risk-return relation implied by the yield curve can be
explored. The second issue is the use of non-aﬃne term structure models for extracting thevolatility. The aﬃne models used by CDGJ (2009), and JK (2009) have restrictions thatare necessary for obtaining a closed-form expression of bond prices and yields. Removing
such restrictions may improve the model’s ability of capturing the volatility.
To address the first issue, we perform the volatility forecasting regression with particular attention to long horizons and out-of-sample performance. We use daily data toconstruct weekly data on a realized yield variance measure over 4-, 8-, 16-, and 32-weekhorizons. It is then regressed on a linear combination of the current yields. The possibility of a nonlinear relation between them is also considered in a relatively simple way. We
make nonlinear transformations of the dependent variable, specifically the square-root andlogarithmic transformations, while leaving the independent variables linear. The results
are that the in-sample predictive power of the yield curve is not as bad as it is reported: it
improves the longer the forecasting horizon. On the other hand, the results of the out-ofsample predictive power of the yield curve are mixed. We prepare two approaches for theout-of-sample analysis. In the first approach where the parameter values of the regression
models are fixed at the in-sample estimates throughout the out-of-sample period, the predictive power deteriorates. In the second approach where the parameters are estimatedevery time the prediction is made, the predictive power still remains high. Taken together,at least some portion of volatility information seems to be contained in the yield curve,
but a simple linear combination of yields fails to produce a long-run predictive relation.
This motivates us to use dynamic term structure models for predicting the volatility.
To address the second issue, we attempt to develop term structure models without
sacrificing the goodness-of-fit to the cross-section of yields or the predictive power to3the level of yields. Such models may naturally be outside a class of models having aclosed-form expression of bond prices and yields, which will be a major obstacle to thisresearch. We overcome this obstacle by relying on an analytical approximation proposedby Shoji (2002), and Takamizawa and Shoji (2009). This method approximates a vector
of conditional moments as the solution to an ordinary diﬀerential equation (ODE). Since
the price of a zero-coupon bond is derived as the conditional expectation of the stochastic
discount factor, this method can be directly applied. The accuracy of the approximation,
which is carefully investigated also in this paper, seems to be maintained when realistic
values of parameter and state vectors are provided.
Now, what specification is appropriate for the volatility prediction without sacrificing
the desirable properties of the existing models? Clues to the answer can be found in theexisting models as well as in the results of the volatility forecasting regression. It is wellknown that the three-factor Gaussian model can capture many important features of theactual data, such as the failure of the EHT and a hump shape of unconditional volatilities
of changes in logarithmic yields; see Dai and Singleton (2003). An obvious drawback ofthis model is constant volatility despite the fact that time-varying volatility is one of thedecisive features of the data. We then overcome this diﬃculty by making the volatility inthis model level-dependent, leading to a model in which all factors potentially contributeto time-varying volatility.
To introduce the level-dependent volatility most eﬀectively, we carefully choose a baseline specification of the Gaussian model, which has numerous specifications depending on
factor rotation. One candidate is a classical specification, which is composed of the instantaneous risk-free rate, a stochastic central tendency factor, and an additional factor;see, e.g., Andersen and Lund (1997b), Balduzzi, Das, and Foresi (1998), Balduzzi, Das,
Foresi, and Sundaram (1996), and Bikbov and Chernov (2011). The former two factors
are expected to be highly correlated with short-term and long-term yields, respectively,
as is indeed the case shown later. Therefore, we can take advantage of the results ofthe volatility forecasting regression that a combination of yields is more or less useful for
predicting the volatility.
Another care in the introduction of the level-dependent volatility is to keep the instantaneous covariance matrix of the factors positive definite. One simple way to meet thisrequirement is to model eigenvalues of the covariance matrix as functions taking positive4

values. We propose two models for the time-varying eigenvalues: one given by quadratic
functions and the other by exponential functions. Apart from keeping the volatility positive, the positive definiteness of the covariance matrix has a significant merit. Since byconstruction the inverse of the covariance matrix is well defied, the market prices of risks
can be modeled as flexibly as those for the original Gaussian model, which is the so-called
essentially aﬃne specification; see Duﬀee (2002). Therefore, many of the advantages of theGaussian model are expected to be inherited into the proposed models, which is indeed
the case as shown later.
In comparing the volatility forecasting performance of the proposed models, we consider both aﬃne term structure models and time-series models. The latter include the
GARCH(1,1) model, which is estimated from lower frequency data. In addition, a variantof the AR model employed by AB (2010) and the mixed data sampling (MIDAS) approachdeveloped by Ghysels, Santa-Clara, and Valkanov (2005, 2006) are also considered, bothof which are estimated from higher frequency data. The results of the performance comparison are summarized as follows. First, compared with the aﬃne models, the proposedmodels exhibit a superior performance in most of the cases. Second, compared with thetime-series models, the proposed models are outperformed for the four-week forecastinghorizon. By extending the horizon up to 32 weeks, while the high performance of thetime-series models continues at short maturities, the relative performance of the proposedmodels increases. In fact, there are cases at middle to long maturities in which the predictive power of the proposed models is comparable, or even superior, to that of thetime-series models in both the in-sample and out-of-sample periods. Third, by combining
the forecast of the proposed models with that of the high-frequency regression model, thepredictive power can further increase at middle to long maturities. We also investigate
whether the improvement of the volatility prediction is attributable to the imposition of
non-arbitrage.
In spite of the improvement, the proposed models have a number of drawbacks. Apartfrom the computational complexity due to the lack of closed-form expressions for bondand derivative prices, we particularly point out two drawbacks. First, since the proposedmodels are based on the Gaussian model, negative interest rates occur with positive probability. Second, the proposed models do not accommodate factors that aﬀect derivative
prices but not underlying bond prices, in spite of the earlier work pointing out the exis-5tence and significance of such unspanned factors; see, e.g., Collin-Dufresne and Goldstein(2002), Han (2007), Heidari and Wu (2003, 2009), Jarrow, Li, and Zhao (2007), and Li
and Zhao (2006). We also discuss how to overcome these diﬃculties.
Section 2 explains the data and how to construct a realized volatility measure. Section
3 performs the volatility forecasting regression. Section 4 introduces our models, andexplains pricing and estimation methods. Section 5 examines the volatility forecasting
performance of the proposed models. Section 6 provides concluding remarks. Appendicescollect supplementary analyses as well as technical arguments regarding the derivation andaccuracy of the approximation method.

2

Data and realized variance measure

2.1

DatasetWe use data on U.S. dollar LIBOR and swap rates, covering the period from January 4,
1991 to May 27, 2009. There are mainly two reasons for choosing this dataset. First,
we can make the purpose of predicting the volatility through term structure models morechallenging. As mentioned in Introduction and exhibited in Figure 1, parallel movementsin the level and volatility of a particular yield disappear in the recent data. Second, we canfocus on the volatility prediction without introducing an additional complexity of regime
switching. As documented by Dai, Singleton, and Yang (2007), this sample period can be
regarded as a single regime when viewed from the history of U.S. interest rates.
The LIBOR and swap rates are transformed to zero-coupon bond yields on a continuously compounded basis using a bootstrap method with linear interpolation applied to
discount functions. The maturities of the zero yields used for the analysis are 0.5, 1, 2, 3,
5, and 10 years.

2.2

Realized measure

A realized measure of the conditional variance of changes in yields is constructed from thedaily data. We generate weekly, Wednesday, realized variance, where it is assumed that a
week begins with Thursday and ends with Wednesday. We choose Wednesday because welater estimate term structure models using weekly, Wednesday, data. We have 960 weekly
observations of the realized variance in total, among which first 640 observations (two-6thirds of the sample), ending on April 9, 2003, are used for the in-sample estimation andthe remaining 320 observations are used for the out-of-sample prediction. This division
is aimed at incorporating information on the lowest range of interest rates into model
estimation as well as reserving suﬃcient out-of-sample observations. Other ratios, threeto-one and four-to-one, are also tried, and the diﬀerences in the results are not large
enough to change the conclusions of this paper.
A realized measure of the one-week ahead conditional variance of a τ -maturity yieldis computed as

mt+∆

RVt,t+∆,τ =

�
i=1

(yt+

∆
i,τ
mt+∆

− yt+

∆
(i−1),τ
mt+∆

)2 ,

(1)

where
• yt,τ : a zero-coupon bond yield at time t with τ years to maturity.
• ∆: a week interval set to 1/52.
• mt : the number of observations during a week ending at time t (usually mt = 5).
A realized measure of the h-week ahead conditional variance is computed as
RVt,t+h∆,τ =h�RVt,t+j∆,τ .

(2)

j=1

The annualized realized variance is obtained by dividing RVt,t+h∆,τ by h∆. We consider
h = {4, 8, 16, 32}. The forecasting horizons of h > 4 are longer than previously considered.
The realized measure based on the daily data is crude. Discretization error arisingfrom approximating the integral of the quadratic variation may be large. This crude
measure, however, is more or less unavoidable because of the limited availability of intradaily data on interest rates; see CDGJ (2009), and JK (2009), where daily data are also
used to construct realized measures. Another concern is longer forecasting horizons. In
equations (1) and (2), the variation in the unpredictable component alone is reflected: thecovariation between the predictable and unpredictable components and the variation in
the predictable component are omitted, which are of smaller order than the variation in
the unpredictable component; see, e.g., Andersen, Bollerslev, Diebold, and Labys (2003).
The omission can be justified for a short horizon, say day or week, but it may not for a longhorizon considered in this paper. In particular, if the mean reversion, which is capturedby the predictable component, is crucial, the realized measure based on equations (1) and7

(2) will overestimate the true conditional variance since the mean reversion reduces theyield variation over a finite horizon.
In Appendix A, we verify the impact of omitting the mean reversion on the realizedmeasure by applying to each yield the Cox, Ingersoll, and Ross (CIR) (1985) model, which
has a closed-form expression of the conditional variance for a finite horizon. Specifically,letting β be a parameter in the drift controlling the speed of mean reversion, we compare
the conditional variance with β set at a point estimate from the actual data to that with β
set at zero. The summary of the results is as follows. The longer the maturity, the speedof mean reversion tends to increase, so does the resulting overestimation of the realizedmeasure. However, since the estimate of β is not significantly diﬀerent from zero, as isoften the case for interest rate data, the actual impact of regarding β to be zero is not as
serious as it appears.

3

Volatility Forecasting Regression

The purpose of this section is to examine if the current yield curve contains useful information for predicting the volatility of changes in yields. We perform the volatility forecasting
regression with particular attention to the predictability for longer forecasting horizons
and out-of-sample.

3.1

Setup

We regress the h-week ahead realized measure computed over time t and t + h∆ on linear
combinations of yields observed at time t:
Y-PART:f�

RVt,t+h∆,τh∆�= ah,τ + b�h,τ Ytp + upt+h∆,τ ,

(3)

where
Ytp = (yt,0.5
and
Y-ALL:f�

RVt,t+h∆,τ
h∆

yt,2
�

yt,10 )� ,= ah,τ + b�h,τ Yta + uat+h∆,τ ,where
Yta = (yt,0.5

yt,1

yt,2

8

yt,3

yt,5

yt,10 )� .

(4)

By considering a nonlinear function, f , the possibility of a nonlinear relation between
the realized measure and current yields can be addressed in a simple way. We actually√consider f (x) = x and f (x) = ln x, as well as f (x) = x as in the earlier studies taking
account of an implication of the aﬃne models. To conserve space, we do not report allof the results with diﬀerent f . But not surprisingly, the econometric specification is more
adequate on the standard deviation basis or the log variance basis than it is on the variance√basis. Below, we limit our attention to the standard deviation basis, i.e., f (x) = x, whichallows for a more intuitive interpretation of the results.
For comparison, we also consider the following regression models based on the timeseries of the realized measure:
AR-RV:

HAR-RV:f�

RVt,t+h∆,τ
h∆

�

f

�

RVt,t+h∆,τh∆�= ah,τ + bh,τ f= ah,τ +��

RVt−h∆,t,τ
h∆

�

+ uar
t+h∆,τ ,(5)bh,τ,i f��+ uhart+h∆,τ (. 6)i={4,8,16,32}

RVt−i∆,t,τi∆

The HAR-RV model is employed by AB (2010) and shown to have a superior performance.
It is noted that the data of the first 32 weeks are not used for the regression of the HARRV model. To make the comparison equal, these data are not used for the regression of
the other models.

3.2

In-sample results

We estimate equations (3)–(6) using the in-sample data. The in-sample forecasting performance is evaluated based on the adjusted R-squared coeﬃcient, R̄2 . Table 1 presents theresults. 1 First, for the four-week forecasting horizon, h = 4, the R̄2 coeﬃcients are highat short maturities but decrease with maturity. The result that the volatility of long-termyields is more diﬃcult to explain than the volatility of short-term yields is consistent withthe previous results; see AB(2010), and CDGJ (2009). Also of note is a small diﬀerencein the forecasting performance between the yield-based and RV-based regression models,although the HAR-RV model exhibits a slightly superior performance at both short andlong maturities. Second, the increase in h has diﬀerent impacts on the R̄2 coeﬃcients
between the yield-based and RV-based models. It is generally favorable for the former
1

Extending the in sample to the full sample does not materially change the results on the relative

ranking of the models.

9

but not for the latter, which is most evident at long maturities. At short maturities, theforecasting performance remains high for both types of models.

3.3

Out-of-sample results

High performance in-sample does not necessarily lead to high performance out-of-sample.
We perform an out-of-sample analysis in two approaches. In the first approach, the parameter values of the regression models are fixed at the in-sample estimates throughout
the out-of-sample period. This is aimed at examining whether the regression models can
produce a long-run predictive relation. In the second approach, the parameters are estimated every time the prediction is made in a rolling window fashion, where the size ofthe data for the estimation is fixed at the same as the in-sample data. This is aimed at
examining whether the current yield curve contains useful information for the volatility
prediction. In both approaches, the out-of-sample forecasting performance is evaluatedbased on the root mean squared error (RMSE).
Table 2 presents the RMSEs expressed in basis points (bps, 1bp = 0.01%) in the
fixed-parameter approach. Generally, the RV-based regression models, especially HARRV, have smaller RMSEs than the yield-based regression models. This is most evident at
τ ≤ 2: the diﬀerences in the RMSE between the Y-ALL and HAR-RV models reach 12.4
bps for h = 32 at τ = 1. The performance gap tends to narrow the longer the maturity
and forecasting horizon. Nevertheless, the performance deterioration of the yield-basedmodels is obvious, taking much higher R̄2 coeﬃcients of the in-sample regression into
consideration.2 The results indicate overfitting of the yield-based models to the in-sample
data. Once the parameter values are fixed, the models cannot respond to changes in theforecasting environment.
Table 3 presents the RMSEs in the varying-parameter approach. For h = 4, the HARRV model again exhibits the best performance, followed by the Y-ALL model. The RMSEs
for the Y-PART model are the largest, however, the diﬀerences in the RMSE between the2

The high R̄2 coeﬃcients of the yield-based models for long horizons seem to be caused by thehigh persistence of both the regressor and regressand. Interest rates on the right-hand side are known
to be persistent. The realized volatility on the left-hand side becomes more persistent by extending
the horizon since more overlapping daily observations are used in successive weekly observations of therealized volatility. For example, for h = 32, RVt,t+32∆ and RVt+∆,t+33∆ have overlapping observations
over 31 weeks.

10

Y-PART and HAR-RV models do not exceed 5 bps. The extension of the forecastinghorizon is again favorable for the yield-based models. Although the high performanceof the HAR-RV model continues at τ ≤ 2, the performance of the yield-based models
becomes comparable or superior at τ = {3, 5, 10}.
Taken both the in-sample and out-of-sample results together, the current yield curve
seems to fit the volatility realized in the future. But a simple linear combination of yields
fails to produce a long-run predictive relation, which after all seems diﬃcult to obtain without information on dynamic properties of the data. The time-varying parameter approachpossibly allows a simple linear combination of yields to incorporate such information, however, it lacks consistency. This motivates us to use dynamic term structure models forpredicting the volatility.

4

Model

The existing term structure models are known to have diﬃculties in matching both timeseries and cross-sectional properties of the data. Therefore, the purpose of predicting
the volatility through term structure models, while maintaining the cross-sectional fit,
is originally diﬃcult to achieve, and some breakthrough is required. In this section, we
first review the earlier work attempting to capture the volatility in order to see why
this purpose is diﬃcult and how our solution is diﬀerent. We then propose appropriate
models. Since the models do not have a closed-form expression of bond prices, we rely on
an approximation method.

4.1

A review of the volatility prediction using term structure models

A standard practice to extract volatility factors from the yield curve is to equate a term
structure model with a particular set of observed yields and to solve this equation for
latent factors including the volatility factors. In this inversion method, aﬃne term structure models are frequently used since they provide a one-to-one correspondence between
observed and latent variables in closed-form. To be more specific on the aﬃne models, weuse a notation introduced by Dai and Singleton (2000): Am (n) denotes an n-factor aﬃneterm structure model in which m factors aﬀect the instantaneous covariance matrix of then factors. To model stochastic volatility, m = 1 seems to be preferred as long as n ≤ 4.

11

This is because the dynamics of the m factors are restricted to be in a non-negative region,
which, however, restricts the correlation structure of the factors and the market prices ofrisks at the same time. In short, the more m, the more restrictions are necessary, which
deteriorates the model performance.
Nevertheless, the combination of the inversion method and the A1 (n) model does notwork well. CDGJ (2009) and JK (2009) report that the volatility implied by the A1 (3)model does not behave consistently with typical volatility measures in time series. This
inconsistency indicates that it is diﬃcult for the volatility factor to explain both the crosssection of yield levels and the time-series of yield changes. It is an outcome of the volatility
factor playing as a cross-sectional factor, rather than a time-series factor. The increase inn or m does not seem to be a solution, as long as the inversion method is employed withthe aﬃne models.
A solution taken by the earlier work for reducing the tension on the volatility factor
and avoiding the resulting inconsistency is to abandon the inversion method. Bikbov andChernov (2011) use the extended Kalman filter to obtain the time series of the volatility
factor. Thompson (2008) uses an estimation method where the volatility factor is integrated out from the likelihood function. A more extreme solution is to develop a term
structure model where the volatility factors do not appear in the cross-section of bondprices (but do appear in the cross-section of derivative prices). Collin-Dufresne and Goldstein (2002) propose the so-called unspanned stochastic volatility (USV) model, which isa special case of the Am (n) model by placing certain parameter constraints.
Our solution is diﬀerent from the previous ones. We maintain the inversion method,which has advantages of facilitating the factor identification as well as maintaining thecross-sectional fit with a relatively small number of factors. Since the combination of theinversion method and the existing models does not work well, we naturally abandon theexisting models. Our view is that the tension on the volatility arises because some factors
are presumed to be volatility factors. If instead such prespecification is not made, thetension and the resulting inconsistency will be naturally avoided. Our solution is thus toabandon the prespecification of any factor as the volatility factor. As a result, all factors
are supposed to explain the cross-section of yields as the A0 (n) model. This does not seemto be a bad idea, given a number of advantages of the A0 (n) model reported by manystudies. To overcome a major drawback of the A0 (n) model that the volatility is constant,12

we make the volatility of each factor level-dependent.

4.2

A baseline A0 (3) model

We build our models based on the A0 (3) model, and hence begin with its specification.
Among numerous specifications of the A0 (3) model depending on factor rotations, we
adopt a classical one characterized by the mean-reverting behavior of the instantaneousrisk-free rate; see, e.g., Andersen and Lund (1997b), Balduzzi et al. (1996, 1998), andBikbov and Cherenov (2005). Let Xt = (rt θt �t )� be a three-dimensional state vector, andthe risk-neutral distribution of the instantaneous change in Xt is given by κ1 (θt − rt ) + �t dXt ∼ N  κ2 (θ̄ − θt ) dt ,−κ3 �t

where









Σ dt

 σ11 σ12 σ13 





Σ=
 σ12 σ22 σ23  .

σ13 σ23 σ33





 ,



(7)

(8)

In this specification, θt is presumed to be the central tendency, toward which the instantaneous risk-free rate rt reverts. �t is presumed to be a shock that makes the mean-reverting
behavior a bit irregular. In reality, however, we have to place the following constraint a
priori, as pointed out by CDGJ (2008) and Joslin (2006), for each process to be interpreted
as above: 0 < κ2 < κ1 < κ3 . That is, θt also mean-reverts to some constant θ̄, but its
speed is slower than that of rt reverting to θt . The speed of mean reversion of �t is assume
to be the fastest, aiming at making the shock temporary.
This classical specification is beneficial in introducing the level-dependent volatility
in the next step. We can expect rt and θt to be highly correlated with short-term andlong-term yields, respectively, as is indeed the case ex post. In other words, the dynamics
of both ends of the yield curve can be more naturally described. Because at least twoof the factors are directly related to some yields and because a combination of yields has
some information useful for the volatility prediction as seen in Section 3, the introductionof the level-dependent volatility becomes more eﬀective within this classical specification.
To obtain the physical distribution of the instantaneous change in Xt for the purposesof estimation and prediction, the market prices of risks, Λ(Xt ), need to be specified.

13

Following the essentially aﬃne specification proposed by Duﬀee (2002),
Λ(Xt ) = Σ−1 (Λ0 + Λ1 Xt ) .

(9)

Then,
 
 
 

0
1

κ
(θ
−
r
)
+
�
λ

1
t
t
t


  1   λ11

 
 
 

 +  λ0  +  λ 1
dXt ∼ N  κ2 (θ̄ − θt )
  2   21


 
 


0
1


−κ3 �t

4.3



λ122




















λ123   θt  dt , Σ dt .(10)λ31 λ132 λ133

λ3




λ112 λ113   rt Proposed models

�tUsing A0 (3) as a baseline model, we make the instantaneous covariance matrix of changes
in Xt time-varying. But, first of all, it must be positive definite. We wish to meet thisrequirement without heavy parameterization. One possible way, similar in spirit to Fan,
Gupta, and Ritchken (2003), Han (2007), Jarrow et al. (2007), and Longstaﬀ, Santa-Clara,
and Schwartz (2001), is to take the spectral decomposition of the covariance matrix, andthen to specify the eigenvalues as functions of the state vector while holding the corresponding eigenvectors fixed. Specifically, the risk-neutral distribution of the instantaneouschange in Xt is given by






 κ1 (θt − rt ) + �t  dt , Σt dt ,







dXt ∼ N  κ2 (θ̄ − θt )

−κ3 �twhere Σt is decomposed as

Σt = P L t P � .

(11)

(12)

Lt is a diagonal matrix consisting of the eigenvalues, which are functions of the statevector: L1 (Xt )
Lt = 



0

0

0

L2 (Xt )000

L3 (Xt )





 .



(13)

P is an orthogonal matrix having the corresponding eigenvectors of unit length in its
columns. By the conditions of the orthogonality and unit length, the number of free14

parameters in P is actually three. We then parameterize P using three rotation matrices:


 1


0

P
P =
 0 cos ϕ3





P
  cos ϕ2

0




− sin ϕP3 

0 sin ϕP3

cos ϕP3



0 − sin ϕP2   cos ϕP1

0

sin ϕP2

1

00

cos ϕP2  sin ϕP1− sin ϕP1
cos ϕP1

0

0



0 


0 
 . (14)1This rotation is called the yaw, pitch, and roll rotation. The parameters to be estimatedare sin ϕPi (i = 1, 2, 3). For identification, we restrict ϕPi ∈ [−π/2, π/2], so that cos ϕPi =�

1 − sin2 ϕPi .

For Σt to be positive definite, all eigenvalues must be positive, i.e., Li (Xt ) > 0 for any

Xt (i = 1, 2, 3). Taking this into consideration, we propose two models for the time-varying
eigenvalues.
The first model, abbreviated as SV-Q (Stochastic Volatility with Quadratic specification), specifies Li (Xt ) as
Li (Xt ) = ci + Xt� Γi Xt(i = 1, 2, 3) ,(15)

where Γi is either a positive definite matrix if ci ≥ 0 or a nonnegative definite matrix if

ci > 0. In the estimation, we impose the latter restriction on Γi and ci as this can leadto a more parsimonious specification. For example, Γi = 0 is possible as long as the data
support. Similar to Σt , the non-negative definite matrix Γi is parameterized based on the
spectral decomposition:
Γi = Qi M i Qi�

(i = 1, 2, 3) ,

(16)

where


i
 m1


M =
 0

i

0

and


 1


0

0 sin ϕQ
3

mi2
0

i

0 

0 
 ,

mi3

 

0

Q
Qi = 
 0 cos ϕ3



0



i

cos ϕQ
3

i

i



Q
  cos ϕ2

 

i
 
− sin ϕQ
3  

 

0 ≤ mi1 ≤ mi2 ≤ mi3 ,

with

i

0

sin ϕQ
2

i

 

i

(17)

i

Q
0 − sin ϕQ
2   cos ϕ1   sin ϕQi
1
0
 
1

cos ϕQ
2

0

i

i

 

0

− sin ϕQ
1
cos ϕQ
1
0

i

i



0 


0 
 ,(18)1

Qwith ϕQj ∈ [−π/2, π/2] (j = 1, 2, 3). It is noted that sin ϕj cannot be identified for somei

mij . For example, when mij = 0 for all j, sin ϕQ
j cannot be identified for all j. In suchicases, we set sin ϕQ
j = 0.

15

The second model, abbreviated as SV-E (Stochastic Volatility with Exponential specification), specifies Li (Xt ) as
�

Li (Xt ) = exp ci + γ i� Xt

�(i = 1, 2, 3) .

(19)

No parameter restriction is required for the SV-E model. This exponential specificationthat naturally avoids negative volatility is very popular in time series analysis; see, e,g.,Andersen and Lund (1997a, b), Ball and Torous (1999), and Gallant and Tauchen (1998).
Since Σt is designed to be positive definite, the inverse matrix, Σ−1t , is well defined.
This is extremely important for modeling the market prices of risks. That is, they can be
specified in the same way as those for the A0 (3) model, i.e., equation (9) with Σ replaced byΣt . Likewise, the physical distribution of the instantaneous change in Xt is also obtainedby substituting Σt for Σ in (10). Therefore, the entire drift parameters can virtually be
replaced with new ones in changing the risk-neutral measure to the physical measure asin the A0 (3) model. Consequently, we expect that the proposed models inherit a numberof desirable properties from the A0 (3) model, which is verified in the next section.
At the same time, there are also a number of drawbacks in the proposed models.
Most critically, the models do not accommodate such factors as aﬀecting derivative prices
but not underlying bond prices. We discuss this issue and present a possible remedy in
Section 5.1. Another drawback is that the models generate negative interest rates withpositive probability as the original A0 (3) model. This problem, however, can actually be
avoided. As originally proposed by Aı̈t-Sahalia (1996), rt , starting from r0 > 0, does notreach zero in finite time by adding a reciprocal term to the drift of rt , i.e., by considering
a−1 /rt + κ1 (θt − rt ) + �t with a−1 > 0. If a suﬃciently small value is given to a−1 , say 10−8 ,the impact of the reciprocal term is negligible in the observed range of the data. Therefore,
while the empirical results shown below are based on the models without the reciprocal
term, virtually identical results will be obtained using the models with the reciprocal term.

4.4

Bond price and its approximation

Let P (Xt , τ ) be the price at time t of a zero-coupon bond with τ years to maturity. Then,by the absence of arbitrage, it is given by
P (Xt , τ ) = EtQ

�

�

exp −

16

� t+τ
t

ru du

��

,

(20)

where EtQ [.] stands for the conditional expectation under the risk-neutral measure. The
yield to maturity of the τ -year zero-coupon bond is given by Y (Xt , τ ) = − τ1 ln P (Xt , τ ).
The proposed models do not have a closed-form expression of P (Xt , τ ) due to thevolatility specifications. To make the estimation of the models and the extraction of thelatent state variables feasible, we approximate P (Xt , τ ) by relying on a method proposedby Shoji (2002), and Takamizawa and Shoji (2009). This method approximates a vector of
conditional moments as a solution to the ODE. Since the zero-coupon bond price is given
as the conditional expectation, this method can be directly applied. We actually use thethird-order approximation. The outline and application of this method are presented inAppendix B, and the accuracy of the approximation is checked in Appendix C. In short,the method seems to work at least for pricing bonds with maturities of up to ten years
when realistic values of parameter and state vectors are provided.

4.5

Estimation method

We estimate the models with the quasi-maximum likelihood method, which is one of thestandard methods for estimating term structure models; see, e.g., CDGJ (2008), Chen and
Scott (1993), Duﬀee (2002), and Pearson and Sun (1994). The zero-yields with maturities
of 0.5, 1, 2, 3, 5, and 10 years are used. Among them, those with maturities of 0.5, 2,
and 10 years are assumed to be explained exactly by the models in order to back out thelatent state variables. Specifically, let Θ be a parameter vector to be estimated, and let
Υ(Xt ; Θ) be a three-dimensional vector consisting of the approximate yield functions:�Υ(Xt ; Θ) = Ỹ (Xt , 0.5; Θ) Ỹ (Xt , 2; Θ) Ỹ (Xt , 10; Θ)

��

,

(21)

where Ỹ is for clarifying the approximation. Then, Xt is obtained by solving Ytp =Υ(Xt ; Θ). This is done numerically, however, only a few iterations are suﬃcient if a goodinitial value of Xt is given: it is actually the value of Xt implied by the A0 (3) model.
The rest of the yields, denoted as Yte = (yt,1denoted as Ut = (ut,1

ut,3

yt,3

yt,5 )� , are measured with error,

ut,5 )� : it is assumed to be independent from Xs for any s andfollow
Ut ∼ i.i.d.N (0, ς 2 I) .

(22)

The reason for assuming such a simple distribution is to let the models explain various
features of the data as much as possible.

17

The joint density function at time t conditioned on time t − ∆ can be written anddeveloped aspe; Θ)f (Ytp , Yte |Yt−∆, Yt−∆==��� dΥ(Xt ; Θ) �−1��f (Xt , Ut |Xt−∆ ; Θ) �
�
dXt�
�
�
� dΥ(Xt ; Θ) �−1��fT (Xt |Xt−∆ ; Θ)fC (Ut ; Θ) �
.
�
dX �

(23)

t

The first equality is from the changes of variables from Ytp to Xt , by which the Jacobian
term appears, and from Yte to Ut . The second equality is from the independence between
Ut and Xs for any s. For both the SV-Q and SV-E models, the transition density, fT , has
no analytical expression for finite ∆. We approximate it by the multivariate normal density
function, which might be justified by a relatively short interval, ∆ = 1/52. The conditional
first and second moments to be substituted are computed with the same approximation
method. It is noted that for the SV-Q model, these moments can be computed exactly
since the drift terms are linear and instantaneous (co)variances are quadratic in the statevariables. The Jacobian term is already computed in backing out the latent state variables.
On the other hand, fC is the multivariate normal density function from (22). The objectivefunction for estimating Θ is then
��
t

5

�
��
� dΥ(Xt ; Θ) �
� .
�
dX �

ln fT (Xt |Xt−∆ ; Θ) + ln fC (Ut ; Θ) − ln ��Main Empirical Analysis

(24)

t

Parameter estimates and time series of implied state variables are first reported. Wethen examine whether the proposed models can reproduce some important properties of
the actual data: the failure of the EHT and a hump shape of unconditional volatilities
of changes in logarithmic yields. Furthermore, we examine the predictive power of the
proposed models to the level of yields. After confirming that the proposed models have atleast a comparable performance to the A0 (3) model in these dimensions, we finally proceed
to examining whether the predictive power to the volatility is improved.

5.1

Parameter estimates and implied state variables

We first report the parameter estimates for the baseline A0 (3) model in Table 4. It is notedthat when the full parameters are estimated, some parameters, particularly those of themarket prices of risks, are insignificant. We set such parameters to zero, and re-estimate

18

the model until all t-values exceed 1.5 in absolute value. As a result, the instantaneous
covariance between θt and �t is zero: σ23 = 0. Also, θt alone contributes to the variationin the market prices of risks: the first and third columns of Λ1 are all zeros. This result
is consistent with Joslin, Singleton, and Zhu (2011, Table 5), where the constraint ofrank(Λ1 ) = 1 is not rejected. The speed of mean reversion, captured by the estimate of
κi , is by far the lowest for θt . It is similar between rt and �t , suggesting that the threefactors are suﬃcient for describing the cross-section of yields. This is also suggested by arelatively small standard deviation of the measurement errors, 6.1 bps.
Figure 2 exhibits the time series of the state variables implied by the A0 (3) model.
As is indeed anticipated, the dynamics of rt and θt are very closely tracked by those of
the six-month and ten-year yields, respectively. It is so even in the out-of-sample period,
implying that there is no structural break between the samples. �t ranges between −0.07
and 0.04, and looks more volatile. This is confirmed by the estimate of σ33 shown in Table
4, which is much larger than the estimates of σ11 and σ22 . In sum, the A0 (3) model isshown to be appropriately specified in that the ex post dynamics of the state variables are
consistent with the ex ante roles.
Tables 5 and 6 present the parameter estimates for the SV-Q and SV-E models, respectively. Similar to the estimation of the A0 (3) model, insignificant parameters are set
to zero in the final estimation.3 It is also noted that only the same parameters of Λ(Xt )as those for the A0 (3) model are estimated, which is aimed at distinguishing the models
solely by the specification of the covariance matrix. We notice that the estimates of therisk-neutral drift are very similar across the three models. For both the SV-Q and SV-E
models, the estimate of sin ϕP2 alone is significant in the eigenvector matrix P . As a result,the instantaneous covariance matrix for these models isΣt = 

L1 (Xt ) cos2 ϕP2 + L3 (Xt ) sin2 ϕP20{L1 (Xt ) − L3 (Xt )} sin ϕP2 cos ϕP20

L2 (Xt )0{L1 (Xt ) − L3 (Xt )} sin ϕP2 cos ϕP2

0

L1 (Xt ) sin2 ϕP2 + L3 (Xt ) cos2 ϕP2 .(25)

It turns out that L2 (Xt ) is the instantaneous variance of θt . In addition, because sin ϕP2 <cos ϕP2 , L1 (Xt ) is more closely related to the instantaneous variance of rt , whereas L3 (Xt )is to the instantaneous variance of �t . Among the instantaneous covariances, only that3

As exceptions, for the SV-Q model, c1 is set at a very small number, 10−8 , due to the constraint of
2

2

c1 > 0, and sin ϕQ= − sin ϕQ
3
2 .

19

between rt and �t is significant.
The estimates of ς, the standard deviation of the measurement errors, for the SV-Q
and SV-E models are 6.1 bps, which are little changed from that for the A0 (3) model.
The result reconfirms that as far as the description of the cross-section of yields, theA0 (3) model has already done a suﬃcient job. Time-varying volatility does not seem tobe eﬀective in this dimension. It is, however, in the time-series dimension of the data.
Indeed, compared with the A0 (3) model, the log-likelihood value for the SV-Q model isincreased by 152 with additional 10 parameters and that for the SV-E model is increasedby 95 with additional 7 parameters.
Figure 3 exhibits the time series of the state variables implied by the SV-Q (leftpanels) and SV-E (right panels) models. For each state variable, two implied series arevisually indistinguishable. Furthermore, they are almost identical to that implied by theA0 (3) model shown in Figure 2. The result also indicates that the specification of theinstantaneous covariance matrix little aﬀects the behavior of the state variables extracted
from the cross-section of yields, given the identical specification of the risk-neutral drift.
In other words, it is the risk-neutral drift that matters with the identification of the statevariables.
But if this is indeed the case, the following question naturally arises: can a statevariable be identified from the cross-section of yields that appears in the instantaneous
covariance matrix but not in the risk-neutral drift? If not, this state variable would
virtually act as an unspanned factor. Indeed, a similar argument is made by Joslin (2010)
using aﬃne term structure models. 4 In addition, in CDGJ (2009, p.51), the suﬃcientconditions to obtain the A1 (4) USV model from the A1 (4) model are presented, one ofwhich indicates that the coeﬃcient of a volatility factor in the risk-neutral drift term of aconditionally Gaussian factor should be very small. In fact, their Table 2A shows that theestimate of this coeﬃcient for the A1 (4) USV model is −0.260, which is in sharp contrast
4

Specifically, the convexity arising from the Jensen’s inequality can be defined as
1
CV (Xt , τ ) ≡
τ

� t+τ
t

EtQ [r(Xu )]du − Y (Xt , τ ) .

But Joslin (2010) showed that both the level and variation of CV (Xt , τ ) are small even at τ = 10. This
indicates that Y (Xt , τ ) is mostly determined by

� t+τ
t

EtQ [r(Xu )]du/τ and that this relation is stable overtime. Then, if the volatility factors are not included in EtQ [r(Xu )], they cannot have an impact on Y (Xt , τ )large enough to be backed out from Y (Xt , τ ) with some precision.

20

with the corresponding estimate of −809.5 for the A1 (4) model obtained by the inversion
method. Other parameters in the risk-neutral drift terms do not exhibit such a dramaticchange. The result indicates that by this restriction alone the volatility factor becomes
nearly, if not completely, unspanned. A similar result is reported by Thompson (2008,
Tables 3 and 4). Taken together, the introduction of a nearly unspanned factor does notseem as diﬃcult as previously thought: simply exclude a factor from the risk-neutral drift
terms.

5.2

Ability of the models to replicate real properties of the data

We have seen many similarities between the A0 (3) and proposed models. We then expectthat the advantages of the A0 (3) model are inherited into the SV-Q and SV-E models.
We verify this expectation. Specifically, we examine whether these models can replicate
the failure of the EHT and the hump shape of unconditional volatilities of changes inlogarithmic yields.
5.2.1

Failure of the EHT

Following Campbell and Shiller (1991) but using a slightly diﬀerent notation, we consider
the following regression model:
yt+τ1 ,τ2 −τ1 − yt,τ2 = φ0τ2 + φ1τ2

τ1(yt,τ2 − yt,τ1 ) + ut+τ1 ,τ2 − τ 1(26)

where τ1 and τ2 are measured in years and τ2 /τ1 is integer. Since we use weekly data,
we set τ1 = ∆ (= 1/52). But, this requires the one-week yield. Instead of generating theone-week yield from the original LIBOR and swap curves by extrapolation, we use as aproxy the six-month yield. This simplification might be justified because our purpose is
to examine whether the proposed models can replicate an empirical pattern rather thanto test the EHT rigorously. Also, we replace yt+∆,τ2 −∆ on the left-hand side of equation(26) with yt+∆,τ2 for convenience. We consider τ2 = {1, 2, 3, 5, 10}.

It is know that the actual estimates of φ1τ2 tend to be negative, and more so withincreasing τ2 , contrary to the EHT that φ1τ2 = 1 holds for all τ2 . We estimate equation(26) using both actual and simulated data. For each model, we simulate 1,000 sets of datawith the same sample size as the actual in-sample data. Precisely, the initial realization
is drawn randomly from the actual data. The subsequent realizations are generated from21the physical distribution of changes in Xt given by (10): for the SV-Q and SV-E models,Σ is replaced by Σt . For the SV-Q and SV-E models, since analytical expressions of
the distribution over a finite interval are unavailable, we use a discretized version of the
distribution, where dt is replaced by ∆t = 1/1, 000. We then obtain 1,000 estimates of
φ1τ2 , and record the average, and the 5th and 95th percentiles of these estimates.
The left panels of Figure 4 exhibit the actual (the solid line) and model-implied (themarked line) estimates of φ1τ2 together with the 90% confidence band (the dotted line).
While the extent of the failure of the EHT is weaker for the simulated data than for theactual data, all three models have a similar pattern, indicating that the proposed models
can regenerate the failure at least as well as the A0 (3) model.
5.2.2

Hump shape of unconditional volatilities

The right panels of Figure 4 exhibit the actual (the solid line) and model-implied (themarked line) estimates of unconditional volatilities (annualized standard deviations) ofchanges in logarithmic yields with the 90% confidence band (the dotted line). Again, allmodels can regenerate a hump shape that is close to the actual one.

5.3

Predictive power to the level of yields

Duﬀee (2002) shows that the essentially A0 (3) model having the market prices of risks
given by equation (9) has a superior predictive power to the level of yields than the otheressentially Am (3) models with m > 0. We examine whether the proposed models haveat least a comparable ability in this dimension. To compute the conditional expectationof the h-week ahead τ -maturity yield for the proposed models, Et [Ỹ (Xt+h∆ , τ )], we also
employ the approximation method of conditional moments: the details of the application
are provided in Appendix B and the accuracy analysis is in Appendix C. To match theforecasting horizons with those for the volatility, we also consider h = {4, 8, 16, 32}.
Table 7 presents the in-sample RMSEs in basis points, where the prediction error
is defined as the diﬀerence between predicted and realized yields. As a benchmark, we
provide the results for random walk (RW), whose forecast is the current yield for any h.
The forecasts of the three term structure models are generally better than those of RW,
except for h = {8, 16} at τ = 1. Among the term structure models, the performance issimilar, though the RMSEs for the SV-Q model tend to be slightly larger.

22

Table 8 presents the out-of-sample RMSEs in basis points. Generally, the three term
structure models have smaller RMSEs than RW. At τ = 10, however, RW is comparable
for h ≤ 16 and superior for h = 32 to the term structure models. The result indicatesthe diﬃculty of predicting the ten-year yield in terms of both the level and volatility.
Among the term structure models, while the proposed models are slightly better at shortmaturities, the performance is generally similar, confirming that the ability of the proposedmodels to predict the level of yields is maintained.

5.4

Predictive power to the volatility of yields

We consider several competing models. First, we select the A0 (3) and A1 (3) models.
Since, given h and τ , the A0 (3) model produces a constant forecast independently fromthe current state, it serves as a benchmark, as does RW for the level prediction. As forthe specification of the A1 (3) model, we select one proposed by CDGJ (2008, p.762). The
A1 (3) model is estimated in the same way as the other term structure models.
Second, we consider a time-series model applied to the same weekly data. We selectthe GARCH(1,1) model. Although there are a number of variants of the GARCH model,
it is selected because more complicated models do not necessary beat the simplest one in
the out-of-sample period; see Hansen and Lunde (2005). The model is fitted to each yield,
yt,τ :
�yt+∆,τ

= α + βyt,τ +Vt+∆,τ zt+∆,τ ,Vt+∆,τ

2= ω0 + ω1 Vt,τ zt,τ+ ω2 Vt,τ ,

(27)
(28)

where it is assumed for simplicity that zt,τ ∼ i.i.d.N (0, 1). The h-week ahead conditionalvariance, vart [yt+h∆,τ ], is obtained by iteration, an explanation of which is provided inAppendix D.
Third, we consider time-series models applied directly to the daily data. We select the
HAR-RV model presented in Section 3 and the mixed data sampling (MIDAS) approachdeveloped by Ghysels et al. (2005, 2006). The MIDAS regression in this paper is asfollows:
MIDAS:

�

RVt,t+h∆,τ= ah,τ
h∆
+bh,τ

��

nt
2
i=1 (yt−L∆(i−1)/nt ,τ − yt−L∆i/nt ,τ ) wh,τ (i)(L∆/nt )23+ um
t+h∆,τ , (29)

where following Ghysels et al. (2005, equation (3)), the weighting function, wh,τ (i), is
given byexp{wh,τ,1 i + wh,τ,2 i2 }wh,τ (i) = �nt,2
k=1 exp{wh,τ,1 k + wh,τ,2 k }(30)

and nt stands for the number of daily observations over the past L weeks: usually nt =L × 5. We set L = 32 to make the amount of information equal between HAR-RV and
MIDAS. The diﬀerence between the two is in the weighting structure to past observations.
The predictive power to the volatility is evaluated by the RMSE criterion using boththe in-sample and out-of-sample data. The RMSE for the term structure and GARCH
models is computed based on the following equations:��

RVt,t+h∆,τh∆= f orecastt,h,τ + u1t+h∆,τ ,(31)

RVt,t+h∆,τh∆= ah,τ + bh,τ f orecastt,h,τ + u2t+h∆,τ ,(32)

where f orecastt,h,τ is a model forecast at time t for the h-week ahead volatility of theτ -maturity yield. Some explanation is needed as to how to compute f orecastt,h,τ forthe proposed models. First, while we assumed that the yields not used for backing out
the factors are measured with error in estimating the models, we omit the variance of
this error in computing the conditional variance of these yields, as it is a small constant.
That is, we assume here that yt,τ = Ỹ (Xt , τ ) for all τ . Second, since f orecastt,h,τ =�vart [Ỹ (Xt , τ )] has no closed-form expression, we employ the approximation method of

conditional moments again; see Appendices B and C for the application and the accuracy,
respectively.
For the HAR-RV and MIDAS regressions, the RMSE is computed from equations (6)and (29), respectively.
As in Section 3, the out-of-sample predictive power is evaluated in both the fixed- andvarying-parameter approaches. The out-of-sample RMSE in the fixed-parameter approachis computed with the parameter values held fixed at the in-sample estimates throughout
the out-of-sample period. The out-of-sample RMSE in the varying-parameter approachis computed as follows. First, for the term structure and GARCH models, the model
parameters are fixed at the in-sample estimates. The parameters that vary every time theprediction is made are (ah,τ , bh,τ ) in equation (32). Second, for the MIDAS regression,(wh,τ,1 , wh,τ,2 ) in equation (30) are fixed at the in-sample estimates, while (ah,τ , bh,τ ) in24

equation (29) are re-estimated. Finally, for the HAR-RV regression (6), all parameters are
re-estimated as before.
It is noted that the above comparison scheme is challenging to the term structure
models. There is an information gap between the term structure and GARCH models.
While each yield is used for estimating the GARCH model, only three yields, the sixmonth, two-year, and ten-year yields, are used for estimating time-series properties of theterm structure models. The information gap widens further against the HAR-RV and
MIDAS regressions, which are constructed directly from the realized volatility.
5.4.1

In-sample results

To obtain an intuition about the model performance, we first look at Figures 5–8, where
the time series of the four-week ahead forecast of the volatility (annualized standard deviation) of the six-month, two-year, and ten-year yields are displayed together with thecorresponding realized measure. In Figure 5, the GARCH forecast tracks the realized measure well with the correlation coeﬃcients being around 0.6, though it tends to overestimate
the level of the volatility of the six-month yield and underestimate the variation of thevolatility of the ten-year yield. In Figure 6, the A1 (3) forecast fluctuates little. Besides,
the correlations between the predicted and realized series are low: for the ten-year yield,
it is −0.19. The results are consistent with those reported by CDGJ (2008, 2009), and JK(2009). In Figure 7, the SV-Q forecast is apparently better than the A1 (3) forecast withthe correlation coeﬃcients being around 0.5, yet it does not truck the realized measure
as closely as the GARCH forecast. More specifically, the forecast is occasionally wrong
especially during 2002–04, however, it can remarkably capture the sharp rise and fall in therealized measure observed during 2008–09 even though this period is out-of-sample. At
τ = 10, while the correlation remains reasonably high, the variability of the model forecast
is much smaller than that of the realized measure. Finally, in Figure 8, the SV-E forecast
is also improved from the A1 (3) forecast, but it does not track the realized measure as
closely as the GARCH forecast, either. Compared with the SV-Q forecast, the correlation
at τ = 0.5 is slightly higher but those at τ = {2, 10} are lower. Still, the SV-E forecast at
τ = 10 appears to track the trend of the realized measure better than the SV-Q forecast.
More generally, the SV-E forecast appears to be more persistent. It produces neither large
deviations nor intensive variations. In fact, the serial correlations for the SV-E forecast

25

are all around 0.99, which are larger than those for the SV-Q forecast, around 0.95.
Having these figures in mind, we compare the in-sample model performance. Table 9
presents the RMSEs in basis points. For the term structure and GARCH models, the namealone indicates the results based on equation (31), whereas the name “+ Reg” indicatesthe results based on equation (32). Note that there is no A0 (3) + Reg, as f orecastt,h,τ
is constant. In reporting the results, we focus on the comparison first without then withthe forecasting regression (32).
For the four-week forecasting horizon, h = 4, the proposed models predict the volatility
better than the aﬃne models, except at τ = 10, where the benchmark A0 (3) model
has a smaller RMSE. Compared with the GARCH model, the RMSEs for the proposedmodels are larger, except for the SV-E model at τ = 0.5. The proposed models arealso outperformed by the HAR-RV and MIDAS regressions especially at short maturitiesτ ≤ 1. However, the performance gap narrows with increasing τ . Even with the forecasting
regression (32), the RMSEs for the proposed models are larger than those for the HAR-RV
and MIDAS regressions. But the diﬀerences in the RMSE are reduced to around 1 bp
at middle to long maturities. The performance of the GARCH + Reg model is almost
indistinguishable from that of the HAR-RV and MIDAS regressions.
By extending the forecasting horizon, the relative performance of the proposed models
improves. While at τ ≤ 1, the HAR-RV and MIDAS regressions continuously exhibit
a high performance, the performance of the proposed models becomes comparable forh = {16, 32} at τ = {2, 3, 5}. With the forecasting regression, this holds at all τ . In fact,
the RMSEs for the SV-Q + Reg model are the smallest for h = {8, 16, 32} at τ = {2, 3},
and those for the SV-E + Reg model are the smallest for h = {16, 32} at τ = {5, 10}.
5.4.2

Out-of-sample results

We compare the out-of-sample model performance first in the fixed-parameter approach,where the parameter values are fixed at the in-sample estimates throughout the out-ofsample period. Table 10 presents the RMSEs in basis points. For h = 4, the SV-E model
outperforms the aﬃne models for all τ , whereas the SV-Q model fails to outperform theA1 (3) model at τ ≤ 1. However, the SV-Q model improves with increasing τ . In fact, it has
the smallest RMSE at τ = 10 among all models. Compared with the GARCH model, theproposed models have a comparable or superior performance at both ends of the maturity

26

spectrum but are worse at in-between maturities. They are outperformed by the HAR-RV
and MIDAS regressions overwhelmingly at short maturities but become comparable atτ = 10. With the forecasting regression, the performance of the proposed models becomes
better at τ ≤ 1 but worse at τ = {3, 5, 10} than it is without the forecasting regression.
Since the decrease in the RMSE at τ ≤ 1 is insuﬃcient, the proposed models with theforecasting regression are outperformed by the HAR-RV and MIDAS regressions at all τ .
By extending the horizon, while the HAR-RV and MIDAS regressions continuously
exhibit a superior performance at τ ≤ 1, the relative performance of the proposed models
improves at middle to long maturities, as is the case for the in-sample prediction. For
h = 8, the SV-Q model is comparable or superior to the HAR-RV and MIDAS regressionsat τ = {3, 5, 10}, although it is worse than the A1 (3) model at τ ≤ 1. The SV-E model
also outperform the high frequency regressions at τ = 10. This pattern is maintained,
but somewhat strengthened, for h = {16, 32}. As a result, the SV-Q model exhibits
a comparable performance in both sample periods to the time-series models for h ={8, 16, 32} at τ = {3, 5, 10}. The forecasting regression again increases the out-of-sample
RMSEs for the proposed models at middle to long maturities. But the degree of increase
is relatively minor. As a result, the SV-Q + Reg model is comparable in both sample
periods to the time-series models for h = {16, 32} at τ = {2, 3, 5, 10}. Likewise, the SV-E
+ Reg model is comparable for h = {16, 32} at τ = 10.
Next, we compare the out-of-sample model performance in the varying-parameter approach, where the regression parameters are re-estimated every time the prediction ismade. Table 11 presents the RMSEs in basis points. For h = 4, the SV-Q model outperforms the A1 (3) model at all τ , whereas the SV-E model fails to do so at τ ≤ 1. The
proposed models are outperformed by the time-series models, however. In contrast withthe previous performance criteria, the relative performance of the proposed models does
not much improve by extending the horizon. This is partly because the performance ofthe HAR-RV regression remains high due to the advantage of more varying parameters: it
has five varying parameters, not two like the other models. Still, there are cases in which
the diﬀerences in the RMSE between the SV-Q and HAR-RV models are within 2 bps: for
h = 8 at τ = {3, 10} and for h = {16, 32} at τ = 3. Compared with the GARCH model,
the proposed models are outperformed but catch up at τ = 10.

27

5.4.3

A possible economic reasoning for the results

Why is the predictive power of the time-series models especially high for short forecastinghorizons at short maturities and why does the predictive power of the proposed models
relatively improve for long forecasting horizons at middle to long maturities?
The short-term yields are closely related to policy interest rates controlled by thecentral bank. Usually, the policy rates do not change frequently in a short period oftime. Even if they do, they rather move in an anticipated direction. Therefore, theshort-term yields can be more easily predicted, in terms of both the level and volatility,
using information on the past state of the economy. The time-series models are good at
exploiting such information, and thus exhibit a superior forecasting performance for shorthorizons at short maturities.
On the other hand, the role of the expectation about the future state of the economy
increases not only in making longer-term forecasts but also in pricing longer-term bonds.
By the specifications of the covariance matrix, Σt , the volatility forecast of the proposedmodels depends on the current yield curve, which by nature contains forward-looking
information. Thus, an adequately specified term structure model has potential to have arelatively high predictive power for long forecasting horizons at middle to long maturities.
These results suggest that a combination of the two diﬀerent types of forecasts, one
from the time-series models and the other from the term structure models, may further
improve the predictive power. This possibility is explored in the next subsection. Also ofinterest is whether the non-arbitrage constraint, given specifically by equation (20), hassome contribution to the predictive power, beyond the specifications of Σt . This questionis addressed in Section 5.6.

5.5

Combined forecasts

We make a combined forecast based on the HAR-RV model, which has a simple structureand exhibits a superior performance especially at short maturities:��
RVt,t+h∆,τ= ah,τ +bh,τ,ih∆
i={4,8,16,32}

�

RVt−i∆,t,τ
+ ch,τ f orecastt,h,τ + uct+h∆,τ , (33)i∆where f orecastt,h,τ is taken from the term structure and GARCH models. Then, we
examine the sign, magnitude, and significance of the estimate of ch,τ using the in-sample

28

data. Furthermore, we examine the extent to which the combined forecast reduces theout-of-sample RMSE relative to the HAR-RV forecast alone.
Table 12 presents the estimates of ch,τ in equation (33) together with the standard
errors in parenthesis computed using the method of Newey and West (1987) with 32 lags.
First, for the four-week forecasting horizon, h = 4, the slope estimates for the SV-Q forecast lie in a relatively narrow range from 0.36 to 0.47, which are all significantly diﬀerent
from zero but not from 0.5. Those for the SV-E forecast are also all significant with therange of estimates somewhat wider. The GARCH forecast has significant estimates except
at τ = {0.5, 2}. The result indicates that incremental information of the GARCH forecast
is relatively small. This is not surprising given the similarity between the HAR-RV and
GARCH models in the way of using time-series information to make forecasts. The slope
estimates for the A1 (3) forecast are significant only at τ ≤ 1. Though insignificant, theestimates at τ = {5, 10} are negative. A similar pattern is observed for longer horizons,
with a notable exception that the estimates for the proposed models at τ = 10 tend to beless significant.
Table 13 presents the % changes in the out-of-sample RMSE in the fixed-parameter
approach. Of particular attention is the success of the HAR-RV + SV-Q forecast. This
combined forecast reduces the out-of-sample RMSEs at middle to long maturities for all h,confirming an incremental value of the SV-Q forecast. Adding the GARCH forecast also
reduces the RMSEs at middle to long maturities, however, the magnitude of reductionis less than that achieved by the SV-Q forecast. On the other hand, adding the SV-E
forecast worsens the out-of-sample performance except at τ = 10. Adding the A1 (3)forecast generally deteriorates the performance with the increase in the RMSE reaching
30% for h = 32 at τ = 1.
Table 14 presents the % changes in the out-of-sample RMSE in the varying-parameter
approach. The overall picture agrees with that in Table 13. Again, the performance ofthe HAR-RV + SV-Q forecast is remarkable. In all cases except for h = 32 at τ = 0.5, it
reduces the RMSEs by up to 16%. Adding the SV-E forecast, however, leads to the worst
out-of-sample performance in this criterion: the RMSEs increase in all cases.
Taken together, the combined forecast based on the HAR-RV model works especiallyat middle to long maturities, as expected. But the choice of the other model matters. In
our sample, the SV-Q model is an appropriate one. The GARCH model also contributes

29

to the reduction in the RMSE, but not as much as the SV-Q model.

5.6

Is the non-arbitrage relevant with the volatility prediction?

Recently, Joslin et al. (2011) have proposed a canonical representation of the Gaussian
term structure model in which linear combinations of observed yields serve as underlying
factors. Since there is no need to back out the factors from the cross-section of yields, thephysical drift terms of the factors, which involve the prediction of the level of factors, arefree of any cross-sectional constraint including the non-arbitrage. In predicting the level ofthe other yields not used for making factors, some cross-sectional constraint relating thefactors to these yields is required. But Duﬀee (2011) shows that a simple linear regression
is suﬃcient for constructing this relation, so that the more complicated non-arbitrage
condition is unnecessary for the level prediction of both yields used and not used formaking factors.
The above studies have shed new light on the relation between the role of non-arbitrage
and the level prediction of yields, using the Gaussian term structure model. However, lessclear is whether this relation also holds for the prediction of the volatility that changes
stochastically. Since the covariance matrix, Σt , is invariant to measure changes, it isrelated to both time-series and cross-sectional properties of the yield curve once the nonarbitrage is imposed, though the relation to the latter may be weak as argued in Section
5.1. Therefore, the estimate of Σt with the non-arbitrage might diﬀer from that without the
non-arbitrage obtained using time-series properties of the data alone, which would produce
the diﬀerence in the volatility prediction between with and without the non-arbitrage.5
To examine the impact of non-arbitrage on the volatility prediction, we first clarify therole of non-arbitrage in this particular context:(i) description of the cross-sectional relation between yields,(ii) estimation of Σt ,(iii) identification of factors.
According to these roles, we control the degree of impact in the following three stages.
In the first stage, we take (ii) and (iii) as given, and remove (i). Specifically, we use the5This is the case for options data on stocks and currencies; see, e.g., Bakshi, Cao, and Chen (1996),

and Bates (1996, 2000).

30

time-series of the state vector, Xt , extracted through the SV-Q and SV-E models withthe parameter values fixed at the estimates presented in Tables 5 and 6, respectively. The
cross-sectional relation is described by a linear regression as�
yt,τ = δ0,τ + δ1,τXt + ucross
.
t,τ

(34)

� ), are estimated by OLS using the in-sample data. These estiThe coeﬃcients, (δ0,τ δ1,τmates are fixed throughout the out-of-sample period to maintain the consistency with theprevious analysis in Section 5.4. But it is noted that even though they are changed every
time the prediction is made in the out-of-sample analysis, the results are not changed
much, indicating that the cross-sectional relation is stable over the whole sample period.
In computing the h-week ahead conditional variance, we omit the variance of ucross
as in
t,τ
Section 5.4:�vart [yt+h∆,τ ] = δ1,τvart [Xt+h∆ ]δ1,τ ,(35)

It is noted that although the constant term in equation (34), δ0,τ , does not appear inequation (35), it may aﬀect the volatility prediction by aﬀecting the estimate of δ1,τ .
But this impact is negligible. The results of the volatility prediction with δ0,τ excluded
preliminarily from equation (34) little change from those with δ0,τ included, which arepresented below. Equations (34) and (35) are common also in the next two stages.
In the second stage, we take only (iii) as given, and remove (i) and (ii). The physicaldistribution of changes in Xt with the same SV-Q and SV-E specifications of Σt is estimated using the in-sample data. These estimated parameters are fixed throughout theout-of-sample period.
In the last stage, we remove all of (i)–(iii). Thus, the underlying factors need to be
identified first. To avoid arbitrariness, we use the conventional level (levt ), slope (slot ),and curvature (curt ) factors, which are given by












 levt 

 0 0

 
 slo  =  −1 0
t 



curt









1   yt,0.5 




1 
  yt,2  .

−1 2 −1

yt,10

(36)

The rest of the process is the same as the second stage. It is noted that since Xt is
composed of the observed yields with τ = {0.5, 2, 10}, there is no residual in equation(34) for these maturities. Also, though the same SV-Q and SV-E specifications are used,significant parameters remaining in Σt are diﬀerent because diﬀerent factors are used.

31

Table 15 presents the in-sample RMSEs in basis point, where panels A and B present
the results based on equations (31) and (32), respectively. The rows labeled SV-Q and
SV-E show the original results, taken from Table 9, with the non-arbitrage imposed. NCi(i = 1, 2, 3) indicates “No Constraint”, beyond the linear regression, on the cross-sectional
relation with the number corresponding to the stage where the non-arbitrage is involved.
We limit our attention to h = 4 since the results for other forecasting horizons are basicallythe same: the imposition of non-arbitrage seems irrelevant to the forecasting horizon.
First, in panel A, the RMSEs for NC1 are little changed from the original ones for boththe SV-Q and SV-E specifications, indicating that once the factor and parameter values
are provided, there is virtually no diﬀerence between using the non-arbitrage condition andthe linear regression. Also, the result implies that although the SV-Q and SV-E models
are outside a class of aﬃne models, the cross-sectional relation between the factors andyields is close to linear. The RMSEs for NC2 are somewhat changed from the original ones,
however, the diﬀerences do not seem large. The result is attributed to similar estimates ofΣt between with and without the non-arbitrage, indicating that Σt is estimated mostly bytime-series properties of the data even though the non-arbitrage is imposed. The impact
of non-arbitrage, if any, seems neutral for the volatility forecasting performance: it isfavorable for SV-Q but not for SV-E. The RMSEs for NC3 are not much changed althoughsome reduction is observed for the SV-E specification. These results together indicate thatthe key to the volatility forecasting performance is neither the imposition of non-arbitrage
nor the choice of factors associated with this imposition, but the specification of Σt .
Next, looking at panel B, we notice that there is no diﬀerence at any stage, with aminor exception for NC3 of the SV-Q specification at τ = {0.5, 1}, where the RMSEs area bit increased. Thus, the volatility forecasting regression further decreases the impact of
non-arbitrage to a level where it is almost negligible. It is not surprising that these insample results are inherited into the out-of-sample results, so they are omitted for brevity.
In sum, the non-arbitrage condition has little eﬀect on the volatility prediction, as is thecase for the level prediction shown by Duﬀee (2011).

6

Concluding remarks

Andersen and Benzoni (2010, p.644) pointed out, after rejecting the aﬃne spanning conditions associated with the yield variance, that “further extensions to the term structure

32modeling framework are warranted.” This paper has shown one such extension. We improved the predictive power of term structure models with respect to the volatility. Basedon the three-factor Gaussian model, we made the instantaneous covariance matrix ofthe factors level-dependent. Specifically, the eigenvalues of this matrix are specified byquadratic (model SV-Q) and exponential (model SV-E) functions of the factors. The proposed models indeed predicted the volatility more accurately than the aﬃne modes. The
key to the success is not to identify the volatility factors directly from the yield curve,
but to specify the volatility functions with the yield curve factors and estimate themusing time-series properties of the data. The non-arbitrage constraint on cross-sectional
properties of the data is of little relevance to this success, however.
Compared with the time-series models, GARCH(1,1) and high frequency regressions,the volatility forecasting performance of the proposed models is low for the four-week
forecasting horizon especially at short maturities. However, by extending the horizons up
to 32 weeks, the relative performance of the proposed models increases. In fact, there aresome cases at middle to long maturities in which the performance is comparable or evensuperior to that of the time-series models in both the in-sample and out-of-sample periods.
We further considered combined forecasts, where a forecast of the proposed models wasadded to the high-frequency regression model as an explanatory variable. The predictive
power indeed increased when the forecast of the SV-Q model was included, which, however,
did not hold for the SV-E model.
Level-dependent volatility works from a statistical perspective. But less clear is whether
it does from an economic perspective. One approach to measure the economic significanceof level-dependent volatility, as developed by Fleming, Kirby, and Ostdiek (2001), is toconstruct bond portfolios based on model-implied first and second moments of returns andto examine their profitability. In this analysis, the predictability of the covariances, which
is not a main focus of this paper, has a crucial role. It may be the case that adequate models diﬀer between statistical and economic perspectives. Another question is, as argued
in Section 5.1, whether volatility-specific factors that do not appear in the risk-neutral
drift terms can further increase the predictive power to the time series of interest rates, orcontribute to the description of the cross-section of interest-rate option prices. Addressingthese questions is left for future research.

33

Appendix A: Impact of omitting the mean reversion on therealized volatility
We consider a CIR-type square-root process for each yield, yt,τ :√dyt,τ = (α + βyt,τ )dt + σ yt,τ dWt ,(37)

where Wt is Brownian motion in the physical measure. The variance of yt+h∆,τ conditioned
on yt,τ isσ 2 2βh∆ασ 2 βh∆(e− eβh∆ )yt,τ +(e− 1)2 .
β2β 2vart [yt+h∆,τ ; α, β, σ] =(38)

By setting β = 0, it is simplified asvart [yτ,t+h∆ ; σ] = σ 2 h∆yt,τ .

(39)Then, we compute the time series of % diﬀerence between the two conditional standard
deviations:
100

��

�

vart [yt+h∆,τ ; σ]−1vart [yt+h∆,τ ; α, β, σ].

We evaluate the impact of omitting the mean reversion by mean and standard deviationof the above series. The parameters are estimated using the in-sample data on eachyield. The quasi-maximum likelihood method is employed, where the conditional mean
and variance are substituted into the normal density function.
In panel A of Table A, the parameter estimates (standard errors) are presented. The
speed of mean reversion tends to increase the longer the maturity, as the estimates of β
become more negative. However, none of them is significant. In panel B of Table A, themeans (standard deviations) of the time series of % diﬀerence are presented: in computing
vart [yt+h∆,τ ; σ], we do not re-estimate σ but simply use the estimates shown in panel A.
For h = 32, the mean diﬀerence reaches nearly 15% at τ = 10. However, due to theinsignificant estimates of β, the actual impact of omitting the mean reversion is not as
serious as it appears.

34

Appendix B: An approximation method of conditional moments and its application to the pricing of bonds
B1. Outline of the method
The method is originally developed by Shoji (2002) and applied to the pricing of bonds
by Takamizawa and Shoji (2009). The method generally allows for the computation of up
to n-th conditional moments, if they exist, for a d-dimensional diﬀusion process. To ease
the explanation, we limit our attention to the case of (n, d) = (2, 2), i.e., the conditional
first and second moments of a two-dimensional diﬀusion process. As seen below, n can beconsidered as the order of approximation.
Let Xt = (xt,1 xt,2 )� be a two-dimensional diﬀusion process, which evolves according
to the following stochastic diﬀerential equation (SDE):
dxt,i = fi (Xt )dt + ξi (Xt )� dWt(i = 1, 2) ,(40)

where Wt is two-dimensional Brownian motion, and the drift and diﬀusion functions, fi
and ξi (i = 1, 2), satisfy certain technical conditions for the solution to equation (40) toexist for an arbitrary X0 .
Let Ψs,t be a vector consisting of the first and second moments of an increment of Xt
conditioned on time s < t:�Ψ�s,t = Es xt,1 − xs,1

xt,2 − xs,2

(xt,1 − xs,1 )2(xt,2 − xs,2 )2

�(xt,1 − xs,1 )(xt,2 − xs,2 ) .

The goal is to obtain an approximation of Ψs,t , which will turn out to be the solution to
an ODE.
By integrating equation (40) and taking the conditional expectation,Es [xt,i − xs,i ] = Es

�� t
s

�

fi (Xu )du .

(41)

By applying the Taylor expansion to fi (Xu ) around Xs up to the second orderfi (Xu ) = fi (Xs )
(1,0)

+fi

(0,1)

(Xs )(xu,1 − xs,1 ) + fi

1 (2,0)(Xs )(xu,2 − xs,2 ) + fi (Xs )(xu,1 − xs,1 )2
2

1 (0,2)(1,1)+ fi (Xs )(xu,2 − xs,2 )2 + fi (Xs )(xu,1 − xs,1 )(xu,2 − xs,2 ) + ei ,
2

35

(42)

k+l

∂
f
where f (k,l) = ∂x
By substituting equation (42) intok ∂xl , and ei is a residual term.
12equation (41) and expressing the resulting equation in a vector formEs [xt,i − xs,i ] = fi (t − s)�
1 (2,0)
(1,0)
(0,1)
+ fi
fi
f
2 i

1 (0,2)
f
2 i

(1,1)
fi

�� t
s

Ψs,u du + Ri ,(43)

where Xs is omitted for notational convenience, and Ri = Es [ei ].
Next, by applying the Ito formula to (xt,1 − xs,1 )2 and taking the conditional expectation,
2

Es [(xt,1 − xs,1 ) ] = Es

�� t
s

�

{2f1 (Xu )(xu,1 − xs,1 ) + g11 (Xu )}du ,(44)

where g11 = ξ1� ξ1 . By applying the Taylor expansion to f1 (Xu ) and g11 (Xu ) around Xs
up to the first and second orders, respectively, the integrand of equation (44) becomes
2f1 (Xu )(xu,1 − xs,1 ) + g11 (Xu )
(1,0)

(0,1)

= g11 (Xs ) + {2f1 (Xs ) + g11 (Xs )}(xu,1 − xs,1 ) + g11 (Xs )(xu,2 − xs,2 )
1 (2,0)1 (0,2)(1,0)+{2f1 (Xs ) + g11 (Xs )}(xu,1 − xs,1 )2 + g11 (Xs )(xu,2 − xs,2 )2
2
2
(0,1)
(1,1)
+{2f1 (Xs ) + g11 (Xs )}(xu,1 − xs,1 )(xu,2 − xs,2 ) + e11 ,(45)

where g (k,l) is defined analogously with f (k,l) , and e11 is a residual term. By substituting
equation (45) into equation (44),Es [(xt,1 − xs,1 )2 ] = g11 (t − s)
+
×

�

(1,0)
2f1 + g11

� t
s

(0,1)
g11

1 (2,0)
(1,0)
2f1
+ g11
2

1 (0,2)
g
2 11

(0,1)
(1,1)
2f1
+ g11

Ψs,u du + R11 ,

�

(46)

where R11 = Es [e11 ]. A similar manipulation is applied to Es [(xt,2 − xs,2 )2 ] and Es [(xt,1 −xs,1 )(xt,2 − xs,2 )]. Expressing the resulting equations together in a vector form leads toΨs,t = A(Xs )
where

A =



(1,0)
f1


(1,0)

f2



(1,0)
 2f1 + g11


(1,0)

g22


(1,0)
f2 + g12

(0,1)
f1
(0,1)
f2
(0,1)
g11
(0,1)
2f2 + g22
(0,1)
f1 + g12� ts

Ψs,u du + b(Xs )(t − s) + R ,1 (2,0)2 f1
1 (2,0)
2 f2
(1,0)
(2,0)
2f1
+ 12 g11
1 (2,0)
2 g22
(1,0)
(2,0)
f2
+ 12 g1236

1 (0,2)2 f1
1 (0,2)2 f2
1 (0,2)
2 g11
(0,1)
(0,2)
2f2
+ 12 g22
(0,1)
(0,2)
f1
+ 12 g12

(47)


(1,1)
f1


(1,1)

f2



(0,1)
(1,1)
 ,
2f1
+ g11


(1,0)
(1,1)

2f2
+ g22


(1,0)
(0,1)
(1,1)
f1
+ f2
+ g12

b = (f1f2

g11

R = (R1

R2

R11

g22

g12 )� ,R22

R12 )� .

Equation (47) can be solved asΨs,t =� tseA(Xs )(t−u) b(Xs )du + R̂ .

(48)

If, in addition, A is invertible, we obtain
Ψs,t = A−1 (Xs ){eA(Xs )(t−s) − I}b(Xs ) + R̂ .

(49)

It is noted that equations (47)–(49) hold for any (n, d) with modification to A(Xs ) andb(Xs ). In general, Ψs,t consists of

�n+d�
n

− 1 = (n + d)!/(n!d!) − 1 elements when up to n-th

conditional moments for a d-dimensional diﬀusion process are computed. Correspondingly,
up to n-th derivatives of fi and gij (i, j = 1, ..., d) are taken to compute the elementsof A(Xs ). Omitting the residual vector, R or R̂, leads to the approximation formula.
According to Shoji (2002), both R and R̂ have order of O((t − s)(n+3)/2 ). Thus, n can beconsidered as the order of approximation. In pricing bonds, we consider n = 3.
It is also noted that R contains the conditional expectation of derivatives of fi higher
than the first order and derivatives of gij higher than the second order. Then, if fi andgij are linear and quadratic in Xs , respectively, there is no residual term. In other words,
the conditional moments computed by the formula are exact. Even in such a case, the useof this formula may be beneficial when the derivation of closed-form conditional moments
is demanding.

B2. Application of the method
To apply the approximation method to the pricing of bonds, define�zs,t = exp −� tsr(Xu )du

�

,

(50)

and the price of a discount bond at time t maturing at time T is equal to the conditional
first moment of zt,T under the risk-neutral measure. This (actually EtQ [zt,T − zt,t ]) is
computed as one of the elements of the moment vector, Ψt (T ). Specifically, we firstextend a state vector as X̂t = (Xt� zs,t )� , where Xt is a d-dimensional diﬀusion process and37

zs,t is treated as the (d + 1)-th process. By the Ito formula,dzs,t = −r(Xt )zs,t dt ,zs,s = 1 ,(51)

and we have
fd+1 (X̂t ) = −r(Xt )zs,t ,(52)

gi d+1 (X̂t ) = 0(53)(i = 1, ..., d + 1) .

Then, the elements of A(X̂s ) can be readily computed by taking appropriate derivativesof fi (the risk-neutral drift functions here) and gij (i, j = 1, ..., d + 1). The accuracy of theapproximation to EtQ [zt,T ] is investigated in the next appendix.
The approximation method is also applied to the computation of conditional first andsecond moments of a model-implied yield, Es [Ỹ (Xt , τ )] and Es [Ỹ (Xt , τ )2 ]. Similar to thecase in zs,t , we first extend a state vector as X̂t = (Xt�Ỹ (Xt , τ ))� , and then derive theSDE for Ỹ (Xt , τ ):
�dỸ (Xt , τ ) = fd+1 (Xt , τ )dt + ξd+1
(Xt , τ )dWt ,

(54)

where
fd+1 (Xt , τ ) =�∂ Ỹ (Xt , τ )1∂ 2 Ỹ (Xt , τ )µt + trΣt�
∂Xt
2∂Xt ∂Xt�ξd+1 (Xt , τ ) = Σ0.5�
t∂ Ỹ (Xt , τ )
,
∂Xt

�

,

(55)
(56)

and where Wt is d-dimensional Brownian motion, and µt and Σt are the physical drift
vector and the instantaneous covariance matrix of dXt , respectively. gi d+1 is obtained byξi� ξd+1 (i = 1, ..., d + 1).
In computing A(X̂s ), the derivatives of fd+1 and gi d+1 are required. But fd+1 and gi d+1
already contain the derivatives of Ỹ (Xt , τ ) up to the second order, which complicates thecalculation. To avoid the tedious calculation, the derivatives of Ỹ (Xt , τ ) higher than thesecond order are omitted, while the derivatives of µt and Σt are taken as many times
as necessary. The accuracy of the approximation to Et [Ỹ (Xt+h∆ , τ )n ] (n = 1, 2) is also
investigated in the next appendix.

38

Appendix C: Accuracy of the approximation
The purpose of this appendix is to let the cost of using the approximation be known. By
construction of the method, the accuracy becomes worse the longer the time interval, t−s.
Here, the interval is up to ten years for pricing bonds, which raises concern with the application of this method. To check the accuracy of the approximation, we consider two caseswith and without a closed-form solution for bond prices. This is aimed at examining how
close the approximate solution is to the closed-form and numerical solutions, respectively.

C1. Comparison to the closed-form solution
We treat A0 (3) as the true model. We divide the analysis into three steps according to thedegrees of approximation involved. Let Θ0 be the parameter vector of the A0 (3) model,
the elements of which are set at the estimates presented in Table 4 (also re-exhibited inTable C2). The extracted state vector can then be expressed as X(Ytp ; Θ0 ).
In the first step, we examine the impact of the approximation on the pricing of bonds
alone. Specifically, both Θ0 and X(Ytp ; Θ0 ) are given as input for the approximation
method. Then, we compare
Y (X(Ytp ; Θ0 ), τ ; Θ0 ) v.s.

Ỹ (X(Ytp ; Θ0 ), τ ; Θ0 )(τ = {0.5, 1, 2, 3, 5, 10}) ,where Y (·) and Ỹ (·) stand for the closed-form and approximate functions, respectively.
In the second step, we examine the impact of the approximation on the extractionof state variables as well as on the pricing of bonds. Here, only Θ0 is given. Using theapproximation, the state vector is first backed out, which is denoted as X̃(Ytp ; Θ0 ), andthe rest of the yields are computed. Then, we compare:X(Ytp ; Θ0 )v.s. X̃(Ytp ; Θ0 ) ,Y (X(Ytp ; Θ0 ), τ ; Θ0 ) v.s. Ỹ (X̃(Ytp ; Θ0 ), τ ; Θ0 ) (τ = {1, 3, 5}) .
Note that at τ = {0.5, 2, 10}, both Y (Xt , τ ; Θ0 ) and Ỹ (X̃t , τ ; Θ0 ) are equal to the observed
yields by construction of the inversion method.
In the last step, we examine the impact of the approximation on the estimation ofmodel parameters as well as on the pricing of bonds and the extraction of state variables.
Here, no prior information on the true values of the parameter and state vectors is given.
Using the approximation, the parameter vector of the A0 (3) model is first estimated;39

denote it as Θ̃0 . Next, the state vector is backed out; denote it as X̃(Ytp ; Θ̃0 ). Finally, therest of the yields are computed. Then, we compare:Θ0v.s. Θ̃0

X(Ytp ; Θ0 )v.s. X̃(Ytp ; Θ̃0 ) ,Y (X(Ytp ; Θ0 ), τ ; Θ0 ) v.s. Ỹ (X̃(Ytp ; Θ̃0 ), τ ; Θ̃0 ) (τ = {1, 3, 5}) .
It is noted that the accuracy in the third step, which is a more realistic setup, is not
examined by Takamizawa and Shoji (2009).
Apart from the parameter vector, the key input for these comparisons is Ytp . We use
the actual data on Ytp . To condense the analysis, we pick up nine observations from theentire sample in the following way. First, Ytp is transformed to the conventional level(levt ), slope (slot ), and curvature (curt ) factors by equation (36). Then, we choose threedates in which levt takes the minimum, median, or maximum value. Likewise, the threedates are chosen for each of the other proxies, leading to nine dates in total. In this way,
the accuracy of the approximation is evaluated at not only the usual times but also theunusual times.
Table C1 presents the diﬀerences, expressed in basis points, between the approximate
and true yields or state variables. Panel A presents the results for the first step comparison, where the true values of the parameter and state vectors are given as input for theapproximation. For maturities of up to five years, the approximation errors are negligibly
small at all states. Even for the ten-year maturity, the error exceeds 2 bps only at the
maximum-spread state. The values of the curvature, on the other hand, little aﬀects theaccuracy.
Panel B of Table C1 presents the results for the second-step comparison, where only
the true value of the parameter vector is given. A systematic pattern is found in theapproximation errors for the state vector. Specifically, both r and θ are undervalued,
which is compensated by the overvaluation of �. The diﬃculty of the approximation arises
again at the maximum-spread state. On the other hand, the approximation errors for theremaining yields are small.
Panel C of Table C1 presents the results for the third-step comparison, where noprior information is given. Before discussing the results, we first look at Table C2, where
the parameter estimates (standard errors) obtained by the approximation method arepresented. Overall, the estimates do not seem to diﬀer much from the corresponding true40

values. Except for θ̄, the diﬀerences in absolute value are smaller than the standard errors.
Now, looking back at panel C of Table C1, we find, nevertheless, that the magnitude of
the approximation errors generally increases. This is the reality. Here, an error pattern isless clear, but it is seen that r tends to be overvalued and θ tends to undervalued. Also,
the diﬃculty of the approximation is not limited to the maximum-spread state. In fact,
the largest errors in absolute value for r and θ appear at the minimum-curvature state andthe minimum-level state, respectively. The remaining yields are accurately computed.

C2. Comparison to the numerical solution
We employ the Monte Carlo (MC) method to evaluate the accuracy. The models we use
are SV-Q and SV-E. Let Θi (i = {Q, E}) be the true parameter vectors, the elements ofwhich are set at the estimates presented in Tables 5 and 6, respectively. Further, the statevector is extracted by the approximation method, not by the MC method; denote it asX̃(Ytp ; Θi ). Then, we compare
Y (X̃(Ytp ; Θi ), τ ; Θi ) v.s.

Ỹ (X̃(Ytp ; Θi ), τ ; Θi )(τ = {0.5, 1, 2, 3, 5, 10}) .

In the MC simulations, {Xs }t+τis generated from (11) (the risk-neutral distribution),twhere dt is replaced by ∆t. We set ∆t = 1/1, 000, an interval corresponding roughly
to four observations per day. The number of repetition is set at 10,000 with antithetic
variates.
Table C3 presents the diﬀerences, expressed in basis points, between the approximate
and MC yields. Generally, the error pattern is similar to that for the first-step comparison
to the closed-form solution. For maturities of up to five years, the approximation errors
are within 1 bp at all states for both models. For the ten-year yield, the approximation
error exceeds 2 bps only at the maximum-spread state. However, this comparison does nottake the approximation errors in both the parameter and state vectors into consideration.
In reality, therefore, the approximation errors for the resulting yields would be larger, as
seen in the third-step comparison to the closed-form solution.

41

C3. Accuracy to conditional first and second moments of a model-implied
yield
Based on the SV-Q and SV-E models, we again employ the MC method to evaluate theaccuracy at the nine states. In the MC simulations, {Xs }t+h∆
is generated from (10)t(the physical distribution), where Σ is replaced by Σt and dt by ∆t (= 1/1, 000). For
simplicity, once Xt+h∆ is obtained, a model-implied τ -maturity yield (and its squared
value) is computed by the approximation method. This procedure is repeated 10,000
times with antithetic variates. Then, we compare
v.s. Etap [Ỹ (Xt+h∆ , τ )] ,Et [Ỹ (Xt+h∆ , τ )]�vart [Ỹ (Xt+h∆ , τ )]/h∆ v.s.

�

varap
t [Ỹ (Xt+h∆ , τ )]/h∆ ,where the superscript “ap” is for clarifying the approximation. We set h = 32 as this is
the longest horizon we consider and thus the accuracy is the worst.
Table C4 presents the diﬀerences, expressed in basis points, in the conditional mean.
The errors in absolute value are all within 1 bp for both models. Table C5 presents thediﬀerences, expressed in basis points, in the annualized conditional standard deviation.
Although the accuracy to the second moment becomes worse than that to the first moment,
the errors in absolute value are all within 2 bps.
It is noted that the more rigorous analysis of the accuracy requires to compare
Et [Y (Xt+h∆ , τ )n ]v.s.

Etap [Ỹ (Xt+h∆ , τ )n ](n = 1, 2) .

That is, no approximation is involved in the former in computing the τ -maturity yield attime t + h∆: but this analysis is computationally very demanding. In reality, therefore,
the approximation errors would be larger than presented in Tables C4 and C5.

42

Appendix D: Computation of the conditional variance for the
GARCH (1,1) model
The variance of yt+∆,τ conditioned on time t is simply Vt+∆,τ , which is observed at timet. The variance of yt+k∆,τ (k = 2, ..., h) conditioned on time t is computed iteratively
as follows. In equation (27), by substituting t + (k − 1)∆ for t and taking the varianceconditioned on time t,vart [yt+k∆,τ ] = vart [α + βyt+(k−1)∆,τ +�

Vt+k∆,τ zt+k∆,τ ]�

= β 2 vart [yt+(k−1)∆,τ ] + vart [ Vt+k∆,τ zt+k∆,τ ]+2β covt [yt+(k−1)∆,τ ,�

Vt+k∆,τ zt+k∆,τ ]= β 2 vart [yt+(k−1)∆,τ ] + Et [Vt+k∆ ](k = 2, ..., h) .

(57)

On the other hand, in equation (28), by substituting t + (k − 1)∆ for t and taking the
expectation conditioned on time t,
2
Et [Vt+k∆ ] = Et [ω0 + ω1 Vt+(k−1)∆,τ zt+(k−1)∆,τ+ ω2 Vt+(k−1)∆,τ ]= ω0 + (ω1 + ω2 )Et [Vt+(k−1)∆,τ ](k = 2, ..., h) .

(58)

Then, vart [yt+h∆,τ ] is obtained by iteratively solving equations (57) and (58) starting fromvart [yt+∆,τ ] = Vt+∆,τ .

43

ReferencesAı̈t-Sahalia, Y., 1996, Testing continuous-time models of the spot interest rate, Review of
Financial Studies 9, 385-426.
Andersen, T. G., and L. Benzoni, 2010, Do bonds span volatility risk in the U.S. Treasury
market? A specification test for aﬃne term structure models, Journal of Finance 65,603-653.
Andersen, T. G., T. Bollerslev, F. X. Diebold, and P. Labys, 2003, Modeling and forecasting realized volatility, Econometrica 71, 579-625.
Andersen, T. G., and J. Lund, 1997a, Estimating continuous-time stochastic volatility
models of the short-term interest rate, Journal of Econometrics 77, 343-377.
Andersen, T. G., and J. Lund, 1997b, Stochastic volatility and mean drift in the short
rate diﬀusion: Sources of steepness, level and curvature in the yield curve, Working paper.
Bakshi, G., C. Cao, and Z. Chen, 1997, Empirical performance of alternative option pricing
models, Journal of Finance 52, 2003-2049.
Balduzzi, P., S. R. Das, and S. Foresi, 1998, The central tendency: A second factor inbond yields, Review of Economics and Statistics 80, 62-72.
Balduzzi, P., S. R. Das, S. Foresi, and R. Sundaram, 1996, A simple approach to threefactor aﬃne term structure models, Journal of Fixed Income 6, 43-53.
Ball, C. A., and W. N. Torous, 1999, The stochastic volatility of short-term interest rates:
Some international evidence, Journal of Finance 54, 2339-2359.
Bates, D. S., 1996, Jumps and stochastic volatility: Exchange rate processes implicit in
Deutsche mark options, Review of Financial Studies 9, 69-107.
Bates, D. S., 2000, Post-’87 crash fears in the S&P 500 futures option market, Journal of
Econometrics 94, 181-238.
Bekaert, G., and R. J. Hodrick, 2001, Expectations hypotheses tests, Journal of Finance
56, 1357-1394.
Bikbov, R., and M. Chernov, 2011, Yield curve and volatility: Lessons from Eurodollar
futures and options, Journal of Financial Econometrics 9, 66-105.

44

Brown, R. H., and S. M. Schaefer, 1994, Interest rate volatility and the shape of the term
structure, Philosophical Transactions: Physical Sciences and Engineering 347, 563-576.
Campbell, J. Y., and R. J. Shiller, 1991, Yield spreads and interest rate movements: A
bird’s eye view, Review of Economic Studies 58, 495-514.
Chan, K. C., G. A. Karolyi, F. A. Longstaﬀ, and A. B. Sanders, 1992, An empiricalcomparison of alternative models of the short-term interest rate, Journal of Finance 47,
1209-1227.
Chen, R., and L. Scott, 1993, Maximum likelihood estimation for a multifactor equilibrium
model of the term structure of interest rates, Journal of Fixed Income 3, 14-31.
Christiansen, C., and J. Lund, 2005, Revisiting the shape of the yield curve: The eﬀect ofinterest rate volatility, Working paper, University of Aarhus.
Collin-Dufresne, P., and R. S. Goldstein, 2002, Do bonds span the fixed income markets?
Theory and evidence for unspanned stochastic volatility, Journal of Finance 57, 1685-1730.
Collin-Dufresne, P., R. S. Goldstein, and C. S. Jones, 2008, Identification of maximal aﬃneterm structure models, Journal of Finance 63, 743-795.
Collin-Dufresne, P., R. S. Goldstein, and C. S. Jones, 2009, Can interest rate volatility
be extracted from the cross section of bond yields? Journal of Financial Economics 94,
47-66.
Cox, J. C., J. E. Ingersoll, Jr., and S. A. Ross, 1985, A theory of the term structure ofinterest rates, Econometrica, 53, 385-407.
Dai, Q., and K. J. Singleton, 2000, Specification analysis of aﬃne term structure models,
Journal of Finance 55, 1943-1978.
Dai, Q., and K. J. Singleton, 2003, Term structure dynamics in theory and reality, Review
of Financial Studies 16, 631-678.
Dai, Q., K. J. Singleton, and W. Yang, 2007, Regime shifts in a dynamic term structure
model of U.S. Treasury bond yields, Review of Financial Studies 20, 1669-1706.
Duﬀee, G. R., 2002, Term premia and interest rate forecasts in aﬃne models, Journal of
Finance 57, 405-443.
Duﬀee, G. R., 2011, Forecasting with the term structure: The role of no-arbitrage restrictions, Working Paper, Johns Hopkins University.

45

Durham, G. B., 2003, Likelihood-based specification analysis of continuous-time modelsof the short-term interest rate, Journal of Financial Economics 70, 463-487.
Fama E. F., and R. R. Bliss, 1987, The information in long-maturity forward rates, American Economic Review 77, 680-692.
Fan, R., A. Gupta, and P. Ritchken, 2003, Hedging in the possible presence of unspanned
stochastic volatility: Evidence from swaption markets, Journal of Finance 58, 2219-2248.
Fleming, J., C. Kirby, and B. Ostdiek, 2001, The economic value of volatility timing,
Journal of Finance 56, 329-352.
Gallant, A. R., and G. Tauchen, 1998, Reprojecting partially observed systems with application to interest rate diﬀusions, Journal of the American Statistical Association 93,
10-24.
Ghysels, E., P. Santa-Clara, and R. Valkanov, 2005, There is a risk-return trade-oﬀ after
all, Journal of Financial Economics 76, 509-548.
Ghysels, E., P. Santa-Clara, and R. Valkanov, 2006, Predicting volatility: Getting themost out of return data sampled at diﬀerent frequencies, Journal of Econometrics 131,
59-95.
Han, B., 2007, Stochastic volatilities and correlations of bond yields, Journal of Finance
62, 1491-1524.
Hansen, P. R., and A. Lunde, 2005, A forecast comparison of volatility models: Does
anything beat a GARCH(1,1)? Journal of Applied Econometrics 20, 873-889.
Heidari, M., and L. Wu, 2003, Are interest rate derivatives spanned by the term structureof interest rates? Journal of Fixed Income 13, 75-86.
Heidari, M., and L. Wu, 2009, A Joint framework for consistently pricing interest rates
and interest rate derivatives, Journal of Financial and Quantitative Analysis 44, 517-550.
Jacobs, K., and L. Karoui, 2009, Conditional volatility in aﬃne term-structure models:
Evidence from Treasury and swap markets, Journal of Financial Economics 91, 288-318.
Jarrow, R. A., H. Li, and F. Zhao, 2007, Interest rate caps smile too! But can the LIBOR
market models capture the smile? Journal of Finance 62, 345-382.
Joslin, S., 2006, Can unspanned stochastic volatility models explain the cross section ofbond volatilities? Working Paper, MIT Sloan School of Management.

46

Joslin, S., 2010, Pricing and hedging volatility in fixed income markets, Working Paper,
MIT Sloan School of Management.
Joslin, S., K. J. Singleton, and H. Zhu, 2011, A new perspective on Gaussian dynamic
term structure models, Review of Financial Studies 24, 926-970.
Li, H., and F. Zhao, 2006, Unspanned stochastic volatility: Evidence from hedging interest
rate derivatives, Journal of Finance 61, 341-378.
Litterman, R., J. Scheinkman, and L. Weiss, 1991, Volatility and the yield curve, Journalof Fixed Income 1, 49-53.
Longstaﬀ, F. A., P. Santa-Clara, and E. S. Schwartz, 2001, The relative valuation of caps
and swaptions: Theory and empirical evidence, Journal of Finance 56, 2067-2109.
Newey, W. K., and K. D. West, 1987, A simple, positive semi-definite, heteroskedasticity
and autocorrelation consistent covariance matrix, Econometrica 55, 703-708.
Pearson, N. D., and T. Sun, 1994, Exploiting the conditional density in estimating theterm structure: An application to the Cox, Ingersoll, and Ross model, Journal of Finance
49, 1279-1304.
Phoa, W., 1997, Can you drive market volatility forecasts from the observed yield curve
convexity bias? Journal of Fixed Income 7, 43-54.
Shoji, I., 2002, Approximation of conditional moments of diﬀusion processes, International
Journal of Computational and Numerical Analysis and Applications 1, 163-190.
Takamizawa, H., and I. Shoji, 2009, Modeling the term structure of interest rates withgeneral diﬀusion processes: A moment approximation approach, Journal of Economic
Dynamics and Control 33, 65-77.
Thompson, S., 2008, Identifying term structure volatility from the LIBOR-swap curve,
Review of Financial Studies 21, 819-854.

47

Maturity
h=4
Y-PART
Y-ALL
AR-RV
HAR-RV

0.5

1

2

3

5

10

0.34
0.36
0.39
0.43

0.30
0.33
0.35
0.39

0.17
0.23
0.22
0.24

0.13
0.16
0.13
0.14

0.13
0.15
0.14
0.15

0.07
0.08
0.12
0.12

h=8
Y-PART
Y-ALL
AR-RV
HAR-RV

0.39
0.42
0.40
0.45

0.35
0.40
0.40
0.42

0.23
0.29
0.21
0.23

0.18
0.21
0.11
0.12

0.17
0.21
0.13
0.14

0.10
0.11
0.07
0.09

h = 16
Y-PART
Y-ALL
AR-RV
HAR-RV

0.42
0.45
0.40
0.44

0.39
0.44
0.34
0.40

0.28
0.35
0.15
0.22

0.24
0.27
0.06
0.13

0.22
0.25
0.06
0.12

0.16
0.17
0.00
0.07

h = 32
Y-PART
Y-ALL
AR-RV
HAR-RV

0.40
0.43
0.33
0.38

0.40
0.44
0.33
0.37

0.31
0.40
0.17
0.20

0.31
0.35
0.11
0.12

0.27
0.32
0.06
0.07

0.28
0.29
0.03
0.03

Table 1: Adjusted R2 coeﬃcients for the volatility forecasting regression√
The regression models are given by equations (3)–(6), where f (x) = x and h = {4, 8, 16, 32}.
The best and the second best results in each column of each panel are displayed in bold
and italic, respectively. The in-sample data from January 4, 1991 to April 9, 2003 areused for the estimation.

48

Maturity
h=4
Y-PART
Y-ALL
AR-RV
HAR-RV

0.5

1

2

3

5

10

45.0
40.4
31.5
30.7

41.9
38.3
28.6
27.3

42.9
42.2
32.2
30.5

41.4
39.8
34.0
32.2

39.9
39.7
34.3
32.4

40.6
39.7
35.5
35.7

h=8
Y-PART
Y-ALL
AR-RV
HAR-RV

42.4
36.6
30.7
29.3

39.3
34.7
25.0
24.4

39.8
38.0
29.1
28.1

38.6
35.3
32.2
30.6

37.4
35.5
32.0
31.2

37.7
36.0
35.2
35.3

h = 16
Y-PART
Y-ALL
AR-RV
HAR-RV

40.3
35.1
29.4
27.7

37.5
33.8
23.8
23.2

37.7
36.3
30.1
28.4

37.0
33.7
33.6
30.6

36.5
34.1
34.4
31.7

36.3
34.1
38.3
34.9

h = 32
Y-PART
Y-ALL
AR-RV
HAR-RV

37.3
35.1
25.6
25.6

35.1
34.8
23.4
22.3

36.6
35.9
29.9
28.4

36.7
34.4
31.5
30.6

35.4
33.2
32.7
32.0

34.5
33.0
31.5
31.4

Table 2: Out-of-sample RMSEs in the fixed-parameter approach
Out-of-sample forecasts are generated from the regression equations (3)–(6), where f (x) =√x and h = {4, 8, 16, 32}. The parameters are fixed at the in-sample estimates throughout
the out-of-sample period from April 16, 2003 to May 27, 2009. The best and the second
best results in each column of each panel are displayed in bold and italic, respectively.

49

Maturity
h=4
Y-PART
Y-ALL
AR-RV
HAR-RV

0.5

1

2

3

5

10

34.3
31.4
32.0
31.1

32.4
30.6
29.1
27.7

34.8
34.0
31.7
30.1

33.9
32.9
32.3
30.9

32.3
30.7
31.8
29.9

33.5
28.8
33.1
32.2

h=8
Y-PART
Y-ALL
AR-RV
HAR-RV

32.1
29.3
31.0
29.8

29.3
27.7
25.3
24.9

31.3
30.6
27.9
27.5

30.4
29.2
29.4
28.8

28.4
26.3
28.2
27.2

29.2
23.7
31.7
30.0

h = 16
Y-PART
Y-ALL
AR-RV
HAR-RV

31.2
28.2
29.8
27.5

27.7
26.0
24.3
23.6

29.8
29.4
29.0
27.8

29.0
28.3
30.6
28.9

27.4
25.9
29.5
26.7

28.0
23.3
33.4
28.6

h = 32
Y-PART
Y-ALL
AR-RV
HAR-RV

29.4
27.6
25.9
25.9

25.8
24.9
23.8
22.2

28.0
27.6
28.5
26.3

28.0
27.2
29.3
27.5

27.4
26.0
27.6
26.0

27.0
23.5
27.3
27.0

Table 3: Out-of-sample RMSEs in the varying-parameter approach
Out-of-sample forecasts are generated from the regression equations (3)–(6), where f (x) =√x and h = {4, 8, 16, 32}. The parameters are re-estimated in a rolling-window fashion
with the sample size fixed at the same as the in-sample data. The best and the second
best results in each column of each panel are displayed in bold and italic, respectively.
The out-of-sample period is from April 16, 2003 to May 27, 2009.

50

Parameter / Index
κi
θ̄

i=1
0.947 (0.087)

i=2
0.017 (0.002)
0.146 (0.008)

σ1i × 104
σ2i × 104
σ3i × 104

0.722

0.095
1.220

λ0i

0.000

0.030

λ11i
λ12i
λ13i

0.000
0.000
0.000

(0.030)

ς × 102LogL

i=30.957 (0.090)−0.9850.000
8.036(0.115)(0.019)−0.013(0.008)−0.195−0.5140.000(0.037)(0.280)

0.000
0.000
0.000

0.061
22263

(0.001)(0.035)(0.063)(0.438)

Table 4: Parameter estimates (standard errors) for the A0 (3) model
The risk-neutral and physical distributions of the instantaneous change in Xt for the A0 (3)model are given by (7) and (10), respectively. ς is the standard deviation of the measurement errors, ut,τ = yt,τ − Ỹ (Xt , τ ) with τ = {1, 3, 5}. LogL stands for the maximum
log-likelihood value. Standard errors in parenthesis are computed by the outer productof the gradient vector of the log-likelihood function. The in-sample data from January 4,
1991 to April 9, 2003 are used for the estimation.

51

Parameter / Index
κi
θ̄

i=1
1.034 (0.073)

i=2
0.019 (0.002)
0.137 (0.007)

i=30.917 (0.063)

sin ϕPi

0.000

0.126

(0.010)0.000

ci × 103

1e−5

0.111

(0.008)0.149(0.107)

m1i

0.000

0.022

(0.002)0.260(0.034)0.631(0.021)−0.705(0.042)0.097(0.054)

1
sin ϕQ
i

0.250

m2i

(0.047)

0.000

0.000

2
sin ϕQ
i

0.000

0.547

(0.105)−0.547

m3i

0.000

0.228

(0.061)

2.153(0.627)0.551(0.058)−0.616(0.087)(0.008)3
sin ϕQ
i

0.273

λ0i

0.000

0.020

(0.020)−0.011

λ11i
λ12i
λ13i

0.000
0.000
0.000

−0.150−0.3570.000(0.031)(0.301)

0.000
0.000
0.000

0.061
22415

(0.001)(0.134)

ς × 102
LogL

Table 5: Parameter estimates (standard errors) for the SV-Q model
The risk-neutral and physical distributions of the instantaneous change in Xt for the SV-Q
model are given by (7) and (10), respectively, with Σ replaced by Σt . Σt is decomposedas Σt = P Lt P � , where Lt is the diagonal eigenvalue matrix, and P is the orthogonal
eigenvector matrix parameterized in equation (14). The i-th diagonal element of Lt is
specified as Li (Xt ) = ci +Xt� Γi Xt (i = 1, 2, 3), where ci > 0 and Γi is a non-negative definite
matrix. Γi is also parameterized based on the spectral decomposition as Γi = Qi M i Qi�
(i = 1, 2, 3), where M i is the diagonal eigenvalue matrix with its elements satisfying
0 ≤ mi1 ≤ mi2 ≤ mi3 , and Qi is the orthogonal eigenvector matrix parameterized inequation (18). ς is the standard deviation of the measurement errors, ut,τ = yt,τ − Ỹ (Xt , τ )with τ = {1, 3, 5}. LogL stands for the maximum log-likelihood value. Standard errors inparenthesis are computed by the outer product of the gradient vector of the log-likelihood
function. The in-sample data from January 4, 1991 to April 9, 2003 are used for theestimation.

52

Parameter / Index
κi
θ̄

i=1
0.955 (0.070)

i=20.023 (0.003)
0.118 (0.006)

i=30.981 (0.078)

sin ϕPi

0.000

0.123

(0.013)0.000

ci−12.942(0.263)−8.693(0.214)−8.574(0.379)

γi1
γi2
γi3

10.538−8.828−13.489(2.235)(3.697)(4.665)35.1100.000
29.079(3.476)−17.267−10.551−9.665(2.431)(2.943)(2.733)

λ0i

0.000

0.036

(0.022)−0.014(0.008)

λ11i
λ12i
λ13i

0.000
0.000
0.000

−0.171−0.5950.000(0.038)(0.318)

0.000
0.000
0.000

0.061
22358

(0.001)

ς × 102
LogL(4.460)

Table 6: Parameter estimates (standard errors) for the SV-E model
The risk-neutral and physical distributions of the instantaneous change in Xt for the SV-E
model are given by (7) and (10), respectively, with Σ replaced by Σt . Σt is decomposed as
Σt = P Lt P � , where Lt is the diagonal eigenvalue matrix, and P is the orthogonal eigenvector matrix parameterized in equation (14). The i-th diagonal element of Lt is specified as
Li (Xt ) = exp{ci + γ i � Xt } (i = 1, 2, 3). ς is the standard deviation of the measurement errors, ut,τ = yt,τ − Ỹ (Xt , τ ) with τ = {1, 3, 5}. LogL stands for the maximum log-likelihood
value. Standard errors in parenthesis are computed by the outer product of the gradientvector of the log-likelihood function. The in-sample data from January 4, 1991 to April
9, 2003 are used for the estimation.

53

Maturity
h=4
RW
A0 (3)
SV-Q
SV-E

0.5

1

2

3

5

10

23.1
22.0
22.1
22.0

26.8
28.2
28.7
28.3

30.3
30.1
30.2
30.1

30.1
29.7
29.8
29.7

28.4
28.1
28.2
28.1

26.1
25.7
25.8
25.7

h=8
RW
A0 (3)
SV-Q
SV-E

37.8
35.9
36.2
35.9

42.9
44.0
44.8
44.1

47.7
47.4
47.6
47.4

46.9
46.1
46.3
46.1

43.6
42.5
42.7
42.5

37.8
36.9
37.0
36.9

h = 16RW
A0 (3)
SV-Q
SV-E

63.4
60.7
61.4
60.6

69.2
70.5
71.8
70.6

74.3
73.5
74.2
73.4

72.1
70.0
70.7
69.9

66.2
63.2
63.7
63.1

55.6
52.8
53.2
52.8

h = 32RW
A0 (3)
SV-Q
SV-E

107.0
102.7
104.3
102.4

109.9
109.0
111.6
108.8

110.8
106.1
108.3
105.7

105.4
98.4
100.6
97.8

96.1
86.9
89.0
86.3

82.4
73.3
74.8
72.9

Table 7: In-sample RMSEs for the level prediction
Root mean squared errors (RMSEs) for the in-sample prediction of the h-week ahead yields
are presented in basis points with h = {4, 8, 16, 32}. RW represents random walk, whose
forecast is the current yield for any h. The in-sample period is from January 4, 1991 toApril 9, 2003.

54

Maturity
h=4
RW
A0 (3)
SV-Q
SV-E

0.5

1

2

3

5

10

28.9
25.2
24.8
24.9

28.2
27.8
27.6
27.4

29.4
28.2
28.2
28.2

31.0
28.8
29.0
28.9

32.4
30.7
30.8
30.9

31.5
31.5
31.4
31.7

h=8
RW
A0 (3)
SV-Q
SV-E

44.5
35.9
35.0
35.1

42.7
37.2
36.9
36.4

45.1
42.1
42.2
42.0

46.5
42.2
42.5
42.4

47.3
43.7
44.0
44.1

45.5
45.4
45.3
45.9

h = 16RW
A0 (3)
SV-Q
SV-E

70.3
55.3
53.4
53.6

66.6
55.4
54.6
53.9

70.1
63.7
64.0
63.7

68.4
61.2
61.6
61.7

63.8
59.1
59.1
60.1

58.2
59.0
58.2
60.6

h = 32RW
A0 (3)
SV-Q
SV-E

112.3
91.1
88.0
88.0

104.8
85.0
83.6
82.6

102.6
88.8
89.4
89.1

93.8
82.9
83.5
84.6

81.0
78.2
77.3
81.4

69.4
77.2
73.8
81.8

Table 8: Out-of-sample RMSEs for the level prediction
Root mean squared errors (RMSEs) for the out-of-sample prediction of the h-week ahead
yields are presented in basis points with h = {4, 8, 16, 32}. RW represents random walk,whose forecast is the current yield for any h. The out-of-sample period is from April 16,
2003 to May 27, 2009.

55

Maturity
h=4
A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

0.5

1

2

3

5

10

39.2
36.1
33.7
32.0
33.1
29.8
28.7
27.8
27.1
25.6
25.6

39.8
37.0
35.0
33.8
31.4
34.1
32.5
31.9
29.1
29.1
28.8

34.8
33.4
32.8
32.8
31.1
33.1
31.1
31.5
30.5
29.8
29.5

34.6
34.1
33.9
33.6
32.3
33.9
32.4
32.3
31.9
31.8
31.6

28.9
29.4
28.3
28.4
27.0
28.7
27.5
27.1
26.6
26.4
26.1

26.6
28.0
26.9
27.8
25.0
26.4
26.3
25.8
24.7
24.9
24.7

h=8
A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

37.5
34.0
32.2
30.1
31.6
27.0
25.9
24.7
24.6
23.3
23.5

37.4
33.9
32.0
30.8
27.8
30.6
28.7
28.1
25.5
26.1
25.9

30.6
28.8
27.7
28.0
26.9
28.7
26.0
26.7
26.3
26.2
26.0

29.5
28.8
27.9
28.1
27.5
28.8
26.7
26.9
27.1
27.4
27.2

24.2
24.9
23.4
23.7
22.9
24.2
22.7
22.4
22.5
22.4
22.1

22.1
23.7
22.4
23.2
20.9
21.9
21.8
21.2
20.9
21.0
20.8

Table 9: In-sample RMSEs for the volatility prediction

56

Maturity
h = 16A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

0.5

1

2

3

5

10

37.8
33.1
34.1
31.1
32.0
24.2
23.7
22.2
22.8
21.6
21.6

36.2
31.4
31.4
29.7
26.6
27.3
25.5
25.2
23.9
24.4
24.4

26.9
24.9
24.1
24.0
24.0
24.7
21.9
22.6
23.4
23.1
23.4

24.8
24.3
22.5
22.9
23.4
24.0
21.5
21.8
22.9
23.0
23.4

20.5
21.6
19.6
20.1
19.9
20.4
19.0
18.5
19.3
19.1
19.0

18.6
20.7
18.7
19.8
17.7
17.5
17.7
17.1
17.4
17.2
17.3

h = 32A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

42.2
34.3
41.8
37.8
35.0
21.4
21.7
20.6
20.9
20.2
20.1

36.2
28.8
33.9
31.1
24.9
23.8
22.5
22.7
21.5
22.4
22.3

22.8
21.1
21.5
20.6
20.0
20.5
18.3
18.8
19.7
19.8
19.7

20.0
21.3
17.9
18.7
19.3
18.8
16.8
17.0
18.2
18.5
18.6

17.7
20.5
15.9
18.2
17.5
16.2
15.0
14.7
15.6
15.6
15.7

17.2
20.1
15.5
19.1
16.4
12.5
12.9
12.4
12.9
12.7
13.0

Table 9 (continued): In-sample RMSEs for the volatility prediction
Root mean squared errors (RMSEs) for the in-sample prediction of the h-week ahead
yield volatilities (annualized standard deviations) are presented in basis points with h ={4, 8, 16, 32}. For the term structure and GARCH models, the name alone indicates theresults based on equation (31), whereas the name “+ Reg” indicates the results based onequation (32). The best and the second best results in each column of each panel are
displayed in bold and italic, respectively. The in-sample period is from January 4, 1991
to April 9, 2003.

57

Maturity
h=4
A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

0.5

1

2

3

5

10

50.0
43.9
47.5
39.9
46.7
51.5
36.1
39.3
34.1
30.7
30.3

42.3
37.1
42.5
35.9
34.3
47.5
33.8
35.7
30.7
27.3
26.8

42.9
44.2
37.7
37.9
31.5
46.8
34.3
38.8
31.8
30.5
30.3

43.4
46.7
34.5
37.2
31.0
46.1
34.6
39.1
32.3
32.2
32.5

44.3
50.5
32.5
35.7
30.0
45.3
35.9
39.3
32.5
32.4
31.8

44.3
49.9
34.1
35.5
34.7
43.6
39.8
40.2
34.3
35.7
34.5

h=8
A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

48.8
42.0
46.8
39.0
47.2
49.4
34.6
37.9
33.6
29.3
29.6

40.8
34.6
40.5
34.1
31.9
45.8
31.0
33.5
28.5
24.4
24.3

40.4
42.1
33.5
35.2
27.7
44.7
30.6
36.1
28.6
28.1
28.4

41.1
45.2
30.5
34.9
27.9
44.2
30.9
36.5
29.9
30.6
31.6

42.3
49.4
28.8
33.5
27.6
42.9
32.6
36.8
31.1
31.2
30.9

42.2
48.6
30.8
32.6
33.8
40.9
37.0
37.5
33.8
35.3
34.0

Table 10: Out-of-sample RMSEs for the volatility prediction in the fixedparameter approach58

Maturity
h = 16A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

0.5

1

2

3

5

10

49.8
40.9
50.8
40.6
49.2
46.6
35.1
37.2
33.6
27.7
28.1

41.2
32.9
42.4
34.5
30.9
44.0
30.6
32.4
27.0
23.2
23.3

38.8
40.8
32.1
33.9
27.3
43.0
29.5
34.5
28.7
28.4
29.0

39.5
44.5
29.4
34.0
28.2
42.7
29.9
34.9
30.3
30.6
32.7

41.3
48.8
29.1
33.3
29.4
41.2
32.4
35.5
32.5
31.7
32.7

40.5
47.2
30.7
31.6
35.2
38.6
36.3
35.4
35.0
34.9
35.3

h = 32A0 (3)A1 (3)
SV-Q
SV-E
GARCH
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

54.6
41.4
56.7
46.2
39.4
42.0
34.2
35.8
27.9
25.6
25.4

43.6
31.7
44.8
37.5
25.9
40.5
29.4
31.5
22.4
22.3
22.1

37.2
39.2
30.4
33.6
26.5
40.4
29.2
33.7
28.3
28.4
28.5

37.6
43.4
28.7
34.1
29.1
40.1
30.1
34.0
29.8
30.6
31.0

39.4
47.1
30.2
34.0
32.2
37.9
32.1
33.8
32.5
32.0
32.7

36.6
42.9
30.0
30.3
35.1
34.6
32.8
31.6
33.3
31.4
32.7

Table 10 (continued): Out-of-sample RMSEs for the volatility prediction in
the fixed-parameter approach
Root mean squared errors (RMSEs) for the out-of-sample prediction of the h-week ahead
yield volatilities (annualized standard deviations) are presented in basis points with h ={4, 8, 16, 32}. For the term structure and GARCH models, the name alone indicates theresults based on equation (31), whereas the name “+ Reg” indicates the results based onequation (32). All parameters are fixed at the in-sample estimates throughout the out-ofsample period from April 16, 2003 to May 27, 2009. The best and the second best resultsin each column of each panel are displayed in bold and italic, respectively.

59

Maturity
h=4
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

0.5

1

2

3

5

10

38.8
37.1
39.8
32.3
31.1
30.6

35.4
33.2
36.4
28.9
27.7
27.1

42.1
34.5
38.1
31.7
30.1
29.9

43.2
34.1
37.7
31.4
30.9
30.7

43.4
33.6
35.6
30.3
29.9
29.7

42.9
34.4
36.6
32.8
32.2
32.2

h=8
A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

36.9
35.9
38.4
31.5
29.8
29.8

32.9
30.6
34.3
26.5
24.9
24.6

39.9
30.9
35.3
28.2
27.5
27.6

41.1
30.5
35.0
28.3
28.8
29.2

41.2
29.4
32.7
27.7
27.2
27.6

40.3
30.2
33.5
31.3
30.0
30.7

h = 16A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

35.0
36.3
37.5
31.6
27.5
28.5

31.2
30.0
33.1
25.7
23.6
23.7

38.3
29.8
34.0
28.1
27.8
28.0

39.6
30.0
33.8
28.4
28.9
30.2

39.9
29.5
31.6
28.6
26.7
28.9

38.4
30.6
31.9
32.3
28.6
31.8

h = 32A1 (3) + Reg
SV-Q + Reg
SV-E + Reg
GARCH + Reg
HAR-RV
MIDAS

33.0
36.0
36.4
29.3
25.9
26.2

29.6
29.7
32.2
23.4
22.2
22.4

36.4
28.7
33.5
25.7
26.3
26.5

37.4
29.2
33.5
25.7
27.5
27.4

37.5
29.0
31.1
26.3
26.0
26.7

34.4
29.3
29.6
30.9
27.0
29.2

Table 11: Out-of-sample RMSEs for the volatility prediction in the varyingparameter approach
RMSEs for the out-of-sample prediction of the h-week ahead yield volatilities (annualized standard deviations) are presented in basis points with h = {4, 8, 16, 32}. For the term structure and
GARCH models, the model parameters are fixed at the in-sample estimates and the regression parameters in equation (32) are re-estimated in a rolling-window fashion. For MIDAS, the weighting
parameters in equation (30) are fixed at the in-sample estimates and (ah,τ , bh,τ ) in equation (29)are re-estimated. For HAR-RV, all parameters in equation (6) are re-estimated. The best and thesecond best results in each column of each panel are displayed in bold and italic, respectively. The
out-of-sample period is from April 16, 2003 to May 27, 2009.

60

Maturity
h=4+ A1 (3)

0.5

1

2

3

5

10

1.036
(0.437)0.867(0.429)0.482(0.325)0.324(0.314)−0.056(0.259)−0.374(0.267)+ SV-Q

0.360(0.142)0.427(0.157)0.439(0.115)0.461(0.114)0.466(0.126)0.361(0.174)+ SV-E

0.555(0.129)0.516(0.204)0.358(0.137)0.420(0.128)0.373(0.120)0.279(0.116)+ GARCH

0.156(0.120)0.553(0.141)

0.355(0.202)0.528(0.208)0.443(0.197)0.832(0.270)1.164(0.527)1.066(0.535)0.650(0.414)0.482(0.377)−0.022(0.313)−0.418(0.320)+ SV-Q

0.362(0.159)0.459(0.198)0.518(0.129)0.542(0.126)0.506(0.146)

0.329(0.213)+ SV-E

0.615(0.170)0.563(0.248)0.432(0.160)0.482(0.142)0.413(0.141)0.314(0.142)+ GARCH

0.273(0.132)0.718(0.193)0.460(0.229)0.578(0.237)0.422(0.213)0.802(0.305)h=8+ A1 (3)Table 12: Estimates (standard errors) of the slope coeﬃcient on the model
forecast in the combined regression

61

Maturity
h = 16+ A1 (3)

0.5

1

2

3

5

10

1.375
(0.631)1.390(0.681)

0.842(0.527)0.665(0.455)0.020(0.375)−0.566(0.342)+ SV-Q

0.302(0.161)0.463(0.236)0.525(0.161)

0.563(0.147)0.470(0.169)0.225(0.288)+ SV-E

0.613(0.229)0.580(0.315)0.488(0.207)0.534(0.169)0.446(0.169)0.354(0.179)+ GARCH

0.304(0.170)0.759(0.290)0.444(0.302)0.626(0.291)0.443(0.246)0.909(0.541)

1.635(0.739)1.632(0.883)0.972(0.742)0.840(0.595)−0.013(0.515)−0.822(0.357)+ SV-Q

0.313(0.167)0.513(0.247)0.542(0.179)0.565(0.166)0.501(0.209)0.235(0.312)+ SV-E

0.542(0.249)0.513(0.350)0.481(0.232)0.516(0.198)0.442(0.201)

0.343(0.207)+ GARCH

0.289(0.196)0.978(0.303)0.485(0.414)0.661(0.456)0.524(0.311)

0.561(0.913)h = 32+ A1 (3)Table 12 (continued): Estimates (standard errors) of the slope coeﬃcient on
the model forecast in the combined regression
The combined forecasting regression is based on the HAR-RV model:��RVt,t+h∆,τ= ah,τ +bh,τ,ih∆
i={4,8,16,32}

�

RVt−i∆,t,τ
+ ch,τ f orecastt,h,τ + uct+h∆,τ ,i∆where f orecastt,h,τ is taken from the term structure and GARCH models. The table
presents the estimates of ch,τ . The standard errors, computed using the method of Newey
and West (1987) with 32 lags, are in parenthesis. The in-sample data from January 4,
1991 to April 9, 2003 are used for the estimation.

62

Maturity
h=4+ A1 (3)+ SV-Q+ SV-E+ GARCH

0.5

1

2

3

5

10

3.7
−2.4
0.9
1.3

5.0
0.6
2.3
6.2

3.3
−3.7
1.7
−0.4

3.2
−7.8
2.5
−3.6

−0.4
−10.5
0.1
−2.0

−1.1
−10.6
−4.7
−3.5

h=8
+ A1 (3)+ SV-Q+ SV-E+ GARCH

6.0
−2.0
3.5
3.2

10.8
1.3
5.7
11.6

7.4
−5.4
4.0
−1.9

6.3
−11.0
4.0
−5.5

−0.2
−14.3
0.3
−2.5

−2.0
−11.5
−6.2
−3.6

h = 16+ A1 (3)+ SV-Q+ SV-E+ GARCH

12.1
2.4
9.9
4.6

21.4
5.0
10.5
11.9

11.4
−5.5
5.7
−1.4

9.3
−10.2
4.5
−4.6

0.2
−11.6
0.1
−1.7

−2.3
−7.1
−7.3
−2.5

25.0
8.0
18.9
−1.5

30.7
7.3
14.1
0.9

13.0
−6.7
6.5
−1.9

10.1
−9.3
4.0
−4.3

−0.1
−8.4
−0.5
−1.1

0.6
−5.1
−6.1
0.2

h = 32+ A1 (3)+ SV-Q+ SV-E+ GARCH

Table 13: % changes in the out-of-sample RMSE by combined forecasts in the
fixed-parameter approach
A combined forecast is based on the HAR-RV model, where a forecast taken from theterm structure and GARCH models is added. The table presents % changes in the RMSE
from the HAR-RV forecast alone to the combined forecast. A negative (positive) number
indicates that the combined forecast decreases (increases) the RMSE. All parameters are
fixed at the in-sample estimates throughout the out-of-sample period from April 16, 2003
to May 27, 2009.

63

Maturity
h=4+ A1 (3)+ SV-Q+ SV-E+ GARCH

0.5

1

2

3

5

10

0.1
−3.5
0.2
1.0

−0.6
−2.7
0.4
1.1

−0.2
−3.6
1.0
0.3

0.3
−5.5
2.9
−1.8

1.1
−6.4
5.2
−1.2

0.0
−9.8
0.8
−1.7

h=8
+ A1 (3)+ SV-Q+ SV-E+ GARCH

−0.2
−3.5
0.8
1.6

−1.4
−3.5
1.1
2.0

−0.2
−5.9
0.8
−1.1

0.7
−8.3
3.8
−4.4

1.7
−11.0
7.3
−2.2

−0.5
−15.7
1.8
−2.3

h = 16+ A1 (3)+ SV-Q+ SV-E+ GARCH

−0.8
−0.7
1.9
2.9

−2.5
−2.6
1.7
3.0

−0.5
−6.1
1.3
−0.2

0.8
−6.2
5.2
−3.9

1.9
−8.8
8.4
−1.1

−0.8
−12.1
3.2
−0.9

h = 32+ A1 (3)+ SV-Q+ SV-E+ GARCH

−2.3
3.4
1.1
0.8

−4.0
−0.8
1.8
3.7

−2.1
−9.0
2.6
−1.0

−1.9
−6.4
6.6
−6.0

0.2
−6.8
6.8
−2.0

−0.8
−7.7
2.8
0.1

Table 14: % changes in the out-of-sample RMSE by combined forecasts in thevarying-parameter approach
A combined forecast is based on the HAR-RV model, where a forecast taken from theterm structure and GARCH models is added. The table presents % changes in the RMSE
from the HAR-RV forecast alone to the combined forecast. A negative (positive) number
indicates that the combined forecast decreases (increases) the RMSE. For the term structure and GARCH models, the model parameters are fixed at the in-sample estimates andthe parameters in the combined HAR-RV regression are re-estimated in a rolling-window
fashion. The out-of-sample period is from April 16, 2003 to May 27, 2009.

64

Maturity
0.5
1
2
3
5
Panel A: RMSEs based on eq. (31) (without the forecasting regression)

10

SV-Q
NC1
NC2
NC3

33.7
33.7
34.5
35.3

35.0
34.4
35.2
35.0

32.8
32.7
33.7
32.3

33.9
33.7
35.0
33.4

28.3
28.4
30.3
28.2

26.9
26.9
29.1
28.0

SV-E
NC1
NC2
NC3

32.0
32.0
31.7
32.2

33.8
33.3
33.0
33.0

32.8
32.8
31.8
31.6

33.6
33.6
32.6
32.4

28.4
28.6
27.7
27.2

27.8
27.6
26.7
26.7

Panel B: RMSEs based on eq. (32) (with the forecasting regression)
SV-Q + Reg
NC1
NC2
NC3

28.7
28.7
28.9
29.7

32.5
32.5
32.6
33.4

31.1
31.1
31.0
31.3

32.4
32.4
32.4
32.4

27.5
27.5
27.6
27.3

26.3
26.3
26.4
26.1

SV-E + Reg
NC1
NC2
NC3

27.8
27.8
27.8
27.9

31.9
31.9
31.8
31.9

31.5
31.5
31.5
31.6

32.3
32.3
32.3
32.3

27.1
27.1
27.1
26.9

25.8
25.8
25.8
25.7

Table 15: In-sample RMSEs for the 4-week ahead volatility prediction withand without the non-arbitrage
Root mean squared errors (RMSEs) for the in-sample prediction of the 4-week ahead yieldvolatilities (annualized standard deviations) are presented in basis points, where panels
A and B present the results based on equations (31) and (32), respectively. The rows
labeled SV-Q and SV-E show the original results with the non-arbitrage imposed. NCi(i = 1, 2, 3) indicates “No Constraint”, beyond the linear regression (34), on the crosssectional relation with the number corresponding to the stage where the non-arbitrage isinvolved: it is involved with both identification of the state vector Xt and estimation of
the covariance matrix Σt in stage 1; with only the latter in stage 2; and with none in stage
3, where Xt is given by the conventional level, slope, and curvature factors. The in-sample
period is from January 4, 1991 to April 9, 2003.

65

y0.5
y1
y2
y3
y5
√
Panel A: Parameter estimates for dyt,τ = (α + βyt,τ )dt + σ yt,τ dWt

y10

α × 102−0.270(0.689)−0.278(1.001)0.253(2.404)0.787(1.723)1.503(1.803)2.626(1.992)β−0.043(0.125)−0.040(0.178)−0.133(0.405)−0.215(0.282)−0.311(0.280)−0.446(0.289)σ × 102

3.248(0.091)3.745(0.105)4.229(0.119)4.212(0.118)4.019(0.113)3.774(0.106)

Panel B: Means (standard deviations) of % diﬀerences
between two conditional standard deviations with and without β
h=4

0.38
(0.07)

0.36(0.07)

0.67(0.04)

0.96(0.09)

1.31(0.11)

1.80(0.13)

h=8

0.76(0.15)

0.71(0.13)

1.34(0.08)

1.92(0.18)

2.63(0.22)

3.62(0.26)

h = 161.53(0.30)

1.43(0.27)

2.70(0.16)

3.87(0.36)

5.29(0.46)

7.29(0.53)

h = 32

3.11(0.62)

2.90(0.57)

5.45(0.33)

7.82(0.74)

10.70(0.96)

14.75(1.12)

Table A: Impact of omitting the mean reversion on the realized volatility
In Panel A, the estimates (standard errors) for the CIR-type square-root process fitted toeach yield are presented. In Panel B, the means (standard deviations) for the time seriesof���vart [yt+h∆,τ ; σ]100−1vart [yt+h∆,τ ; α, β, σ]
are presented, where vart [yt+h∆,τ ; α, β, σ] and vart [yt+h∆,τ ; σ] are presented in equations(38) and (39). The in-sample data from January 4, 1991 to April 9, 2003 are used for theestimation.

66

Panel A: the first-step comparison
y0.5

y1

y2

y3

y5

y10

Minimum

Level
Slope
Curvature

0.00
0.00
0.01

0.00
0.01
0.02

0.00
0.00
−0.01

0.00
0.00
−0.02

0.01
0.00
0.00

0.25
0.16
0.39

MedianLevel
Slope
Curvature

0.00
0.00
0.00

0.00
0.00
0.00

0.00
0.00
0.00

0.01
0.01
0.00

0.06
0.05
0.03

1.20
0.78
0.73

Maximum

Level
Slope
Curvature

0.00
0.00
0.00

0.00
0.00
−0.01

0.00
0.00
0.00

0.01
0.01
0.01

0.06
0.14
0.04

1.29
2.20
0.62

θ

�

y1

y3

y5

Panel B: the second-step comparisonrMinimum

Level
Slope
Curvature

−0.09
−0.06
−0.16

−0.39
−0.25
−0.62

0.80
0.52
1.34

0.03
0.02
0.06

−0.06
−0.04
−0.10

−0.14
−0.10
−0.23

MedianLevel
Slope
Curvature

−0.41
−0.27
−0.25

−1.90
−1.24
−1.16

3.89
2.55
2.37

0.13
0.08
0.08

−0.26
−0.17
−0.16

−0.66
−0.42
−0.41

Maximum

Level
Slope
Curvature

−0.44
−0.75
−0.22

−2.03
−3.46
−0.99

4.16
7.09
2.03

0.14
0.23
0.06

−0.28
−0.47
−0.13

−0.71
−1.18
−0.34

Table C1: Comparison to the closed-form solution using the A0 (3) model

67

Panel C: the third-step comparison
r

θ

�

y1

y3

y5

Minimum

Level
Slope
Curvature

0.04
0.71
2.10

−4.98
−1.28
−2.28

5.24
0.34
−8.07

0.02
−0.16
−0.46

−0.13
0.18
0.52

−0.39
0.22
0.64

MedianLevel
Slope
Curvature

0.50
−0.37
0.03

−0.88
2.78
−1.78

−7.17
−4.07
−1.90

−0.10
0.07
0.01

0.06
−0.04
−0.07

−0.10
−0.01
−0.25

Maximum

Level
Slope
Curvature

1.26
0.02
−0.78

3.13
−1.91
2.98

−13.49
−8.58
−1.53

−0.31
0.03
0.16

0.40
−0.17
−0.13

0.52
−0.63
−0.11

Table C1 (continued): Comparison to the closed-form solution using the A0 (3)model
Approximation errors, defined as the diﬀerence between the approximate and true yields
or state variables, are presented in bps. The errors are evaluated at nine states takenfrom the actual data, where the level, slope, and curvature factors take the minimum,
median, or maximum value. Panel A presents the results for the first-step comparison,
where the true values of both the parameter and state vectors are given as input for theapproximation. Panel B presents the results for the second-step comparison, where thetrue value of the parameter vector alone is given. Panel C presents the results for thethird-step comparison, where no prior information is given.

68

κ1
κ2
κ3
θ̄
σ11 × 104
σ22 × 104
σ33 × 104
σ12 × 104
σ13 × 104
λ02
λ03
λ112
λ122
ς × 102

True
Values
0.947
0.017
0.957
0.146
0.722
1.220
8.036
0.095
−0.9850.030−0.013−0.195−0.514
0.061

Approximation
Est.
S.E.
0.966
(0.088)0.019(0.002)0.964(0.089)0.135(0.006)0.730(0.030)1.257(0.066)8.240(0.441)
0.099(0.036)−1.019(0.115)0.030(0.019)−0.014(0.008)−0.191(0.037)−0.518(0.285)0.061(0.001)

Diﬀ.

%Diﬀ.

−0.019−0.002−0.0060.011−0.008−0.037−0.204−0.004
0.034
0.000
0.001
−0.004
0.004
0.000

−2.0
−13.6
−0.7
7.5
−1.1
−3.0
−2.5
−3.8
−3.4
−0.9
−4.9
2.0
−0.8
0.4

Table C2: Parameter estimates of the A0 (3) model by the approximation
methodThe true parameter values are taken from the estimates presented in Table 4. Est. and
S.E. are the estimates and standard errors, respectively, obtained by the approximation
method. Diﬀ. and % Diﬀ. are diﬀerences and % diﬀerences between the true values andestimates, respectively.

69

Panel A: Based on the SV-Q model
y0.5

y1

y2

y3

y5

y10

Minimum

Level
Slope
Curvature

0.11
0.17
0.37

0.04
0.22
0.13

0.12
0.35
−0.15

0.13
0.44
−0.11

−0.03
0.55
−0.04

0.19
0.79
0.99

MedianLevel
Slope
Curvature

−0.03
−0.18
−0.05

−0.02
−0.13
−0.05

0.00
−0.08
−0.08

0.03
−0.02
−0.08

0.10
0.01
−0.16

1.41
0.45
0.14

Maximum

Level
Slope
Curvature

0.06
−0.14
−0.23

0.14
−0.14
−0.21

0.23
−0.19
−0.14

0.27
−0.19
−0.04

0.13
−0.08
−0.05

1.34
2.64
0.83

Panel B: Based on the SV-E model
y0.5

y1

y2

y3

y5

y10

Minimum

Level
Slope
Curvature

0.12
0.19
0.39

0.11
0.16
0.35

0.09
0.12
0.20

0.10
0.07
0.10

0.07
0.04
0.08

0.46
0.42
1.39

MedianLevel
Slope
Curvature

−0.03
−0.20
−0.04

−0.05
−0.16
−0.04

−0.08
−0.09
−0.04

−0.04
−0.06
−0.08

−0.03
−0.04
−0.09

1.38
0.30
0.16

Maximum

Level
Slope
Curvature

0.04
−0.18
−0.25

0.02
−0.18
−0.20

0.00
−0.16
−0.10

−0.06
−0.12
−0.03

−0.01
0.11
0.02

0.85
2.37
0.41

Table C3: Comparison to the MC solution using the SV-Q and SV-E models
Approximation errors, defined as the diﬀerence between the approximate and MC yields,
are presented in bps. The errors are evaluated at nine states taken from the actual data,
where the level, slope, and curvature factors take the minimum, median, or maximum
value. Panels A and B present the results based on the SV-Q and SV-E models, respectively.

70

Maturity
0.5
Panel A: Based on the SV-Q model

1

2

3

5

10

Minimum

Level
Slope
Curvature

−0.16
−0.24
−0.46

−0.16
−0.31
−0.61

−0.17
−0.30
−0.67

−0.15
−0.28
−0.63

−0.13
−0.23
−0.64

−0.12
−0.22
−0.76

MedianLevel
Slope
Curvature

0.00
−0.14
−0.06

−0.06
−0.09
−0.03

−0.06
−0.08
−0.05

−0.06
−0.06
−0.02

−0.05
0.01
−0.01

−0.04
0.10
−0.02

Maximum

Level
Slope
Curvature

0.44
−0.06
−0.01

0.49
−0.12
−0.07

0.51
−0.09
−0.09

0.45
−0.08
−0.05

0.33
−0.06
0.03

0.18
−0.01
0.22

Panel B: Based on the SV-E model

Minimum

Level
Slope
Curvature

0.06
0.05
0.07

0.04
0.03
0.05

0.04
0.01
0.01

0.02
0.01
0.05

0.01
−0.01
0.06

−0.07
−0.07
−0.38

MedianLevel
Slope
Curvature

0.02
−0.09
0.02

0.02
−0.07
0.02

0.02
−0.04
0.01

0.02
−0.02
0.02

0.01
−0.01
0.01

−0.07
0.01
−0.05

Maximum

Level
Slope
Curvature

0.00
−0.04
−0.06

−0.05
−0.01
−0.04

−0.07
−0.02
−0.01

−0.08
−0.03
0.00

−0.08
−0.01
0.02

−0.22
0.01
0.03

Table C4: Accuracy to Et [Ỹ (Xt+h∆ , τ )] for h = 32
Approximation errors, defined as the diﬀerence between the approximate and MC moments, are presented in bps. The errors are evaluated at nine states taken from the actual
data, where the level, slope, and curvature factors take the minimum, median, or maximum value. Panels A and B present the results based on the SV-Q and SV-E models,
respectively.

71

Maturity
0.5
Panel A: Based on the SV-Q model

1

2

3

5

10

Minimum

Level
Slope
Curvature

0.63
0.87
−0.21

0.72
1.36
0.11

0.81
1.53
0.53

0.74
1.41
0.68

0.59
1.00
0.50

0.41
0.49
0.33

MedianLevel
Slope
Curvature

0.56
−0.22
0.64

0.74
0.09
0.73

0.79
0.31
0.58

0.75
0.39
0.33

0.64
0.50
−0.02

0.52
0.46
−0.31

Maximum

Level
Slope
Curvature

−1.52
−0.77
0.05

−1.90
−1.01
−0.34

−1.89
−1.20
−0.79

−1.53
−1.26
−0.90

−0.88
−1.21
−0.85

−0.20
−1.19
−0.82

Panel B: Based on the SV-E model

Minimum

Level
Slope
Curvature

−0.51
0.37
0.98

−0.62
0.33
0.53

−0.62
0.29
0.10

−0.55
0.27
−0.08

−0.43
0.24
−0.23

−0.28
0.23
−1.28

MedianLevel
Slope
Curvature

0.06
−0.23
−0.12

−0.06
0.02
−0.15

−0.16
0.13
−0.12

−0.16
0.10
−0.03

−0.12
0.01
0.12

−0.15
−0.09
0.22

Maximum

Level
Slope
Curvature

1.99
0.08
−0.02

1.76
0.01
−0.17

1.49
−0.12
−0.36

1.28
−0.32
−0.44

1.07
−0.66
−0.52

0.19
−1.37
−0.57

�

Table C5: Accuracy to vart [Ỹ (Xt+h∆ , τ )]/h∆ for h = 32
Approximation errors, defined as the diﬀerence between the approximate and MC moments, are presented in bps. The errors are evaluated at nine states taken from the actual
data, where the level, slope, and curvature factors take the minimum, median, or maximum value. Panels A and B present the results based on the SV-Q and SV-E models,
respectively.

72

0.10

0.030

0.08

0.025

0.020
0.06
0.015
0.04
0.010
0.02
0.005
0.00

0.000
91

93

95

97

99

01

03

05

07

09

0.10

91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

0.030

0.025

0.08

0.020
0.06
0.015
0.04
0.010
0.02

0.005

0.00

0.000
91

93

95

97

99

01

03

05

07

09

0.030

0.10

0.025

0.08

0.020
0.06
0.015
0.04
0.010
0.02

0.005

0.00

0.000
91

93

95

97

99

01

03

05

07

09

0.10
r
0.08

y_0.5

0.06

0.04

0.02

0.00
91

93

95

97

99

01

03

05

07

09

0.10
theta
0.08

y_10

0.06

0.04

0.02

0.00
91

93

95

97

99

01

03

05

07

09

0.06
epsilon
0.04
0.02
0.00
-0.02
-0.04
-0.06
-0.08
91

93

95

97

99

01

03

05

07

09

0.10

0.10

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02

0.00

0.00
91

93

95

97

99

01

03

05

07

09

0.10

0.10

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02

0.00

91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

0.00
91

93

95

97

99

01

03

05

07

09

0.06

0.06

0.04

0.04

0.02

0.02

0.00

0.00

-0.02

-0.02

-0.04

-0.04

-0.06

-0.06

-0.08

-0.08
91

93

95

97

99

01

03

05

07

09

0.030

0.025

0.020

0.015

0.010

0.005

0.000
91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

0.030

0.025

0.020

0.015

0.010

0.005

0.000

0.030

0.025

0.020

0.015

0.010

0.005

0.000

0.030

0.025

0.020

0.015

0.010

0.005

0.000
91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

0.030

0.025

0.020

0.015

0.010

0.005

0.000

0.030

0.025

0.020

0.015

0.010

0.005

0.000

0.030

0.025

0.020

0.015

0.010

0.005

0.000
91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

0.030

0.025

0.020

0.015

0.010

0.005

0.000

0.030

0.025

0.020

0.015

0.010

0.005

0.000

0.030

0.025

0.020

0.015

0.010

0.005

0.000
91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

91

93

95

97

99

01

03

05

07

09

0.030

0.025

0.020

0.015

0.010

0.005

0.000

0.030

0.025

0.020

0.015

0.010

0.005

0.000