Computing Arbitrage-Free Yields in Multi-Factor Gaussian Shadow-Rate Term Structure Models

Finance and Economics Discussion Series Divisions of Research & Statistics and Monetary Affairs Federal Reserve Board, Washington, D.C. Computing Arbitrage-Free Yields in Multi-Factor Gaussian Shadow-Rate Term Structure Models Marcel A. Priebsch 2013-63 NOTE: Staff working papers in the Finance and Economics Discussion Series (FEDS) are preliminary materials circulated to stimulate discussion and critical comment. The analysis and conclusions set forth are those of the authors and do not indicate concurrence by other members of the research staff or the Board of Governors. References in publications to the Finance and Economics Discussion Series (other than acknowledgement) should be cleared with the author(s) to protect the tentative character of these papers.

Computing Arbitrage-Free Yields in Multi-Factor Gaussian Shadow-Rate Term Structure Models Preliminary Draft Marcel Priebsch∗ First Draft: March 11, 2013 This Draft: September 24, 2013 Abstract This paper develops a method to approximate arbitrage-free bond yields withinatermstructuremodelinwhichtheshortratefollowsaGaussianprocess censored at zero (a (cid:16)shadow-rate model(cid:17) as proposed by Black, 1995). The censoring ensures that model-implied yields are constrained to be positive, but it also introduces non-linearity that renders standard bond pricing formulas inapplicable. Inparticular,yieldsarenotlinearfunctionsoftheunderlyingstate vector as they are in a(cid:30)ne term structure models (see Piazzesi, 2010). Existing approaches towards computing yields in shadow-rate models su(cid:27)er from high computational burden or low accuracy. In contrast, I show that the technique proposed in this paper is su(cid:30)ciently fast for single-step estimation of a threefactor shadow-rate term structure model, and su(cid:30)ciently accurate to evaluate yields to within approximately half a basis point. ∗Federal Reserve Board, Washington D.C., marcel.a.priebsch@frb.gov. The analysis and conclusions set forth in this paper are those of the author and do not indicate concurrence by other members of the research sta(cid:27) or the Board of Governors of the Federal Reserve System. 1

1 Introduction In late-2008, short-term nominal interest rates in the U.S. fell to their e(cid:27)ective (cid:16)zero lower bound(cid:17) (see Bernanke et al., 2004). Since standard Gaussian term structure models do not rule out the possibility of negative model-implied yields, they provide a poor approximation to the behavior of nominal yields when the lower bound is binding (Kim and Singleton, 2012; Christensen and Rudebusch, 2013; Bauer and Rudebusch, 2013). Kim and Singleton (2012) (cid:28)nd that shadow-rate models in the spirit of Black (1995) successfully capture yield-curve properties observed near the zero lower bound. However, arbitrage-free multi-factor versions of these models tend to be computationally intractable (Christensen and Rudebusch, 2013). Gorovoi and Linetsky (2004) show that bond prices in a one-factor shadow-rate model can be computed analytically by an eigenfunction expansion, but their approach does not generalize to multiple dimensions. Kim and Singleton (2012) and Ichiue and Ueno (2007) successfully estimate shadow-rate models with up to two factors, but they compute bond prices using discretization schemes that are subject to the curse of dimensionality. Christensen and Rudebusch (2013) use a yield formula proposed by Krippner (2012) to estimate shadow-rate Nelson-Siegel models with up to three factors, but Krippner’s derivation deviates from the usual no-arbitrage approach. Bauer and Rudebusch (2013) evaluate bond prices by Monte Carlo simulation for given model parameters from an unconstrained Gaussian term structure model, but they do not estimate a shadow-rate version of the model due to the computational burden. This paper develops and applies a new technique for fast and accurate approximation of arbitrage-free zero-coupon bond yields in multi-factor Gaussian shadow-rate models of the term structure of interest rates. The computational complexity of the 2

method does not increase with the number of yield curve factors, and, empirically, it produces yields that are accurate to within about half a basis point. The method is su(cid:30)cientlyfasttoestimatea(cid:29)exible,arbitrage-free,three-factortermstructuremodel in which the shadow rate follows a Gaussian process. For illustration purposes, I estimate such a model by quasi-maximum likelihood on a sample of U.S. Treasury yields, using the unscented Kalman (cid:28)lter to account for the non-linear mapping between factors and yields. 2 Model Consider(cid:28)rstthestandard, continuous-timeN-factorGaussiantermstructuremodel. In particular, let WP be N-dimensional standard Brownian motion on a complete t probability space (Ω,F,P) with canonical (cid:28)ltration {F } . Assume there is a prict t≥0 ing measure Q on (Ω,F) that is equivalent to P , and denote by W Q Brownian motion t Q under as derived from Girsanov’s Theorem (Karatzas and Shreve, 1991). Suppose N latent factors (or states) representing uncertainty underlying term-structure securities follow the multivariate Ornstein-Uhlenbeck process dX = (Kµ +KµX )dt+ΣdWµ (1) t 0 1 t t were µ ∈ {P,Q}. Let the short rate be r = ρ +ρ ·X . (2) t 0 1 t Then by de(cid:28)nition, the arbitrage-free time t price of a zero-coupon bond maturing at time T is given by ˆ (cid:20) (cid:18) T (cid:19)(cid:21) PT = E Q exp − r ds (3) t t s t 3

with associated zero-coupon bond yield logPT yT = − t . (4) t T −t Bond prices (and hence yields) can equivalently be de(cid:28)ned in terms of forward rates: ˆ (cid:18) T (cid:19) d PT = exp − fsds ⇔ fT = − logPT (5) t t t dT t t where fT denotes the instantaneous time T forward rate e(cid:27)ective at time t. t SinceX isaGaussianprocess(KaratzasandShreve,1991),itfollowsfrom(2)that t the short rate r takes on strictly negative values with strictly positive probability. To t modify the model in a way that accounts for the zero lower bound on nominal yields, Black (1995) proposes to think of r as a shadow short rate (and, analogously, of PT, t t yT, and fT as shadow bond price, shadow yield, and shadow instantaneous forward t t rate, respectively) and de(cid:28)ne the observed short rate as the shadow rate censored at zero: r = max{r ,0}. (6) t t With the observed short rate r in place of the shadow rate r , the observed bond t t price PT, yield yT, and instantaneous forward rate fT are then de(cid:28)ned as in (3)(cid:21)(5). t t t 3 Parameterizing the Lower Bound The theoretical argument for a lower bound at zero on the nominal short rate (and hence on nominal yields) is based on arbitrage between bonds and currency (Black, 1995). In practice, the two assets may not be perfect substitutes for reasons such as convenience, default risk, or legal requirements. This may push the empirical lower 4

bound into slightly negative or slightly positive territory. The derivations in Section 2 are easily modi(cid:28)ed to accommodate a lower bound at r (cid:54)= 0. In particular, suppose min r = max{r ,r } = r +max{r −r ,0}. t t min min t min Then, ˆ 1 (cid:20) (cid:18) T (cid:19)(cid:21) yT = − E Q exp − r ds t T −t t s t ˆ 1 (cid:20) (cid:18) T (cid:19)(cid:21) Q = r − E exp − max{r −r ,0}ds . min T −t t s min t The last term is equal to the expression for the yield when the lower bound is zero, except that r is subtracted from the shadow rate. Therefore, since r = ρ +ρ ·X , min s 0 1 s when the lower bound is nonzero we can compute zero-coupon yields as if the bound were zero, with ρ − r in place of ρ , and then add r to the (cid:28)nal result. The 0 min 0 min lower bound r can be set to a speci(cid:28)c value based on a priori reasoning, or treated min as a free parameter in estimation. 4 Bond Price Computation A central task in term-structure modeling is the (analytical and/or numerical) computation of arbitrage-free bond prices (and hence yields) based on equation (3). Section 4.1 reviews the standard approach using di(cid:27)erential equations formalized by Du(cid:30)e and Kan (1996), best suited to a(cid:30)ne models. While it can be adapted to the shadow-rate framework, it loses much of its analytical tractability, and becomes computationally infeasible as the number of factors increases. Section 4.2 discusses an alternative method proposed by Krippner (2012). It de(cid:28)nes a pseudo(cid:21)forward rate 5

that satis(cid:28)es the lower bound (though di(cid:27)ers from the arbitrage-free forward rate), and uses relationship (5) to approximate bond prices. Finally, Section 4.3 proposes a new approximation technique for yields in the shadow-rate model based on the expansion of a cumulant-generating function. 4.1 Partial Di(cid:27)erential Equation ´ Like any conditional expectation of an F T -measurable random variable, e− 0 trsdsP t T (the time t price of a zero-coupon bond maturing at time T in the unconstrained Gaussian model, discounted to time 0) follows a martingale under Q . This is an immediate consequence of the de(cid:28)nition of a martingale, after an application of the Law of Iterated Expectations. Using It(cid:9)o’s Lemma and the Martingale Representation Theorem, we can therefore represent PT by the function D(X ,T−t), where D solves t t the partial di(cid:27)erential equation (PDE) 1 D (x,τ) = D(cid:62)(x,τ)(K Q +K Q x)+ tr (cid:0) Σ(cid:62)D (x,τ)Σ (cid:1) −(ρ +ρ ·x)D(x,τ) (7) τ x 0 1 2 xx 0 1 with boundary condition D(x,0) = 1. Using a separation-of-variables argument, it can be veri(cid:28)ed that D(x,τ) = eA(τ)+B(τ)·x (8) solves (7), where A and B in turn solve ordinary di(cid:27)erential equations (ODEs) in terms of the model parameters, 1 A(cid:48)(τ) = K Q ·B(τ)+ B(τ)(cid:62)ΣΣ(cid:62)B(τ)−ρ (9) 0 2 0 B(cid:48)(τ) = (K Q )(cid:62)B(τ)−ρ (10) 1 1 6

with A(0) = 0, B(0) = 0. Reducing problem (7) to the system of ODEs (9)(cid:21)(10) simpli(cid:28)es the numerical computation of bond prices substantially as it reduces the dimensionality of the problem from N +1 to 1 (the time dimension). Direct computation reveals that B(τ) = ((K 1 Q )−1)(cid:62)(I N −e(K 1 Q )(cid:62)τ)ρ 1 , (11) Q assuming K is invertible. 1 A PDE analogous to (7) can be set up for the observed bond price in the shadowrate model, PT. The only required modi(cid:28)cation is to replace the expression for the t shadowshortrate, r = ρ +ρ ·x, bythatfortheobservedshortrate, r = max{ρ +ρ · 0 1 0 1 x,0}. Unfortunately,whenthisnon-linearityisintroduced,theseparation-of-variables procedure no longer applies, and no solution as straightforward as (8) is available. It is possible to solve the modi(cid:28)ed version of (7) directly by numerical methods. This is the approach taken by Kim and Singleton (2012). It requires discretizing τ and x on a multidimensional grid, which is computationally intensive and subject to the curse of dimensionality. Kim and Singleton (2012) therefore do not estimate models with more than two factors. 4.2 Forward Rate Approximation Krippner (2012) proposes an alternative approach to computing yields in shadow-rate models, which is implemented empirically by Christensen and Rudebusch (2013). It is based on an approximation to the forward rate fT. Substituting for PT from the t t shadow-rate version of (3), and di(cid:27)erentiating, we obtain (cid:34) ´ (cid:35) fT = E Q e− t Trsds r = E Q T [r ] = E Q T [max{r ,0}] (12) t t PT T t T t T t 7

Q where , de(cid:28)ned by the Radon-Nikodym derivative T ´ dQ e− 0 Trsds T = ´ , dQ EQ (cid:2) e− 0 Trsds (cid:3) is referred to as the (cid:16)T-forward measure.(cid:17) Equation (12) says that today’s time T forward rate is equal to today’s expectation under the T-forward measure of the time T shortrate. Itcanbeveri(cid:28)eddirectlyfrom(12)andtheLawofIteratedExpectations that fT is a martingale under Q (note that r ≡ fT by de(cid:28)nition). Note also that t T T T (12) implies that the forward rate is subject to the same lower bound (6) as the short rate, by monotonicity of the mathematical expectation. Analogously, (cid:34) ´ (cid:35) fT = E Q e− t Trsds r = E Q T [r ] (13) t t PT T t T t expresses the shadow forward rate as the expectation under the shadow T-forward measure of the future shadow short rate. Again, fT is a martingale under Q . Unlike t T the observed forward rate, it is not, however, constrained to be non-negative. The distribution of the shadow forward rate fT under Q can be derived more t T explicitly. First, from (5) and (8), fT = −A(cid:48)(T −t)−B(cid:48)(T −t)·X . t t Therefore, by Ito(cid:9)’s Lemma and the Martingale Representation Theorem, ˆ T fT = fT − B(cid:48)(T −s)(cid:62)ΣdW Q T. (14) T t s t 8

Thus, fT is Gaussian under Q conditional on F , with T T t E Q T[fT] =fT t T t ˆ (cid:34) (cid:35) (cid:18) T (cid:19)2 Var Q T[fT] =E Q T B(cid:48)(T −s)(cid:62)ΣdW Q T ] t T t s t ˆ (cid:18) T (cid:19) =ρ(cid:62) 1 e(K 1 Q )(cid:62)(T−s)ΣΣ(cid:62)eK 1 Q (T−s)ds ρ 1 . (15) t (cid:124) (cid:123)(cid:122) (cid:125) ω(T−t) The (cid:28)nal equality uses the Ito(cid:9) Isommetry and (11).1 Krippner (2012) takes advantage of this distributional property of fT ≡ r . He T T de(cid:28)nes a pseudo(cid:21)forward rate as a hybrid between observed forward rate (12) and shadow forward rate (13): (cid:34) ´ (cid:35) fT = E Q e− t Trsds r = E Q T [r ] = E Q T [max{r ,0}]. (16) t t PT T t T t T ˙ t This is the expectation under the shadow T-forward measure of the observed time T short rate (while the shadow-model-implied forward rate consistent with the absence of arbitrage is given by (12) as the expectation under the observed T-forward measure of the observed time T short rate).2 This rate is, by monotonicity, subject to lower bound (6). It is, moreover, relatively straightforward to compute: Lemma A.1 in 1The integral in (15) has the same form as (A.2) in Appendix A and therefore can be computed analytically as in (A.3). 2Krippner (2012) motivates his derivation in terms of options on shadow bonds. To derive (16) from his equations (12) and (13), replace the call option price by (cid:104) ´ (cid:110) (cid:104) ´ (cid:105) (cid:111)(cid:105) E Q t e− t Trsdsmax E Q T e− T T+δrsds −1,0 . Then interchange the limit operations and expectation, and evaluate. 9

Appendix A implies that (cid:18) fT (cid:19) (cid:18) fT (cid:19) fT = E Q T (cid:2) max{fT,0} (cid:3) = fTΦ t +ω(T −t)φ t . (17) t t T t ω(T −t) ω(T −t) ˙ This is formula (32) in Krippner (2012).3 Note from (17) that fT/fT → 1 as fT t t t ˙ increases or ω(T − t) decreases. That is, as the lower bound becomes less binding, the wedge between fT and the shadow forward rate fT (which is the arbitrage-free t t ˙ forward rate in a Gaussian model without lower bound) shrinks. Zero-coupon bond prices PT and yields yT can be approximated by substituting t t ˙ ˙ fT from (17) into (5) and (4). t ˙ 4.3 Cumulants Since the PDE approach to pricing bonds in shadow-rate models becomes computationally intractable as the number of factors increases, and the approach proposed by Krippner (2012) relies on a forward rate that is not equal to the arbitrage-free forward rate, I propose a new cumulant-based technique to approximating yields in Gaussian shadow-rate models. ´ (cid:104) (cid:16) (cid:17)(cid:105) The quantity logPT = logE Q exp − T r ds appearing in the shadow-rate t t t s version of (4) is the conditional cumulant-generating function4 under Q , evaluated at 3Krippner (2012) uses the parametrization proposed by Chen (1995) and hence obtains his formula (31) for ω as a special case of (15). Similarly, the results derived by Christensen and Rudebusch(2013)are(essentially)specialcasesof (15)and(17)undertheirNelson-Siegelparametrization (where some of the derivations must be modi(cid:28)ed appropriately to account for the fact that their Q matrix K is singular). Don Kim (personal communication) independently derives the general 1 expression for ω in (15) by generalizing Krippner’s (2012) computations directly. 4The cumulant-generating function of a random variable X is de(cid:28)ned as the logarithm of its moment-generating function (for example, see Severini, 2005). 10

´ −1, of the random variable RT ≡ T r ds. It has the series representation t t s logE Q(cid:2) exp (cid:0) −RT (cid:1)(cid:3) = (cid:88) ∞ (−1)j κ Q j (18) t t j! j=1 where κ Q is the jth cumulant of RT under Q . An approximation to the zero-coupon j t yield yT can therefore be computed based on a (cid:28)nite number of terms in the series in t (18). Below, I consider the (cid:28)rst-order approximation ˆ 1 1 (cid:20) T (cid:21) y˜T = κ Q = E Q r ds (19) t τ 1 T −t t s t and the second-order approximation ˆ ˆ 1 (cid:18) 1 (cid:19) 1 (cid:18) (cid:20) T (cid:21) 1 (cid:20) T (cid:21)(cid:19) y˜˜T = κ Q − κ Q = E Q r ds − Var Q r ds (20) t T −t 1 2 2 T −t t s 2 t s t t where I make use of the fact that the (cid:28)rst two cumulants of any random variable coincide with its (cid:28)rst two centered moments.5 The (cid:28)rst-order approximation (19) is equivalent to the method proposed independently and contemporaneously by Ichiue and Ueno (2013). I present it here both for comparison and to assess its relative performance in Section 5 below. I will, however, mostly focus on the second-order approximation (20), which I argue is particularly promising a priori because it is exact in the Gaussian benchmark case.6 It can therefore be expected to perform well both for short maturities (where the higher-order terms in (18) are relatively small), and for long maturities as long as Q [r < 0] is t T small for large T (in which case r = max{r ,0} will behave approximately like a t t 5Higher-orderapproximationsfollowingthesamegenerallogicarepossible,buttheyareincreasingly computationally costly while generating decreasing marginal bene(cid:28)ts in terms of precision. 6The third- and higher-order cumulants of a Gaussian random variable are zero, so that (20) coincides with the usual a(cid:30)ne-Gaussian yield formula in that case. 11

Gaussian process over su(cid:30)ciently long horizons). Indeed, empirically, the secondorder approximation turns out to be highly accurate across maturities both during normal times and when rates are low (see Section 5). 4.3.1 Computation of the First Moment Evaluating the (cid:28)rst- and second-order approximations (19) and (20) to zero-coupon yields requires computation of the (cid:28)rst two cumulants (equivalently, centered mo- ´ ments) of RT = T r ds. This subsection will be concerned with the (cid:28)rst moment. t t s As an initial step, ˆ ˆ (cid:20) T (cid:21) T Q Q E r ds = E [r ] ds (21) t s t s t t byanapplicationofFubini’sTheorem. Sincer = max{r ,0}andr ∼ N(µ ,σ2 ) s s s t→s t→s with known expressions for µ and σ2 in terms of the model parameters (as shown t→s t→s in Appendix A), it follows from Lemma A.1 in Appendix A that ˆ ˆ (cid:20) T (cid:21) T (cid:20) (cid:18) µ (cid:19) (cid:18) µ (cid:19)(cid:21) Q t→s t→s E r ds = µ Φ +σ φ ds (22) t s t→s σ t→s σ t t t→s t→s where φ and Φ denote the standard normal probability density function (pdf) and Q cumulative distribution function (cdf), respectively. That is, we can compute E [r ] t s analytically up to the standard normal cdf, which software such as Matlab is able to evaluate precisely and e(cid:30)ciently. The (cid:28)rst cumulant of RT can then be computed by t Q numerical integration of E [r ] over the time dimension, as in (22). t s 12

4.3.2 Computation of the Second Moment ´ Once we know the (cid:28)rst moment of RT = T r ds, it remains to evaluate t t s ˆ ˆ ˆ (cid:34) (cid:35) (cid:18) T (cid:19)2 T s Q Q E r ds = 2 E [r r ] duds (23) t s t u s t t t where the equality uses Fubini’s Theorem and symmetry of the integrand. Since r = max{r ,0} and r = max{r ,0}, and (r ,r ) are jointly normally distributed u u s s u s with mean (µ ,µ ), variances (σ2 ,σ2 ), and covariance σ (see Appendix t→u t→s t→u t→s t→u×s A), we obtain from Lemma A.2 in Appendix A: ˆ (cid:34) (cid:35) (cid:18) T (cid:19)2 Q E r ds t s t ˆ ˆ (cid:40) T s =2 (µ µ +σ )Φd(−ς ,−ς ;χ ) t→u t→s t→u×s 2 t→u t→s t→u×s t t (cid:32) (cid:33) ς −χ ς t→u t→u×s t→s +σ µ φ(ς )Φ t→s t→u t→s (cid:112) 1−χ2 t→u×s (cid:32) (cid:33) ς −χ ς t→s t→u×s t→u +σ µ φ(ς )Φ t→u t→s t→u (cid:112) 1−χ2 t→u×s (cid:114) (cid:32)(cid:115) (cid:33)(cid:41) 1−χ2 ς2 −2χ ς ς +ς2 +σ σ t→u×sφ t→u t→u×s t→u t→s t→s duds (24) t→u t→s 2π 1−χ2 t→u×s where ς = µt→j, j ∈ {u,s}, χ = σt→u×s , and Φd denotes the decumulative t→j σt→j t→u×s σt→uσt→s 2 bivariate normal cdf. Q That is, we can compute E [r r ] analytically up to the bivariate normal cdf, and t u s we can then integrate this expression numerically over u and s to obtain the second cumulant of RT. t 13

4.3.3 Numerical Implementation The following steps summarize the approximation procedure for zero-coupon yields Q Q for a given set of parameters (K ,K ,ρ ,ρ ,Σ) and state vector X : 0 1 0 1 t 1. Compute the conditional mean and covariance matrix of (r ,r ), for u,s ≥ t, u s using (A.4) and (A.5). Q Q 2. Using the results from step 1, compute E [r ] and E [r r ], for u,s ≥ t, as t s t u s described in 4.3.1 and 4.3.2. 3. Integrate E Q [r ] numerically to obtain E Q(cid:2) RT (cid:3) , and integrate E Q [r r ] numert s t t t u s ically to obtain E Q [ (cid:0) RT (cid:1)2 ]. t t 4. Using the moments computed in step 3, approximate yT by y˜T or y˜˜T as de(cid:28)ned t t t in (19) and (20). In terms of numerical implementation, step 1 is straightforward. Step 2 requires evaluation of the univariate and bivariate normal cdf’s. A high-precision, e(cid:30)cient approximation to the univariate normal cdf is built into most computational software packages, so numerically evaluating the (cid:28)rst integrand does not pose a challenge. For the bivariate normal cdf, I implement a vectorized version of an algorithm proposed by Genz (2004) which achieves double machine precision.7 Step 3 is straightforward in principle, though a favorable tradeo(cid:27) between precision and computational burden requires careful choice of quadrature rule and grid. I use composite Gauss-Legendre and Gauss-Lobatto rules with 2(cid:21)20 points per maturity (and corresponding product rules for the double integral in (24)), selected to evaluate y˜T or y˜˜T to an approximate t t minimum numerical precision of 1/ 100 th of a basis point. With fully vectorized Matlab 7Matlab’sbuilt-infunctionmvncdfusesadaptivenumericalintegrationtocomputethebivariate normal cdf. This is slower by several orders of magnitude than the algorithm I use. 14

code,8 I am able to evaluate a full representative sample of model-implied zero coupon yields (approximately 20 years of monthly data across eight maturities) within a fraction of a second. Note that the complexity of the algorithm does not depend on the number of yield curve factors N, so it is not subject to the curse of dimensionality in the same way that some other methods are.9 For illustration purposes, in Appendix B I apply the second-order approximation method to estimate a three-factor shadow-rate model of the U.S. Treasury term structure. 5 Accuracy of Yield Approximation Methods Section 4.3 argued intuitively that the second-order yield approximation (20) should be relatively precise. This section quanti(cid:28)es that claim. To get an initial sense of the relative numerical accuracy, I consider the stylized model used for illustration by Gorovoi and Linetsky (2004), and replicated for the same purpose in Krippner Q Q (2012). It is a one-factor model with ρ = 0.01, ρ = 1, K = 0, K = −0.1, 0 1 0 1 and Σ = 0.02. Gorovoi and Linetsky (2004) derive model-implied yield curves for states X ∈ {−0.06,−0.02,−0.01,0} corresponding to shadow short rates r = t t {−0.05,−0.01,0,0.01}. The four panels of Figure 1 plot the model-implied yield curves for each initial state. Within each panel, I compare four di(cid:27)erent yield approximation schemes: Solving PDE (7) numerically (which, in a one-factor setting, 8That is, without for-loops that grow with the number of quadrature points, states, dates, or maturities in the sample. 9The general approximation methodology I propose has its own curse of dimensionality in that the second-order approximation is substantially more computationally involved than the (cid:28)rst-order approximation (and any higher-order approximation will be substantially more involved than the second-order approximation). In practice, the second-order approximation appears to strike an acceptable balance between precision and computational complexity for many cases of interest, see Section 5 below. 15

is computationally feasible and can be considered the (cid:16)exact(cid:17) solution for comparison purposes), Krippner’s (2012) approach described in Section 4.2, and the (cid:28)rstand second-order approximations proposed in Section 4.3. As the (cid:28)gure shows, the second-order approximation matches the exact PDE solution most closely, and consistently across states. The yield approximation error is uniformly less than one basis point. The (cid:28)rst-order approximation generally overstates yields (an implication of the alternating nature of series expansion (18)), with approximation errors increasing both in yield maturity and the level of the shadow short rate (in both cases, the (cid:28)rstorder approximation is o(cid:27) by an increasingly large convexity adjustment arising from Jensen’s inequality). Krippner’s (2012) method generally undershoots yields, and is relatively more accurate when the shadow short rate is higher. Why this is the case Q can be seen intuitively by comparing the -measure expressions for the arbitrage-free forward rate (12) and Krippner’s (2012) approximate forward rate (16): While both use the same time T short rate r , Krippner’s (2012) formula discounts by the shadow T rate rather than the observed short rate. This means the discount factor tends to be larger than it should be when r is low, reducing the covariance between discount T factor and r , and thus lowering the expectation of their product, fT. T t ˙ To compare the relative performance of the di(cid:27)erent yield approximations in a more realistic empirical setting, I use the estimated model parameters and smoothed states from Appendix B.3 to compute model-implied yield curves for all dates in the sample.10 Since this is a three-factor model, solving PDE (7) numerically is no longer practicable. I therefore replace this benchmark by a simulated yield yˆT that t consistently estimates the true yield yT based on n = 1,000,000 randomly drawn t 10The sample consists of end-of-month observations from January 1990 through December 2012. 16

0.5% 0.4% 0.3% 0.2% 0.1% 0.0% 0 2 4 6 8 10 dleiY 1.2% PDE Krippner 1.0% First‐order Second‐order 0.8% 0.6% 0.4% 0.2% 0.0% 0 2 4 6 8 10 Maturity (years) (a) Yield curves for r =−5%. t dleiY PDE Krippner First‐order Second‐order Maturity (years) (b) Yield curves for r =−1%. t 1.6% 1.4% 1.2% 1.0% 0.8% 0.6% 0.4% 0.2% 0.0% 0 2 4 6 8 10 dleiY 2.0% 1.8% 1.6% 1.4% PDE 1.2% Krippner First‐order 1.0% Second‐order 0.8% 0 2 4 6 8 10 Maturity (years) (c) Yield curves for r =0%. t dleiY PDE Krippner First‐order Second‐order Maturity (years) (d) Yield curves for r =1%. t Figure 1: Yield curves (zero-coupon yield against maturity in years) implied by a one- Q Q factor shadow-rate model with ρ = 0.01, ρ = 1, K = 0, K = −0.1, and Σ = 0.02. 0 1 0 1 The di(cid:27)erent panels correspond to di(cid:27)erent initial shadow short rates. Within each panel, the yield curve is computed using four methods: Numerical solution of PDE (7), Krippner’s (2012) approach described in Section 4.2, and the (cid:28)rst- and secondorder approximations proposed in Section 4.3. 17

short-rate paths per sample date.11,12 Table 1 shows the mean simulation error of yˆT, and the root-mean-square errors (RMSE) against yˆT of the yield yT computed by t t t ˙ Krippner’s (2012) method, the (cid:28)rst-order approximation y˜T de(cid:28)ned in (19), and the t second-order approximation y˜˜T de(cid:28)ned in (20). The table is divided into two panels. t The top panel shows errors for the sub-sample Jan 1990(cid:21)Nov 2008 (when interest rates were at normal levels), and the bottom panel shows errors for the sub-sample Dec 2008(cid:21)Dec 2012 (when the lower bound on nominal yields was binding at the short end of the yield curve). All methods are generally more precise at shorter maturities. As the (cid:28)rst column in both panels shows, the simulated yields are accurate to within approximately one (cid:28)fth of a basis point at the ten-year maturity point. As shown in the second column of the tables, Krippner’s (2012) method produces ten-year yields that are accurate to about one basis point during normal times, and to within four basis points when the lower bound is binding. While the (cid:28)rst-order method is more accurate when rates are low (the bottom panel), its errors remain substantial at the long end. The second-order method, the (cid:28)nal column in the tables, produces ten-year yields that are accurate to approximately half a basis point, both during normal times and when the lower bound is binding. Tofurtherillustratethetime-varyingrelativeperformanceofthethreeapproximation schemes, Figure 2 plots the di(cid:27)erence over time between the simulated ten-year yield and the three approximated yields. Krippner’s (2012) method and the secondorder approximation appear to be equally precise in the (cid:28)rst few years of the sample, 11I simulate (1) under Q based on moments (A.1)(cid:21)(A.2) on a uniformly-spaced grid with ∆t = ´ 1/360. For each simulated path i, I compute RT(i)= T r (i)ds by the trapezoidal method. I then t t s de(cid:28)nePˆT = 1 (cid:80)n exp(−RT(i))andyˆT =− 1 logPˆT. ThesimulationerrorforPˆT iscomputed t n i=1 t t T−t t t as n−1/2 times the sample standard deviation, and the simulation error for yˆT is derived from the t simulation error for PˆT by the delta method. t 12Thesimulationtakesseveralhourstocompleteforthegivenparameterestimatesandsmoothed states. This approach would not, therefore, be feasible as part of an estimation strategy. 18

Maturity se(yˆT) RMSE(yT) RMSE(y˜T) RMSE(y˜˜T) t t t t ˙ 6m 0.04 0.04 0.05 0.04 1y 0.06 0.06 0.18 0.06 2y 0.09 0.10 0.88 0.10 3y 0.12 0.14 2.26 0.13 4y 0.14 0.18 4.22 0.15 5y 0.16 0.27 6.60 0.17 7y 0.19 0.47 12.22 0.21 10y 0.21 0.93 21.81 0.35 (a) Sub-sample Jan 1990(cid:21)Nov 2008 Maturity se(yˆT) RMSE(yT) RMSE(y˜T) RMSE(y˜˜T) t t t t ˙ 6m 0.01 0.01 0.01 0.01 1y 0.02 0.04 0.04 0.02 2y 0.05 0.19 0.33 0.05 3y 0.07 0.51 1.07 0.07 4y 0.10 0.94 2.31 0.09 5y 0.12 1.42 3.99 0.12 7y 0.15 2.43 8.38 0.23 10y 0.17 3.87 16.63 0.52 (b) Sub-sample Dec 2008(cid:21)Dec 2012 Table 1: The mean standard error of simulated yields yˆT (n = 1,000,000 draws per t sampledate)forthemodelestimatedinAppendixB, andtheroot-mean-squareerrors (RMSE) against the simulated yields of Krippner’s (2012) yield approximation yT, t the (cid:28)rst-order yield approximation y˜T, and the second-order yield approximation y ˙˜˜T. t t All errors are in basis points. 19

with (cid:29)uctuations presumably largely due to simulation error.13 The discrepancy between simulated yield and Krippner’s (2012) method increases over time as the level of yields declines, and exceeds (cid:28)ve basis points by December 2012. The discrepancy between simulated yield and second-order approximation remains small and appears to show little systematic variation over time, perhaps trending up modestly towards the end of the sample. The (cid:28)rst-order approximation has a large negative discrepancy initially, which shrinks over time but remains at a high absolute level even at the end of the sample. Figure 2 also con(cid:28)rms that, just like in the simple one-factor model in Figure 1, the approximation errors under Krippner’s (2012) method and the (cid:28)rst-order scheme are largely systematic (rather than mere noise), in that the (cid:28)rst-order approximation overstates arbitrage-free yields while Krippner’s (2012) method tends to understate them. In sum, the analysis above suggests that, empirically, the second-order yield approximation y˜˜T is accurate to within about one half of a basis point at maturities up t to ten years, both during normal times and when the lower bound is binding. The approximation error is one order of magnitude smaller than both the model-implied observation error in yields (see Table 3) and the next best approximation method proposed by Krippner (2012). In contrast, the (cid:28)rst-order approximation is acceptable at most at the very short end of the yield curve. To add perspective, the approximation error in y˜˜T is no greater than commonly t accepted (cid:28)tting error in the derivation of constant-maturity zero-coupon bond yields from observed coupon-bearing Treasuries (e.g., G(cid:252)rkaynak et al., 2006). This puts the second-order approximation roughly on par with the numerical accuracy achieved 13Recall that when there is no lower bound, both methods produce yields equal to the exact arbitrage-free yield. In the early 1990s, the overall level of yields was su(cid:30)ciently high for the e(cid:27)ect of the lower bound to be negligible. 20

6 30 5 25 4 20 3 15 2 10 1 5 0 0 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 ‐1 ‐5 ‐2 ‐10 ‐3 ‐15 ‐4 ‐20 ‐5 ‐25 ‐6 ‐30 )stniop sisab( noitaiveD Krippner (left scale) Second‐order (left scale) First‐order (right scale) Figure 2: Deviation from simulated ten-year yield yˆt+10 of Krippner’s (2012) yield t approximation yt+10, the (cid:28)rst-order yield approximation y˜t+10, and the second-order t t yield approxima ˙ tion y˜˜t+10. Model parameters and (cid:28)ltered states are taken from Apt pendix B. 21

by (cid:16)exact(cid:17) bond pricing methods in standard a(cid:30)ne models, to the extent that they rely on numerical methods (say, to solve the system of ODEs (9)(cid:21)(10)).14 6 Conclusion This paper develops an approximation to arbitrage-free zero coupon bond yields in Gaussian shadow-rate term structure models. The complexity of the scheme does not dependonthenumberoffactors. Further, Idemonstratethatthemethodiscomputationally feasible by estimating a three-factor shadow-rate model of the U.S. Treasury yield curve. Based on Monte Carlo simulation, I also show that the yield approximation is approximately as precise as conventional approaches that are considered to be (cid:16)exact.(cid:17) 14Empirically, the second-order approximation y˜˜T is exact to an absolute tolerance of approxit mately 5×10−5 (see Table 1). The default tolerance for numerical methods in Matlab is typically between 10−4 and 10−6, depending on the complexity of the method. 22

A Useful Mathematical Results A.1 Moments of X t Consider the continuous time stochastic process X de(cid:28)ned in (1). The following t P Q derivations hold under both the -measure and the -measure, hence for notational simplicity I will suppress dependence of moments and parameters on the measure. Since X is a Gaussian process, all its (cid:28)nite-dimensional distributions are Gaussian t (Karatzas and Shreve, 1991). In particular, for u,s ≥ t, the vectors (X ,X ) are u s jointly conditionally Gaussian, with E [X ] = eK1(u−t)X +(I −eK1(u−t))K−1K (A.1) t u t N 1 0 ˆ u∧s Cov t [X u ,X s ] = eK1(u−v)ΣΣ(cid:62)eK 1 (cid:62)(s−v) dv. (A.2) t If K is invertible, the integral on the right-hand side of (A.2) can be evaluated 1 analytically using integration by parts and formula (10.2.15) in Hamilton (1994), vec(Cov [X ,X ]) (A.3) t u s = (K 1 ⊕K 1 )−1vec(eK1(u−t)ΣΣ(cid:62)eK 1 (cid:62)(s−t) −eK1(u−u∧s)ΣΣ(cid:62)eK 1 (cid:62)(s−u∧s)), where (cid:16)⊕(cid:17) denotes the Kronecker sum. Since r as de(cid:28)ned in (2) is a linear function t of the Gaussian random vector X , it is itself Gaussian, with t µ = E [r ] = ρ +ρ ·E [X ] (A.4) t→u t u 0 1 t u and σ = Cov [r ,r ] = ρ(cid:62)Cov [X ,X ]ρ . (A.5) t→u×s t u s 1 t u s 1 23

A.2 Moments of Censored Gaussian Random Variables This section derives two useful mathematical results involving the moments of censored Gaussian random variables. Lemma A.1. If X ∼ N(µ,σ2), then (cid:16)µ(cid:17) (cid:16)µ(cid:17) E[max{X,0}] = µΦ +σφ σ σ where Φ denotes the standard normal cdf, and φ denotes the standard normal pdf. Proof. First, note that (cid:20) (cid:26) (cid:27)(cid:21) X −µ µ E[max{X,0}] = µ+σE max ,− . (A.6) σ σ Thus,itonlyremainstocomputeE[max{Z,a}]whereZ isastandardnormalrandom variable, for arbitrary a ∈ R . By direct computation of the integral de(cid:28)ning the expectation, ˆ 1 ∞ (cid:18) 1 (cid:19) E[max{Z,a}] = aΦ(a)+ √ zexp − z2 dz, (A.7) 2π 2 a where Φ is the standard normal cdf. Further, ˆ ∞ (cid:18) 1 (cid:19) (cid:20) (cid:18) 1 (cid:19)(cid:21)∞ (cid:18) 1 (cid:19) √ zexp − z2 dz = −exp − z2 = exp − a2 = 2πφ(a). (A.8) 2 2 2 a a Recursively substituting from (A.8) into (A.7), and (cid:28)nally into (A.6), establishes the result. 24

Lemma A.2. If       X µ σ2 σ  1  ∼ N  1 , 1 12        X µ σ σ2 2 2 12 2 then E[max{X ,0}max{X ,0}] =(µ µ +σ )Φd(−ς ,−ς ;χ) 1 2 1 2 12 2 1 2 (cid:32) (cid:33) (cid:32) (cid:33) ς −χς ς −χς 1 2 2 1 +σ µ φ(ς )Φ +σ µ φ(ς )Φ 2 1 2 (cid:112) 1 2 1 (cid:112) 1−χ2 1−χ2 (cid:114) (cid:32)(cid:115) (cid:33) 1−χ2 ς2 −2χς ς +ς2 +σ σ φ 1 1 2 2 (A.9) 1 2 2π 1−χ2 where ς = µj, j ∈ {1,2}, χ = σ12 , φ denotes the univariate standard normal pdf, Φ j σj σ1σ2 denotes the univariate standard normal cdf, φ (z ,z ;χ) denotes the bivariate normal 2 1 2 pdfwhenboth variableshavezeromeans, unitvariances, andcorrelationχ, andΦ and 2 Φd denote the corresponding cumulative and decumulative bivariate Gaussian distri- 2 bution functions, where in particular Φd(z ,z ;χ) = 1−Φ(z )−Φ(z )+Φ (z ,z ;χ). 2 1 2 1 2 2 1 2 Proof. Write E[max{X ,0}max{X ,0}] 1 2 (cid:20) (cid:26) (cid:27) (cid:26) (cid:27)(cid:21) X −µ µ X −µ µ 1 1 1 2 2 2 =σ σ E max ,− max ,− 1 2 σ σ σ σ 1 1 2 2 +µ E[max{X ,0}]+µ E[max{X ,0}]−µ µ . (A.10) 2 1 1 2 1 2 The second and third terms can be evaluated using Lemma A.1. For the (cid:28)rst term, it su(cid:30)ces to be able to compute E[max{Z ,a}max{Z ,b}] for random variables Z 1 2 1 and Z that are bivariate normal with zero means, unit variances, and correlation 2 χ, and for arbitrary a,b ∈ R . Using the properties of φ , this expectation can be 2 25

expanded as follows: E[max{Z ,a}max{Z ,b}] 1 2 ˆ ˆ ∞ ∞ = max{z ,a}max{z ,b}φ (z ,z ;χ)dz dz 1 2 2 1 2 1 2 −∞ˆ −∞ˆ ˆ ˆ b a ∞ ∞ =ab φ (z ,z ;χ)dz dz +a z φ (z ,z ;−χ)dz dz 2 1 2 1 2 2 2 1 2 1 2 ˆ−∞ ˆ−∞ b ˆ −aˆ ∞ ∞ ∞ ∞ +b z φ (z ,z ;−χ)dz dz + z z φ (z ,z ;χ)dz dz . (A.11) 1 2 1 2 1 2 1 2 2 1 2 1 2 −b a b a The (cid:28)rst double integral is simply Φ (a,b;χ), the bivariate normal cdf. The second 2 andthirddoubleintegralscorrespondtoexpectedvaluesoftruncatedbivariatenormal random variables, and the last integral is the expected cross product of a truncated bivariate normal random vector. These expected values are known up to the univariate standard normal cdf and the bivariate normal cdf, respectively (see Rosenbaum, 1961). Using the formulas in Rosenbaum (1961) to evaluate the integrals in (A.11), and substituting into (A.10), we obtain (A.9) after simpli(cid:28)cation. 26

B Empirical Implementation This appendix empirically estimates a three-factor Gaussian shadow-rate term structure model using the yield approximation methodology proposed in Section 4.3. The main purpose is to demonstrate the computational tractability of the method in the contextofarealisticapplication. Formorein-depthdiscussionandempiricalanalysis, see Kim and Priebsch (2013). B.1 Data I use end-of-month zero-coupon U.S. Treasury yields from January 1990 through December 2012, for maturities of 6 months, 1(cid:21)5, 7, and 10 years. I derive the 6month yield from the corresponding T-bill quote, while longer-maturity zero yields are extracted from the CRSP U.S. Treasury Database using the unsmoothed Fama and Bliss (1987) methodology.15 I augment the yield data with survey forecasts from Blue Chip, interpolated to constant horizons of 1(cid:21)4 quarters (available monthly), as well as annually out to 5 years and for 5-to-10 years (available every six months).16 Model-implied survey forecasts are subject to the same lower-bound constraint as yields,17 but their computation is substantially simpler: Forecasters report their expectation of the arithmetic ´ (cid:104) (cid:105) mean of future observed short rates, EP 1 t+τ max{r ,0}ds . This is exactly (19) t τ t s P Q with the data-generating measure in place of pricing measure . Therefore, the (cid:28)rst-order method described in Section 4.3 produces exact model-implied survey forecasts. Intuitively, unlike yields, survey forecasts are not subject to compounding, so 15I am grateful to Anh Le for providing the code for this procedure. 16As discussed by Kim and Orphanides (2005), this potentially leads to more precise estimates of the parameters governing the data-generating distribution P. 17This follows from equivalence of the measures P and Q, and more fundamentally from the absence of arbitrage. 27

there are no higher-order Jensen’s inequality terms to consider. B.2 Filtering and Estimation Since the statistical properties of the term structure model laid out in Section 2 are formulatedintermsofthelatentstatevectorX ,butthedataactuallyobservedbythe t econometrician consist of yields, y , and survey expectations, z (see Appendix B.1), t t I set up a joint estimation and (cid:28)ltering problem to obtain estimates of the model’s parameters θ = (KP,KP,K Q ,K Q ,ρ ,ρ ,Σ).18 When discretely sampled at intervals 0 1 0 1 0 1 ∆t > 0, the state vector X follows a (cid:28)rst-order Gaussian vector autoregression, X = m +m X +ε (B.1) t+∆t 0,∆t 1,∆t t t where ε ∼ N(0,Ω ), and m , m , and Ω can be computed from (A.1) and t ∆t 0,∆t 1,∆t ∆t (A.2). Equation (B.1) represents the transition equation of the (cid:28)ltering problem. Next, denote by H : RN × Θ (cid:55)→ RMY the (non-linear) mapping from states y + X and parameters θ to model-implied yields y, and by H : RN × Θ (cid:55)→ RMZ the z + analogous mapping from states and parameters to model-implied survey forecasts z. For estimation purposes, I compute H through the second-order approximation (20), y and H through the exact (cid:28)rst-order method discussed in Appendix B.1. To simplify z notation, denote the stacked mapping (cid:0) H(cid:62),H(cid:62) (cid:1)(cid:62) by H. If we assume that all yields y z and survey expectations are observed with iid additive Gaussian errors, we obtain the observation equation   y t   = H(X )+e . (B.2)   t t z t 18See Duan and Simonato (1999) for an early reference discussing this approach towards term structure model estimation. 28

Together, equations (B.1) and (B.2) form a non-linear (cid:28)ltering problem. The simple (linear) Kalman (cid:28)lter(cid:22)optimal when measurement and observation equation are linear and all shocks are Gaussian(cid:22)has been modi(cid:28)ed in a number of ways to accommodate nonlinearity as in (B.2). The unscented Kalman (cid:28)lter, proposed by Julier et al. (1995), aims to deliver improved accuracy and numerical stability relative to the more traditional extended Kalman (cid:28)lter, without substantially increasing the computational burden.19,20 The algorithm is described in detail in Wan and van der Merwe (2001). As a by-product of the (cid:28)ltering procedure, it conveniently produces estimates of the mean and covariance matrix of (y ,z ) conditional t t on the econometrician’s information set as of time t − 1. I use these to set up a quasi(cid:21)maximum likelihood function based on (B.2),21 which I maximize numerically to obtain estimates of the parameters θ as well as their asymptotic standard errors (following Bollerslev and Wooldridge, 1992). B.3 Estimation Results To achieve econometric identi(cid:28)cation of the parameters θ in light of invariant transformations resulting in observationally equivalent models with di(cid:27)erent parameters (see Dai and Singleton, 2000), I follow Joslin et al. (2011) and impose the normalizations ρ = (1,...,1)(cid:62), K Q = 0, K Q is diagonal and therefore completely determined 1 0 1 by its ordered eigenvalues λQ , and Σ is lower triangular. I estimate the model on the data set described inAppendix B.1, using the quasi(cid:21) maximum likelihood (QML) procedure discussed in Appendix B.2. Table 2 displays ˆ the estimated model parameters θ, as well as their asymptotic standard errors. Table 3 shows the QML-estimated standard deviations of the measurement errors 19AdetailedtreatmentoftheunscentedKalman(cid:28)lter,andacomparisontotheextendedKalman (cid:28)lter, can be found in Wan and van der Merwe (2001). 20Christo(cid:27)ersen et al. (2012) and Wu (2010) con(cid:28)rm that the unscented Kalman (cid:28)lter performs 29

ρ 0.0738 r 0.0010 0 min (0.0043) (0.0001) λQ −0.1038 Σ 0.0268 (0.0226) (0.0084) −0.3566 −0.0324 0.0416 (0.1177) (0.0110) (0.0295) −0.8574 0.0068 −0.0397 0.0090 (0.2876) (0.0100) (0.0302) (0.0007) KP −0.0193 KP −0.4679 −0.3415 0.3785 0 1 (0.0052) (0.1574) (0.1490) (0.6798) −0.0099 −0.5752 −1.1881 −1.1875 (0.0224) (0.6303) (1.1335) (0.5740) 0.0278 0.8908 1.3060 0.3990 (0.0233) (0.6982) (1.0641) (1.3561) Table 2: Quasi(cid:21)maximum likelihood parameter estimates (asymptotic standard errors) for the three-factor shadow-rate model. in yields and survey variables (e in equation (B.2)). The average yield error is 8 t basis points, and the average error in surveys is 21 basis points. For both yields and surveys, errors follow a U-shaped pattern, being largest at the short and long ends. Figure 3 plots the model-implied shadow short rate r over the sample period, t based on the states implied by the Kalman smoother (that is, incorporating all information up to December 2012, the end of the sample). The shadow rate turned negative in December 2008, after the FOMC established a target federal funds rate range of 0 to 0.25 percent and the e(cid:27)ective lower bound became binding, and has stayed negative through the end of the sample. better than the extended Kalman (cid:28)lter in the speci(cid:28)c setting of term structure model estimation. 21This estimation approach is described and analyzed in Lund (1997). 30

Maturity σ Maturity σ Y Z 6m 0.0017 1q 0.0014 1y 0.0014 2q 0.0002 2y 0.0006 3q 0.0009 3y 0.0003 4q 0.0014 4y 0.0004 2y 0.0028 5y 0.0003 3y 0.0026 7y 0.0006 4y 0.0027 10y 0.0015 5y 0.0031 Average 0.0008 5y(cid:21)10y 0.0034 Average 0.0021 Table3: Estimatedstandarddeviationsofobservationerrorsinyields, σ , andsurvey Y forecasts, σ . Z 8% 6% 4% 2% 0% 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 ‐2% ‐4% Figure 3: Model-implied shadow short rate r based on smoothed states X . t t|T 31

References Bauer, M., and Rudebusch, G. (2013), (cid:16)Monetary Policy Expectations at the Zero Lower Bound,(cid:17) Working Paper, Federal Reserve Bank of San Francisco Bernanke, B., Reinhart, V., and Sack, B. (2004), (cid:16)Monetary Policy AlternativesattheZeroBound: AnEmpiricalAssessment,(cid:17) Brookings Papers on Economic Activity, 2:1(cid:21)100 Black, F. (1995), (cid:16)Interest Rates as Options,(cid:17) Journal of Finance, 50(5):1371(cid:21)1376 Bollerslev, T., and Wooldridge, J. (1992), (cid:16)Quasi(cid:21)Maximum Likelihood Estimation and Inference in Dynamic Models with Time-Varying Covariances,(cid:17) Econometric Reviews, 11(2):143(cid:21)172 Chen, R. (1995), (cid:16)A Two-Factor, Preference-Free Model for Interest Rate Sensivite Claims,(cid:17) Journal of Futures Markets, 15(3):345(cid:21)372 Christensen, J., and Rudebusch, G. (2013), (cid:16)Estimating Shadow-Rate Term Structure Models with Near-Zero Yields,(cid:17) Working Paper, Federal Reserve Bank of San Francisco Christo(cid:27)ersen, P., Dorion, C., Jacobs, K., and Karoui, L. (2012), (cid:16)Nonlinear Kalman Filtering in A(cid:30)ne Term Structure Models,(cid:17) CREATES Research Paper 2012-49, Aarhus University Dai, Q., and Singleton, K. (2000), (cid:16)Speci(cid:28)cation Analysis of A(cid:30)ne Term Structure Models,(cid:17) Journal of Finance, 60(5):1943(cid:21)1978 Duan, J.-C., and Simonato, J.-G. (1999), (cid:16)Estimating and Testing Exponential- 32

A(cid:30)ne Term Structure Models by Kalman Filter,(cid:17) Review of Quantitative Finance and Accounting, 13:111(cid:21)135 Du(cid:30)e, D., and Kan, R. (1996), (cid:16)A Yield-Factor Model of Interest Rates,(cid:17) Mathematical Finance, 6:369(cid:21)406 Fama, E., and Bliss, R. (1987), (cid:16)The Information in Long-Maturity Forward Rates,(cid:17) American Economic Review, 77(4):680(cid:21)692 Genz, A. (2004), (cid:16)Numerical Computation of Rectangular Bivariate and Trivariate Normal and t Probabilities,(cid:17) Statistics and Computing, 14:251(cid:21)260 Gorovoi, V., and Linetsky, V. (2004), (cid:16)Black’s Model of Interest Rates as Options, Eigenfunction Expansions and Japanese Interest Rates,(cid:17) Mathematical Finance, 14(1):49(cid:21)78 Hamilton, J. (1994), The Time Series Analysis, Princeton University Press Ichiue, H., and Ueno, Y. (2007), (cid:16)Equilibrium Interest Rate and the Yield Curve in a Low Interest Environment,(cid:17) Bank of Japan Working Paper (2013), (cid:16)Estimating Term Premia at the Zero Bound: An Analysis of Japanese, U.S., and U.K. Yields,(cid:17) Bank of Japan Working Paper No. 13(cid:21)E(cid:21)8 Joslin, S., Singleton, K., and Zhu, H. (2011), (cid:16)A New Perspective on Gaussian Dynamic Term Structure Models,(cid:17) Review of Financial Studies, 24(3):926(cid:21)970 Julier, S., Uhlmann, J., and Durrant-Whyte, H. (1995), (cid:16)A New Approach for Filtering Nonlinear Systems,(cid:17) in Proceedings of the American Control Conference, pp. 1628(cid:21)1632 33

Karatzas, I., and Shreve, S. (1991), Brownian Motion and Stochastic Calculus, Graduate Texts in Mathematics Series, Springer, London Kim, D., and Orphanides, A. (2005), (cid:16)Term Structure Estimation with Survey Data on Interest Rate Forecasts,(cid:17) Sta(cid:27) Working Paper 2005-48, Federal Reserve Board Kim, D., and Priebsch, M. (2013), (cid:16)Estimation of Multi-Factor Shadow Rate Models,(cid:17) Working Paper in Progress, Federal Reserve Board, Washington D.C. Kim, D., and Singleton, K. (2012), (cid:16)Term Structure Models and the Zero Bound: An Empirical Investigation of Japanese Yields,(cid:17) Journal of Econometrics, 170(1):32(cid:21)49 Krippner, L. (2012), (cid:16)Modifying Gaussian Term Structure Models When Interest Rates are Near the Zero Lower Bound,(cid:17) Reserve Bank of New Zealand Discussion Paper 2012/02 Lund, J.(1997),(cid:16)Non-LinearKalmanFilteringTechniquesforTermStructureModels,(cid:17) Working Paper, Aarhus School of Business Piazzesi, M. (2010), (cid:16)A(cid:30)ne Term Structure Models,(cid:17) in Handbook of Financial Econometrics, ed.byY.A(cid:239)t-Sahalia, andL.P.Hansen.North-Holland, Amsterdam Rosenbaum, S. (1961), (cid:16)Moments of a Truncated Bivariate Normal Distribution,(cid:17) Journal of the Royal Statistical Society, 23(2):405(cid:21)408 Severini, T. (2005), Elements of Distribution Theory, Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press Wan, E., and van der Merwe, R. (2001), (cid:16)The Unscented Kalman Filter,(cid:17) in Kalman Filtering and Neural Networks, ed. by S. Haykin, pp. 221(cid:21)280 34

Wu, S. (2010), (cid:16)Non-Linear Filtering in the Estimation of a Term Structure Model of Interest Rates,(cid:17) WSEAS Transactions on Systems, 9(7):724(cid:21)733 35