Contents lists available at ScienceDirect

Journal of Econometrics journal homepage: www.elsevier.com/locate/jeconom

Estimation of affine term structure models with spanned or unspanned stochastic volatility✩ Drew D. Creal a,∗ , Jing Cynthia Wu a,b a

Chicago Booth, United States

b

NBER, United States

article

info

Article history: Received 3 October 2013 Received in revised form 2 October 2014 Accepted 5 October 2014 Available online 4 November 2014 JEL classification: C32 E32 E43

abstract We develop new procedures for maximum likelihood estimation of affine term structure models with spanned or unspanned stochastic volatility. Our approach uses linear regression to reduce the dimension of the numerical optimization problem yet it produces the same estimator as maximizing the likelihood. It improves the numerical behavior of estimation by eliminating parameters from the objective function that cause problems for conventional methods. We find that spanned models capture the cross-section of yields well but not volatility while unspanned models fit volatility at the expense of fitting the crosssection. © 2014 Elsevier B.V. All rights reserved.

Keywords: Affine term structure models Unspanned stochastic volatility Estimation

1. Introduction We propose new estimation procedures for affine term structure models (ATSMs) with spanned or unspanned stochastic volatility that use linear regression to simplify and stabilize estimation. For spanned models, our procedure recovers the maximum likelihood estimator but only requires numerically optimizing over a lower dimensional parameter space. The stability of our method makes it possible for us to study local maxima, explain why they exist, and their economic implications. We show how our insights from spanned models can be extended to estimate unspanned stochastic volatility (USV) models despite the fact that for

✩ We thank Yacine Ait-Sahalia, Boragan Aruoba, Michael Bauer, Alan Bester, John Cochrane, Frank Diebold, Rob Engle, Jim Hamilton, Chris Hansen, Guido Kuersteiner, Ken Singleton, two anonymous referees, and seminar and conference participants at Chicago Booth, NYU Stern, NBER Summer Institute, Maryland, Bank of Canada, Kansas, UMass, and Chicago Booth Junior Finance Symposium for helpful comments. Drew Creal thanks the William Ladany Faculty Scholar Fund at the University of Chicago Booth School of Business for financial support. Cynthia Wu also gratefully acknowledges financial support from the IBM Faculty Research Fund at the University of Chicago Booth School of Business. This paper was formerly titled ‘‘Estimation of non-Gaussian affine term structure models’’. ∗ Corresponding author. E-mail addresses: [email protected] (D.D. Creal), [email protected] (J.C. Wu).

http://dx.doi.org/10.1016/j.jeconom.2014.10.003 0304-4076/© 2014 Elsevier B.V. All rights reserved.

USV models the likelihood function is not known in closed-form. Estimating a range of popular models, we find that models with spanned volatility fit the cross section of the yield curve better, while those with unspanned volatility fit the volatility better. ATSMs are popular among policy makers, practitioners, and academic researchers for studying bond prices, monetary policy, and the macroeconomic determinants of discount rates; for overviews, see Piazzesi (2010), Duffee (2012), Gürkaynak and Wright (2012), and Diebold and Rudebusch (2013). As the literature on ATSMs has developed over the last decade, there is a consensus that estimation can be challenging; see, e.g. Duffee (2002), Ang and Piazzesi (2003), Kim and Orphanides (2005), and Hamilton and Wu (2012). New procedures for Gaussian ATSMs have made them easier to estimate, further increasing their popularity; see, Joslin et al. (2011), Christensen et al. (2011), Hamilton and Wu (2012), Adrian et al. (2012) and Diez de Los Rios (2013). However, these procedures do not address models with stochastic volatility. Moreover, in USV models as proposed by CollinDufresne and Goldstein (2002) and Collin-Dufresne et al. (2009), the likelihood function is not known in closed-form. Potential solutions to this problem are the closed-form expansions of the likelihood for continuous-time models developed by Aït-Sahalia (2008) and Aït-Sahalia and Kimmel (2010) and the expectation maximization (EM) algorithm. Combining our approach with these could

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

potentially improve estimation; we demonstrate this for the EM algorithm explicitly. Our main contribution are new procedures for estimating ATSMs with spanned or unspanned stochastic volatility. For models with spanned factors where volatility factors price bonds, we propose to maximize a concentrated likelihood that when optimized gives exactly the same estimator as maximizing the original likelihood function. However, it only requires numerically optimizing over a subset of the parameters. The concentrated likelihood function is simple to construct from linear regressions. Using this approach, estimation of spanned models only takes a fraction of a second to several minutes compared to hours when optimizing the original likelihood. For USV models where the volatility factors do not price bonds, the log-likelihood function is not known in closed-form adding another layer of difficulty. Nevertheless, we show how the intuition behind the concentrated likelihood for spanned models can be extended to estimate USV models using the EM algorithm of Dempster et al. (1977). The maximization step of the EM algorithm solves a similar problem as optimizing the likelihood function of a spanned Gaussian ATSM. Consequently, we can construct a concentrated objective function for the EM algorithm using linear regressions just as we did for spanned models. Our method outperforms conventional approaches both in terms of stability of convergence and speed. A study for a 3-factor model with one spanned volatility factor shows that our method guarantees convergence as long as it is locally identified, and it converges to a number of local maxima repeatedly. Aside from being able to find the global maximum, our method helps us to locate and understand the economic implications of different local maxima. Conversely, the conventional method of directly maximizing the original likelihood never converges fully to any of the local maxima, nor does it converge to the same point twice in repeated trials even when it is initialized under the same local mode. This makes it difficult for researchers to differentiate between points near a well-behaved local maximum having the same economic meaning and locations corresponding to local maxima that are economically different. The median time it takes for our new procedure is less than 2 minutes for this model, whereas the conventional approach takes over 2 hours. Using our method, we shed light on how local maxima with different economic implications are created in non-Gaussian spanned models. In Gaussian models, different rotations of the factors (such as re-ordering of the factors) result in equivalent global maxima, with identical economic implications. In non-Gaussian models with spanned factors, rotations can have substantial economic impacts. The non-Gaussian state variables must be positive and enter the conditional variance. This creates an asymmetry between the Gaussian and non-Gaussian factors resulting in local maxima that are not economically equivalent. Another contribution of this paper is to develop a family of discrete-time non-Gaussian ATSMs that encompasses continuoustime models, including both spanned models as in Duffie and Kan (1996), Duffee (2002), Cheridito et al. (2007), and Aït-Sahalia and Kimmel (2010) as well as USV models as proposed by CollinDufresne and Goldstein (2002). Gouriéroux et al. (2002) proposed a one factor discrete-time non-Gaussian model and Le et al. (2010) generalized it to have multiple factors. Our model encompasses any admissible rotation of a multivariate discrete-time Cox et al. (1985) process, allowing the factors to be correlated. The model nests the risk-neutral dynamics of other discrete-time ATSMs.1 In our model, the physical and risk neutral dynamics follow the

1 In this paper, we do not consider the class of non-Gaussian ATSMs built from the non-central Wishart process of Gouriéroux et al. (2009).

61

same stochastic process but with different parameters. The market prices of risk have the extended affine form of Cheridito et al. (2007), which is different than Le et al. (2010). Finally, we also provide the restrictions needed to generate USV in discrete-time versions of the continuous-time models studied by Collin-Dufresne et al. (2009) and Joslin (2010). We apply our estimation method to a range of popular spanned and unspanned models with three and four factors. Judging by the estimated likelihood, a model with three spanned non-Gaussian factors has the highest likelihood followed by one of the USV models. Gaussian and non-Gaussian models with spanned factors fit the cross-section of yields equally well. However, spanned models do not capture the volatility well at any maturity, even for the best fitting model. This is because the non-Gaussian state variables must simultaneously fit the conditional mean and variance. Maximum likelihood places more weight on the first moment. In order to guarantee unspanned volatility factors, USV models place restrictions on the bond loadings. This causes USV models to sacrifice some cross-sectional fit; their pricing errors are larger than spanned models. On the other hand, USV models fit the dynamics of yield curve volatility well. The USV restrictions are not unique and we show that the choice of which USV restrictions are imposed is not inconsequential. This paper continues as follows. In Section 2, we specify a general class of discrete-time, non-Gaussian affine term structure models. In Section 3, we describe our new approach to estimation for both spanned and unspanned models. Section 4 describes the data and parameter restrictions of the models. Section 5 studies a three factor spanned model in depth. In Section 6, we study eight three and four factor spanned and unspanned models. In Section 7, we discuss directions for future research and conclude. 2. Model In this section, we describe a class of discrete-time ATSMs with stochastic volatility that encompass both spanned models, as in Duffie and Kan (1996), Dai and Singleton (2000), Cheridito et al. (2007); and unspanned models, as proposed by Collin-Dufresne and Goldstein (2002). 2.1. Bond prices The model has a G × 1 vector of conditionally Gaussian state variables gt , whose volatilities are captured by an H × 1 vector of positive state variables ht . Under the risk-neutral measure Q, the Gaussian state variables follow a vector autoregression with conditional heteroskedasticity Q gt +1 = µQ g + Φg gt + Φgh ht + Σgh εh,t +1 + εg ,t +1 , Q

Q

Q

Q εgQ,t +1 ∼ N 0, Σg ,t Σg′ ,t ,

Σg ,t Σg′ ,t = Σ0,g Σ0′ ,g +

H

(1)

Σi,g Σi′,g hit ,

i=1

εhQ,t +1 = ht +1 − EQ (ht +1 |It ) where It captures agents’ information set at time t. The volatility factors ht are an affine transformation of the exact discrete-time equivalent of a multivariate Cox et al. (1985) process ht +1 = µh + Σh wt +1

(2)

wi,t +1 ∼ Gamma νhQ,i + ziQ,t +1 , 1 , i = 1, . . . , H Q Q zi,t +1 ∼ Poisson e′i Σh−1 Φh Σh wt , i = 1, . . . , H

(3) (4)

where ei denotes the ith column of the identity matrix IH . We discuss the admissibility restrictions and interpretation of the parameters of the model in Section 2.2.

62

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

The price of a zero-coupon bond with maturity n at time t is the expected price of the same asset at time t + 1 discounted by the short rate rt under the risk neutral measure Q

Ptn = Et

exp (−rt ) Ptn+−11 .

rt = δ0 + δ1′ ,h ht + δ1′ ,g gt .

Σg ,t Σg′ ,t = Σ0,g Σ0′ ,g +

Given the dynamics of gt and ht under Q, bond prices are an exponentially affine function of the state variables

= exp a¯ n + b¯ ′n,h ht + b¯ ′n,g gt .

The loadings a¯ n , b¯ n,h and b¯ n,g can be expressed recursively in matrix notation as Q ′

′¯ a¯ n = −δ0 + a¯ n−1 + µQ g bn−1,g + µh − Φh µh + Σh νh Q

b¯ n−1,h

1 + b¯ ′n−1,g Σ0,g Σ0′ ,g b¯ n−1,g 2 −1 ′ Σh b¯ n−1,gh + µ′h ΦhQ′ Σh−1′ IH − diag ιH − Σh′ b¯ n−1,gh − νhQ′ log ιH − Σh′ b¯ n−1,gh + Σh′ b¯ n−1,gh Q′ Q′ b¯ n,h = −δ1,h + Φgh b¯ n−1,g + Φh b¯ n−1,h

+

1 2

IH ⊗ b¯ ′n−1,g Σg Σg′ ιH ⊗ b¯ n−1,g

(5)

(6)

b¯ n,g = −δ1,g + ΦgQ′ b¯ n−1,g

(7)

with initial values a¯ 1 = −δ0 , b¯ 1,g = −δ1,g and b¯ 1,h = −δ1,h , see Appendix B for a derivation. The matrix Σg Σg′ is a GH × GH block diagonal matrix with elements Σi,g Σi′,g for i = 1, . . . , H and ′ ¯ bn−1,g + b¯ n−1,h . The loadings must satisfy the restricb¯ n−1,gh = Σgh

tion that the ith component of Σh′ b¯ n−1,gh < 1 for i = 1, . . . , H. Bond yields ynt ≡ − 1n log Ptn are linear in the factors ′

′

= an + bn,h ht + bn,g gt

(8)

with an = − ¯ bn,h = − ¯ and bn,g = − ¯ Gouriéroux and Jasiak (2006) built the univariate version of the non-Gaussian process (2)–(4) and Gouriéroux et al. (2002) used it to construct a one factor ATSM. Le et al. (2010) extended the process to allow for multiple factors; their specification under Q is (2)–(4) but with µh = 0 and Σh diagonal. 1 a , n n

1 b n n,h

1 b . n n ,g

2.1.1. Unspanned stochastic volatility models USV models (see Collin-Dufresne and Goldstein (2002)) impose restrictions on the parameters under Q guaranteeing that the bond loadings for the volatility factors bn,h in (8) are zero for all maturities. Yields consequently only depend on the Gaussian factors ynt = an + b′n,g gt . The bond loadings (5)–(7) simplify to ′¯ a¯ n = −δ0 + a¯ n−1 + µQ g bn−1,g +

b¯ n,g = −δ1,g + ΦgQ′ b¯ n−1,g

1 ′ b¯ n−1,g Σ0,g Σ0′ ,g b¯ n−1,g

2

H

(11)

Σi,g Σi′,g hit ,

i =1

εh,t +1 = ht +1 − E (ht +1 |It ) . The Gaussian state variables are a function of the non-Gaussian state variables through both the autoregressive term Φgh ht and the covariance term Σgh εh,t +1 . The parameters controlling the conditional mean are different under P and Q measures, while the scale parameters Σgh and Σi,g for i = 0, . . . , H are the same. The model for ht +1 under the physical measure P is h t + 1 = µ h + Σh w t + 1

−1 ′ − ΦhQ′ Σh−1′ IH − diag ιH − Σh′ b¯ n−1,gh Σh b¯ n−1,gh

ynt

gt +1 = µg + Φg gt + Φgh ht + Σgh εh,t +1 + εg ,t +1 ,

εg ,t +1 ∼ N 0, Σg ,t Σg′ ,t ,

The short rate is a linear function of the state vector

Ptn

under Q. The Gaussian state variables follow a vector autoregression with conditional heteroskedasticity

(9) (10)

which are the same as Gaussian ATSMs. Unlike Gaussian ATSMs, however, USV models constrain some of the Q parameters Q (i.e. elements of Φg ) in order to set the non-Gaussian loadings to zero. In Section 4.2.2, we provide conditions under which discretetime ATSMs exhibit USV as in Collin-Dufresne et al. (2009). 2.2. Physical dynamics

lower bound. A similar set of restrictions must be satisfied under Q. The scale parameters Σh are the same under both probability measures and so is the parameter µh . The latter restriction is required for no-arbitrage. The conditional mean of the volatility factors ht +1 can be written in matrix form as E (ht +1 |It ) = (IH − Φh ) µh + Σh νh + Φh ht .

It is a linear function of its own lag ht , similar to a vector autoregression. The conditional variance is also an affine function of ht V (ht +1 |It ) = Σh diag(νh − 2Σh Φh µh )Σh′ −1

+ Σh diag 2Σh−1 Φh ht Σh′ .

In Appendix A.2, we provide the transition density of ht +1 for any admissible rotation. ′ A nice property of the model (11)–(14) for the vector h′t , gt′ is that any admissible affine transformation remains within the same family of distributions. Proposition 1. Let gt and ht follow the process of (11)–(14) with parameters θ . Consider an admissible affine transformation of the form

c C h˜ t = h + hh cg Cgh g˜t

Chg Cgg

ht gt

.

The new processes g˜t and h˜ t remain in the same family of distributions under updated parameters θ˜ . The parameters νh and Σh−1 Φh Σh are invariant to rotation. Proof. See Appendix C.1. The admissibility restrictions and the relationship between the new and old parameterizations can be found in Appendix C.1. This proposition helps to understand identification in Section 4.2.1. The admissibility constraints ensure that the non-Gaussian state variables always remain positive after applying a transformation from h′t , gt′

Analogous to the popular class of Gaussian ATSMs, we specify the dynamics of gt and ht under P to have the same dynamics as

(12)

wi,t +1 ∼ Gamma νh,i + zi,t +1 , 1 , i = 1, . . . , H (13) ′ −1 zi,t +1 ∼ Poisson ei Σh Φh Σh wt , i = 1, . . . , H (14) where νh = νh,1 , . . . , νh,H are shape parameters, Φh is a matrix controlling the autocorrelation of ht +1 , Σh is a scale matrix, and µh is a vector determining the lower bound of ht +1 . Sufficient conditions for non-negativity of ht are that elements of µh , Σh , and Σh−1 Φh Σh are non-negative. The discrete-time equivalent of the Feller condition νh,i > 1 ensures that the process does not attain its

′

to h˜ ′t , h˜ ′t

′

and that there exists another admissible

rotation to get back to the original factors.

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

63

2.3. Stochastic discount factor

3. Estimation methodology

In this section, we demonstrate how an agent gets compensated for risk exposure when holding a zero-coupon bond under stochastic volatility. Given the dynamics of the state vector under P and Q measures, the market prices of risk have the extended affine form of Cheridito et al. (2007). The full expression is in Appendix D. To provide intuition, the log of the stochastic discount factor (SDF) can be decomposed up to a first order approximation into the risk free rate plus three components describing risk compensation

In this section, we introduce new estimation procedures for spanned and unspanned models, both of which use least square regressions to simplify and stabilize estimation. Our approach is based on the following observations: (i) The parameter vector θ can be separated into those parameters that enter the bond loadings and those that do not (e.g. µg , Φg , Φgh ); (ii) Given the parameters that enter the loadings, we can calculatethe bond loadings and

mt +1 = −rt −

1 ′ λgt λgt − λ′gt ϵg ,t +1 − λ′wt ϵw,t +1 − λ′zt ϵz ,t +1 2

(15)

where ϵi,t +1 are standardized shocks with mean zero and identity covariance matrix, and λit is the price of risk i for each of the three types of shocks in the model.2 In addition to the risk-free rt , the agent gets compensated for being exposed to the Gaussian shock ϵg ,t +1 in (11), the gamma shock ϵw,t +1 in (13), and the Poisson shock ϵz ,t +1 in (14). The prices of these risks are defined as

λgt = V(g t +1 |It , ht +1 , zt +1 )−1/2 × E (gt +1 |It , ht +1 , zt +1 ) − EQ (gt +1 |It , ht +1 , zt +1 )

(1)

1 solve for the spanned factors xt = B− Yt 1

(iii) The P parameters (e.g. µg , Φg , Φgh ) of the Gaussian VAR for the factors gt plus Ω enter the objective function as a quadratic form. The first order conditions for these parameters (µg , Φg , Φgh , Ω ) are linear and can be solved by running (generalized) least squares regressions of (11). Using this basic insight, we show how to eliminate these parameters from the numerical optimization problem. 3.1. Spanned models Given the parameters of the model θ , the likelihood function is

(2)

(1)

(1)

p (Y1:T ; θ) = p Y1:T |Y1:T ; θ p Y1:T ; θ

−1/2 λwt = V(w t +1 |It , zt +1 ) × E(wt +1 |It , zt +1 ) − EQ (wt +1 |It , zt +1 ) , λzt = V(zt +1 |It )−1/2 E(zt +1 |It ) − EQ (zt |It ) .

=

T

(2)

(1)

p Yt |Yt ; θ |J (θ )|−T

t =1

The market prices of risk have an intuitive form as the Sharpe ratio measuring per unit risk compensation. Specifically, they are the difference in the conditional means of each shock under P and Q standardized by a conditional standard deviation. The timevarying quantities of risk are a feature of non-Gaussian models that are not available in Gaussian models. USV models In USV models, the components of the SDF associated with the Gaussian factors (the first three terms in (15)) are the only parts that are directly observable from bond yields. The risk premium for non-Gaussian factors can only be estimated jointly by also observing derivatives because the Q parameters of the volatility process do not enter the likelihood and are unidentified by observing only yields.

×

T

p (gt |ht , It −1 ; θ)

t =1

T H

p (hit |It −1 ; θ )

where J (θ ) is the Jacobian of the transformation from xt = (gt′ , h′t )′ (1)

to Yt , see Appendix E for the log-likelihood ℓ (θ) = log p(Y1:T ; θ ).3 Direct maximization of the log-likelihood is however extremely challenging as interest rates are close to non-stationary, the bond loadings are non-linear functions of the models’ parameters, and the maximization must impose the condition that ht > 0. Our approach to spanned models is a result of the following proposition: Proposition 2. If the model is given by Eq. (11)–(14) and (16)–(17) with all spanned factors, then the maximum likelihood estimator θˆ = argmaxθ ℓ (θ) can be solved by maximizing the concentrated

likelihood maxθm ℓ θˆc (θm ) , θm , where θc = µg , Φg , Φgh , Ω and

2.4. State space representation Define xt as N1 × 1 vector of spanned factors; i.e., xt = (gt′ , h′t )′ in spanned models and xt = gt in USV models. Stacking ynt in order for N different maturities n1 , n2 , . . . , nN gives Yt = A + Bxt where A = (an1 , . . . , anN )′ , B = (b′n1 , . . . , b′nN )′ . If more yields are observed than the number of spanned factors (N > N1 ), not all yields can be priced exactly. Following the literature, we assume (1) that N1 linear combinations of the yields Yt = SY1 Yt are priced without error and the remaining N2 = N − N1 linear combina(2) tions Yt = SY2 Yt are observed with Gaussian measurement errors. Given this assumption, the observation equations are (1) (2)

Yt

= A1 + B1 xt

(16)

= A2 + B2 xt + ηt ηt ∼ N (0, Ω )

(17)

where A1 ≡ SY1 A, A2 ≡ SY2 A, B1 ≡ SY1 B, and B2 ≡ SY2 B. The state space representation of the model is completed using the dynamics of the state variables (11)–(14).

2 This approximation is not used during estimation.

(18)

t =1 i=1

Yt

− A1 using (16);

θm are the remaining parameters of the model. The function θˆc (θm ) is obtained by solving maxθc ℓ (θm , θc ) using the GLS estimates for the P dynamics (11) and the OLS estimates for the variance–covariance matrix Ω in (17). Proof. See Appendix F.1. The proposition raises two points. First, optimizing the concentrated likelihood gives exactly the same solution as maximizing the original likelihood function in (18). However, it only requires numerically optimizing over θm instead of both θm and θc . Second, the concentrated likelihood function can be constructed from linear regressions. The method we propose is an immediate result of Proposition 2.

3 The stationary distribution is only known for special sub-classes of the affine family of models. In this paper, we assume a diffuse initial condition and start from t = 2. This provides an analytical solution for the first order conditions of the likelihood. If a researcher wants to include the stationary distribution as the initial condition, we recommend using our procedure first, and then using these estimates as starting values to optimize the likelihood with the initial condition. While including the initial conditions enforces stationarity, it can also introduce a downward bias in estimates of the autoregressive parameters; see, e.g. Bauer et al. (2012).

64

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

Procedure the concentrated log-likelihood function 1. Maximize

maxθm

ℓ θˆc (θm ) , θm . For a given value of θm , the concentrated

likelihood can be constructed as follows:

p

(i.) Given θm , calculate the bond loadings A and B and the state (1) 1 variables gt and ht from xt = B− Yt − A1 . 1

gt +1 − Σgh εh,t +1 = µg + Φg gt + Φgh ht + Σg ,t εg ,t +1

(19)

ˆ gh (θm ). ˆ g (θm ), Φ to calculate µ ˆ g (θm ) , Φ (iii.) Calculate the covariance matrix T

1

T − 1 t =2

(2)

− A2 − B2 xt

Yt

′

(2)

(20) − A2 − B2 xt . ˆ g (θm ) , Φ ˆ gh (θm ) , Ω ˆ (θm ) (iv.) Substitute θˆc (θm ) = µ ˆ g (θm ) , Φ

× Yt

back into the original likelihood to form the concentrated likelihood. In Appendix F, we also derive the analytical gradients of the concentrated log-likelihood. Our derivation is based on the following proposition. It decomposes the gradient into pieces according to whether a parameter enters the bond loadings, the P dynamics, or both. Proposition 3. The gradient of the concentrated log-likelihood ℓ θˆc (θm ) , θm can be decomposed into three terms:

dℓ θˆc (θm ) , θm , A (θm ) , B (θm ) dθm′ ∂ℓ θˆc , θm , A, B

∂ℓ θˆc , θm , A, B ∂ A (θ ) m + = ∂θm′ ∂ A′ ∂θm′ ∂ℓ θˆc , θm , A, B ∂ vec B (θm )′ + . ∂ vec (B′ )′ ∂θm′ The first term is the partial derivative of the P dynamics and Jacobian with respect to θm . This measures the effect parameters have on the log-likelihood through the time series of the factors. The second and third terms measure the effect parameters have on the log-likelihood through the bond loadings A and B. Proof. See Appendix F.2 The expressions for the gradient can be used for other affine models such as models for defaultable bonds and credit default swaps. 3.2. Unspanned models In a USV model, the pricing equation (16) can be inverted to calculate the Gaussian factors gt conditional on the parameters that enter the bond loadings. However, the volatility factors ht cannot be observed by inverting The likelihood of the the pricing formula. (2)

(1)

(1)

model p (Y1:T ; θ ) = p Y1:T |Y1:T ; θ p Y1:T ; θ is no longer known

(2)

(1)

in closed-form. The first term p Y1:T |Y1:T ; θ as in (18). The second term p

(1) Y1:T ; θ

(1) Y1:T ; θ

= |J (θ )|

−T

...

p (g1:T |h0:T −1 ; θ )

× p (h0:T −1 ; θ) dh0 . . . dhT −1 (1)

(ii.) Given gt and ht , calculate εh,t +1 and Σg ,t . Run a GLS regression

ˆ (θm ) = Ω

dynamics of the factors and is an integral over the path of the latent volatility

remains the same

is associated with the P

where J (θ) is the Jacobian from gt to Yt . This integral does not have a closed-form solution. We use the Monte Carlo Expectation Maximization (MCEM) algorithm to estimate the model, see Wei and Tanner (1990). The EM algorithm consists of two steps: the expectation and maximization steps, which are iterated back and forth until convergence of the algorithm to a stationary point of the likelihood. The first step calculates the expected value of the complete data log-likelihood

T (i) (2) =E log p Yt |gt ; θ − T log |J (θ )| Q θ |θ t =1

+

T

log p (gt |gt −1 , ht −1 ; θ )

t =1

+

T −1

log p (ht |ht −1 ; θ ) + p (h0 ; θ) .

(21)

t =1

This expectation is taken with respect to the posterior distribution p h0:T −1 |Y1:T ; θ (i) , which depends on the parameters θ (i) from

the previous iteration. The function Q θ |θ (i) is known as the intermediate quantity of the EM algorithm and it is a function of θ . In the second step of the EM algorithm, the intermediate quantity is maximized θ (i+1) = argmaxθ Q θ |θ (i) to determine the parameters for the next iteration. For USV models, the intermediate quantity has the same form as the log-likelihood for spanned Gaussian models. This means that maximization of the intermediate quantity at each iteration of the EM algorithm is (essentially) equivalent to estimating a Gaussian ATSM. This is why we can construct a concentrated version of the intermediate quantity from the output of linear regressions. For USV models, we separate the parameters into three groups: θc = µg , Φg , Ω are the parameters that can be concentrated out, θm,h = (νh , Φh , Σh ) are the parameters that govern the dynamics of the volatility, and θm,b are the parameters that enter the bond loadings. We only need to optimize numerically over θm,h and θm,b , as the parameters in θc can be determined analytically as a function of θm,b . Moreover, the intermediate quantity can be additively separated into two pieces Q θ|θ (i) = Q1 θm,b , θc |θ (i)

+ Q2 θm,h |θ (i) . The first component corresponds to the first three terms in (21), and depends only on the parameters θm,b and θc . The remaining terms in (21) are the second component Q2 θm,h |θ (i) , which depend only on the volatility parameters θm,h . Our procedure can be implemented as follows: Procedure 2. The maximum likelihood estimator for USV models can be obtained by iterating over the following two steps: (a.) E-step: compute the expectations in the intermediate quantity Q θ |θ (i) from (21). (b.) M-step: maximize Q θ |θ (i) over θ to determine θ (i+1) . This can be separated into two sub-steps. (b1.) Maximize Q1 θm,b , θc |θ (i) with respect to θm,b and θc to (i+1)

(i+1)

determine θm,b and θc

. This can be solved equivalently

by maximizing Q1 θm,b , θˆc (θm,b )|θ (i) with respect to θm,b , where the concentrated objective function can be constructed as follows:

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

(i.) Given θm,b , calculate the bond loadings A and B and (1) 1 the state variables gt from xt = B− Yt − A1 . 1 (ii.) Given gt , run a GLS regression 1

gt +1 = µg + Φg gt + S¯t2 εg ,t +1

(22)

ˆ g θm,b . to calculate µ ˆ g θm,b , Φ ˆ θm,b as in (20). (iii.) Calculate the covariance matrix Ω ˆ g θm,b , (iv.) Substitute θˆc θm,b = µ ˆ g θm,b , Φ ˆ θm,b back into the intermediate quantity. Ω (b2.) Maximize Q2 θm,h |θ (i) with respect to θm,h to determine θm(i,+h1) . Appendix H contains details of our implementation. This procedure for USV models has many similarities with our earlier procedure for spanned models.4 The difference between the GLS regression in (19) of Procedure 1 and the regression in (22) of Procedure 2 is their covariance matrices. The EM algorithm imputes the latent values of ht by taking their expectations. Finally, we note that with only a few minor modifications, the analytical gradients for the likelihood of the spanned model from Proposition 3 can be used to calculate the gradients of the intermediate quantity in (21). In our experience, it takes only a few iterations of the EM algorithm to approach the maximum when estimating the USV models of Section 4.2.2. The rate of convergence of the EM algorithm near the maximum is known to be slow. Once the EM algorithm approaches the maximum, a researcher can switch to alternative estimation procedures for non-Gaussian state space models.5 For USV models, the expectation in Step 1 cannot be calculated in closed-form, requiring a Monte Carlo version of the EM algorithm. We calculate the expectations using sequential Monte Carlo methods or particle filters; see Creal (2012) for a survey. In particular, we use the particle filtering algorithm of Godsill et al. (2004) to draw paths of h0:T −1 from the joint posterior distribution p (h0:T −1 |Y1:T ; θ ); see Appendix H.2 for details. The particle filter also allows us to calculate filtered (one-sided) estimates of the volatility as well as an estimate of the likelihood function p (Y1:T ; θ). 3.3. Discussion

65

Example #2: Hidden factors Recently, Duffee (2011) argued that more than three factors are needed to explain the time-series dynamics of yields and risk premia, where these additional factors are ‘‘hidden’’ from the cross-section of yields because the factors are not priced. For simplicity, we illustrate the basic ideas here for Gaussian models. Extensions to spanned or unspanned non-Gaussian models are straightforward. The Gaussian state vector can be separated into sub-vectors ′ gt = g1′ ,t , g2′ ,t whose dimensions are G1 × 1 and G2 × 1, respectively. The dynamics under the P measure are g1,t +1 = µg ,1 + Φg ,11 g1t + Φg ,12 g2t + ε1,t +1

ε1,t +1 ∼ N 0, Σ0,g Σ0′ ,g g2,t +1 = µg ,2 + Φg ,21 g1t + Φg ,22 g2t + ε2,t +1 ε2,t +1 ∼ N 0, IG2 .

(23)

(24)

The dynamics of g1,t are the same under the Q measure but with Q the restrictions that Φg ,12 = 0 and the last G2 entries of δ1g are zero. These restrictions imply that only g1,t directly impacts yields as the bond loadings on g2,t are zero by construction. Given the parameters that enter the bond loadings, the factors (1) 1 that price bonds are conditionally observable g1,t = B− Yt − 1

A1 just as in step (i.) of Procedure 1. We can treat g1,t as the observed data and (23) is the new observation equation for a linear, Gaussian state space model. The remaining state variables g2t have transition equation (24) and are serially correlated shocks to the factors g1t that price bonds. We can use the Kalman filter to estimate this model, which is equivalent to a GLS regression where the errors are serially correlated. To concentrate the parameters (µg ,1 , µg ,2 , Φg ,11 , Φg ,21 ) out of the likelihood as in step (ii) of Procedure 1, we can either place these parameters in the state vector or use the augmented Kalman filter of de Jong (1991), see also Chapter 5 of Durbin and Koopman (2012). Example #3: parameter constraints In our approach, a researcher can impose constraints – such as in Ang and Piazzesi (2003), Kim and Wright (2005), Cochrane and Piazzesi (2008) and Bauer (2014) – and still concentrate out parameters by linear regression. We denote the penalized or constrained log-likelihood function ℓp (θ ) as

We discuss how our approach can be applied to a wide range of ATSMs.

ℓp (θ ) = log p (Y1:T ; θ ) + p(θ ),

Example #1: observable macroeconomic variables

where p(θ ) is the penalty term. If the constraints are only on the Q parameters, a researcher can directly apply our Procedures 1 and 2. If the goal is to constrain either the P parameters or the relationship between the P and Q parameters, the penalty term is a vector of Lagrange multipliers times the constraints. Step (ii.) of Procedure 1 or 2 can be replaced by constrained GLS. If the constraints are linear in µg , Φg , Φgh , there is a unique solution. Popular restrictions in the literature are all in this category. If the penalty term p(θ ) is a quadratic function of µg , Φg , Φgh , step (ii.) of Procedure 1 or 2 reduces to ridge regression and the parameters can be shrunk to a pre-specified value similar to a Bayesian VAR. A researcher may want to shrink the P parameters toward the Q parameters, which are often measured more precisely.

Our approach can estimate models with observable macroeconomic variables and with homoskedastic or heteroskedastic shocks; examples with homoskedasticity are Ang and Piazzesi (2003) and Hamilton and Wu (2012). Our procedure works the same as before except the state vector xt now contains the yield factors as well as the observed macroeconomic factors. For step (i.) of Procedure 1 or 2, we back out the latent component of xt condi(1) tional on a subset of the parameters, the yields Yt , and the macro variables. Given xt , we can concentrate a large number of parameters out of the objective function including many of the parameters introduced by adding the macroeconomic variables.

Example #4: serially correlated measurement errors 4 There are several versions of the EM algorithm all of which lead to the MLE; see Meng and Rubin (1993). These authors discuss issues such as concentrating the intermediate quantity as well as sequentially maximizing the intermediate quantity over subsets of θ . 5 USV models are an example of a non-linear, non-Gaussian state space model. General approaches for estimating non-Gaussian state space models include importance sampling (Durbin and Koopman (2012) and Richard and Zhang (2007)), particle filters (see Malik and Pitt (2011)), and MCMC (see Jacquier et al. (2007)).

Although often assumed to be i.i.d. normal in the literature, several authors have found the measurement errors ηt in (17) to be serially correlated; see, e.g. Hamilton and Wu (2014). Procedures 1 and 2 can be extended to cover models with serially correlated errors without increasing the dimension of the numerical optimization. For example, if ηt has VAR(1) dynamics ηt = Φη ηt −1 +κt with κt ∼ N (0, Ω ), then we can concentrate out θc = µg , Φg , Φgh , Ω , Φη

66

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

by least squares. Conditional on θm , we calculate the bond loadings A and B and factors xt and, from these, the values of ηt and ηt −1 which are implicitly a function of θm . To concentrate Φη and Ω out of the objective function, we compute their least squares estimates from the regression ηt = Φη ηt −1 + κt and plug these values back into the log-likelihood. Relation to Joslin et al. (2011) and Hamilton and Wu (2012) Both Joslin et al. (2011) and Hamilton and Wu (2012) use linear regression to estimate some parameters of Gaussian models.6 In the special case of Gaussian models with observable factors (A1 = 0 and B1 = I), our method is identical to the ML estimator of Joslin et al. (2011). Like the procedure in this paper, the approach of Hamilton and Wu (2012) works for a wider range of Gaussian models. Their minimum chi-square estimator is asymptotically equivalent to the ML estimator in this paper. The critical difference is that our procedures are designed for non-Gaussian models. Leveraging the analytical solution for linear regressions has not been explored in this area. For spanned models with both Gaussian and non-Gaussian factors, being able to rotate the factors is important. If a researcher takes an arbitrary basis of yields (such as principal components) and assumes that they can be separated a priori into observable Gaussian and non-Gaussian factors, it will restrict the likelihood. Our approach lets the data decide what linear combination of yields are the factors. 4. Data and parameter restrictions 4.1. Data

4.2.2. USV restrictions We focus on discrete-time USV models similar to the continuous-time models presented in Collin-Dufresne et al. (2009) and Joslin (2010).9 We label these models U1 (4) because they have one unspanned volatility factor and three Gaussian factors. USV restrictions are not unique. We present several models whose restrictions under Q result in non-Gaussian loadings where bn,h = 0 for all maturities. The first model, labeled U1 (4)(φ, φ 2 , ψ), has the following set Q of restrictions: (1) δ1,h = 0 and Σgh = 0. (2) Φg is a diagonal matrix with eigenvalues φ, φ 2 , ψ . (3) All entries of Σ1,g Σ1′ ,g are zero but the (1, 1) element. This entry is Σ12,g ,11 . (4) Φgh,3 = 0; Φgh,1 = Q

δ1,g ,1 Q Σ2 ; Φgh ,2 (1−φ) 1,g ,11

=−

(

)

δ12,g ,1 2(1−φ)2 δ1,g ,2 1−φ 2

Σ12,g ,11 . In this model, only the

Gaussian factor associated with the eigenvalue φ has stochastic volatility as the remaining entries of Σ1,g Σ1′ ,g are zero for all the other Gaussian factors. The USV restrictions force two of the eigenQ values of Φg to be related as φ and φ 2 . These restrictions summarize three different USV models depending on the relative size of φ and ψ . We label the models U1 (4)(φ > φ 2 > ψ), U1 (4)(φ > ψ > φ 2 ), U1 (4)(ψ > φ > φ 2 ). Each one of them is identified after imposing an ordering on the eigenvalues.10 A second model, labeled U1 (4)(φ, φ 2 , φ 4 ), allows two out of the three Gaussian factors to share a common stochastic volatility factor. The restrictions for this model are (1) δ1,h = 0 and Σgh = 0. Q (2) Φg is a diagonal matrix with eigenvalues φ, φ 2 , φ 4 . (3) All entries of Σ1,g Σ1′ ,g are zero but the (1, 1) and (2, 2) elements. These entries are Σ12,g ,11 and Σ12,g ,22 . (4) Φgh,1 = Q

We use the Fama and Bliss (1987) zero coupon bond data available from the Center for Research in Securities Prices (CRSP). The data is monthly and spans from June 1952 through June 2012 for a total of T = 721 observations with maturities of (1, 3, 12, 24, 36, 48, 60) months. For three factor models, the yields (1) measured without error Yt include the (1, 12, 60) month matu(1) rities. In models with four spanned factors, Yt are the (1, 12, 24, 60) month maturities. 4.2. Parameter restrictions 4.2.1. Identifying restrictions for spanned models We impose the following restrictions for identification. For the Q Q Gaussian part, these are: (i) µg = 0; (ii) Φg in ordered Jordan 7 form; (iii) δ1g = ι is a column vector of ones; and (iv) Σi,g is lower Q triangular. For the non-Gaussian part, (i) µh = 0. (ii) Φgh = 0.

(iii) Elements of the vector δ1h = ±1 can take either sign,8 which unlike Gaussian-only models will lead to inequivalent maxima as we explain in Section 5.2; (iv) Σh is diagonal. To guarantee non-negativity and admissibility of the factors, we also impose the discrete-time equivalent of the Feller condition νh,i > 1 and νhQ,i > 1 for i = 1, . . . , H. The matrices Σh , Σh−1 Φh Σh

and Σh−1 Φh Σh must also be non-negative.

Q

δ1,g ,2

(1−φ 2 )

(1−φ 2 )δ2 Q Σ12,g ,22 − 2(1−φ)2 δ1,g ,1 Σ12,g ,11 ; Φgh,3 = 1,g ,2

δ1,g ,1 2 Q Σ1,g ,11 ; Φgh,2 = 1−φ (1−φ 4 )δ12,g ,2 2 Σ1,g ,22 . − 2 2 1−φ 2 δ

(

)

1,g ,3

In this model, two diagonal entries in the covariance matrix Σ1,g Σ1′ ,g are non-zero. The additional flexibility in Σ1,g Σ1′ ,g comes at the expense of another restriction on the diagonal components Q of Φg . This model is unique because φ > φ 2 > φ 4 . These restrictions lead to the following proposition. Proposition 4. If the model has risk neutral dynamics (1)–(4) and satisfies the restrictions of U1 (4)(φ, φ 2 , ψ) or U1 (4)(φ, φ 2 , φ 4 ), then the model exhibits USV where bn,h = 0 ∀ n. Proof. See Appendix G. Intuitively, USV models free up the volatility factors to fit the heteroskedasticity of yields because they still enter the covariance matrix of the Gaussian factors in their P dynamics. The cost of Q adding USV factors comes from the constraints they place on Φg . These constraints can sacrifice the cross-sectional fit of the model. We impose the identifying restrictions of Section 4.2.1 for the Gaussian factors. For the non-Gaussian portion of the model, we need an additional restriction on the scale parameters Σ1,g or Σh , as they enter the likelihood in the same way. We set Σh = 0.01.

Q

6 The work by Adrian et al. (2012) and Diez de Los Rios (2013) are also similar in spirit to these methods. They focus only on Gaussian models. Our discussion here applies to these papers as well. 7 For the case where Φ Q has real distinct eigenvalues, it is a diagonal matrix with g

5. A three factor model In this section, we use a three factor model with one spanned volatility factor A1 (3) (in the Dai and Singleton (2000) notation) to demonstrate the performance of our method and to discuss the local maxima that arise in spanned non-Gaussian models. This

diagonal elements in descending order. 8 In theory, δ = 0 is also admissible and creates additional local maxima. We 1h

would like to thank an anonymous referee for pointing this out. Unlike a Gaussian model, Jensen’s inequality term prevents bn,h from being zero for maturities n ≥ 2 even when δ1h = 0. In practice, this model is not well identified, as the factor loadings coming from Jensen’s inequality term are extremely small and close to 0. We found that the likelihood values under this restriction are significantly smaller in such models.

9 Papers that document the potential existence of USV factors include Heidari and Wu (2003), Li and Zhao (2006), Trolle and Schwartz (2009), and Andersen and Benzoni (2010). 10 We do not consider the case where the eigenvalues may be identical and we restrict our attention to eigenvalues that are less than one.

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

model has been the preferred model by many researchers in the literature.11 For an A1 (3) model, the concentrated likelihood drops the number of parameters by one-third from 24 parameters to 16 parameters.

Table 1 Local maxima in the A1 (3) model.

ht

Φh

Q

5.1. Performance comparison To illustrate the mileage we gain from using our method, we compare our approach to the conventional method that does not concentrate out (µg , Φg , Φgh ) or use analytical gradients. We perform an experiment where we estimate the A1 (3) model on the CRSP dataset 100 times from 100 different starting values using both methods.12 We compare our method and the direct approach along two dimensions: convergence and speed. To measure the former, we use the likelihood ratio (two times the difference in loglikelihoods). The global solution found by our method has a log-likelihood of 36 647.69 (estimates and quasi-maximum likelihood standard errors can be found on the right hand side of Table 2). We achieve an identical value for all 17 random starting values whenever the parameters were initialized in this region or mode of the parameter space.13 Seventeen equals the number of times (one-sixth) that it started in this region. Conversely, the conventional method does not find this log-likelihood once nor does the method reproduce the same (incorrect) estimates for each of these 17 starting values. The highest log-likelihood value found by the standard approach is 36 645.29. The differences between the two methods is significant across these 17 starting values. The conventional method yields log-likelihood values ranging between 36 645.29 and 36 636.82, and it only achieves the former for one starting value. With our method producing the same number repeatedly, we can conclude that it is a maximum. An immediate benefit of the stable behavior of our method is that we are able to find that the A1 (3) model has 6 local modes with three well-behaved local maxima and three regions of the log-likelihood that appear to be locally unidentified. The three well-behaved local maxima are listed in Table 1 and we will discuss the properties of the model that create these local modes in Section 5.2. Our method converges 17/100 times to Local 1, 14/100 times to Local 2, and 17/100 times to Local 3. The median likelihood ratio between our procedure and the unconcentrated log-likelihood with no analytical gradient is 29.5 indicating a substantial difference between the two procedures. The conventional method, even if it gets close to a local maximum, always stops before it fully converges. This makes it difficult for researchers to differentiate between points that are near a wellbehaved local maximum that have the same economic meaning and locations corresponding to local maxima that are economically different. Estimation time is another important dimension along which we compare our approach to the conventional method. The median estimation time for our procedure to estimate from a random starting value is less than 2 min, whereas the conventional approach of directly maximizing the log-likelihood function takes more than 2 h. To perform our study with 100 starting values, it takes our method about 4 h, whereas it takes roughly 9 days to complete the same exercise with the conventional method.

11 This model has been widely considered as the benchmark non-Gaussian ATSM, see Dai and Singleton (2000), Cheridito et al. (2007), Collin-Dufresne et al. (2008), and Aït-Sahalia and Kimmel (2010) for examples. 12 To make the comparison as parallel as we can, we write the likelihood function the same way, impose the same identifying restrictions, and use the same scaling and initial values for the parameters except that the conventional method has additional parameters entering the numerical optimizer. 13 We consider two log-likelihoods to be numerically identical if they agree up to 2 decimal points. In practice, the log likelihood values are identical up to 8 decimal points.

67

Φg

Q

δ1,h LLF

Local 1

Local 2

Local 3

Level 0.9961

Slope 0.9528

Curvature 0.5812

0.9514 0.5358 1

1.0001 0.5881 1

0.9983 0.9389 −1

36 647.69

36 482.56

36 530.19

Estimates of Φh and Φg from the A1 (3) model with corresponding log-likelihood. Each value is a different local maximum depending on the sign of δ1,h and whether the non-Gaussian factor is the level, slope, or curvature. Q

Q

In summary, our method addresses all of the following problems with the conventional method. The conventional method is painfully slow. It does not achieve the global maximum. And, it is extremely hard to assess convergence behavior and the number of local maxima because conventional approaches do not repeatedly find the same local maximum even when started in that region of the parameter space. 5.2. Local maxima Using our approach simplifies estimation and helps uncover some features of the log-likelihood surface that may be obscured by directly maximizing the log-likelihood. In this section, we discuss the characteristics of the model that create local maxima and their economic consequences. In Gaussian models, a change in the sign of δ1g rotates the factors from gt to −gt . This rotation is economically irrelevant because the estimated model switches between two global maximums. Unlike Gaussian models, fixing the sign of δ1h in spanned models is not inconsequential, because the state variable ht is positive by definition. Changing the sign of δ1h does not rotate ht to −ht . Therefore, there can exist inequivalent local maxima for each combination of different signs of δ1h . For each of the local maxima, the estimated latent factor ht is different, which changes the conditional variance of gt and consequently the log-likelihood. Reordering the eigenvalues in Φ Q has completely different implications for spanned non-Gaussian models than for Gaussian models.14 If the eigenvalues are reordered in a multi-factor Gaussian model, it implies equivalent global maxima with the same economic implication. However, with non-Gaussian spanned factors, they can yield inequivalent local maxima. Using the A1 (3) model as an example, the factors are labeled as level, slope and Q Q curvature. Reordering the eigenvalues across Φg and Φh does not generally change the shape of the factors but it does change whether ht is the level, slope, or curvature. Any change in ht from one type of factor (level) to another (slope/curvature) implies a different conditional variance for gt making the likelihood no Q longer equivalent. Changing the order of the eigenvalues within Φg Q and/or Φh results in an equivalent global maximum. The intuition is the same as re-ordering of the factors gt within a Gaussian ATSM. In an ATSM with spanned non-Gaussian factors, it is not clear a priori which local maximum created by these characteristics of the model will be the global maximum. One must intentionally search each region that potentially has a local maximum and compare their likelihood values. To illustrate this idea, we present different local maxima for the A1 (3) model corresponding to different signs of δ1,h and different orderings of the eigenvalues. We report Φ Q , δ1,h , and log-likelihood values in Table 1. In the first column, ht is

14 We collect the autoregressive parameters together in matrices as

Φ=

Φh Φgh

0

Φg

ΦQ =

Φh Q Φgh Q

0

Φ

Q g

.

68

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

Table 2 Maximum likelihood estimates for the A0 (3) and A1 (3) models. G = 3, H = 0

LLF =

37 080.94

µg 6.97e−05 (4.92e−05)

−4.85e−05 (1.12e−04)

−3.37e−04 (7.06e−05)

1.0074 (0.0084)

0.0475 (0.0137)

0.0665 (0.0316)

−0.0115

0.9375 (0.0309) −0.0585 (0.0183)

0.0192 (0.0607) 0.6306 (0.0535)

0 –

0 –

Φg

0.0077 (0.0132) −0.0405 (0.0077)

0.9854 (0.0224) −0.0729 (0.0152)

δ0

Σh νh

µQg

4.09e−05 –

0 –

Φh

Φg

0.9961 (0.0006)

0.9514 (0.0029)

Q

Q

0.9538 (0.0031)

0.5299 (0.0295)

3.99e−04 (2.62e−05)

0 –

0 –

−3.09e−04 (4.54e−05) −4.50e−06 (2.23e−05)

5.09e−04 (3.66e−05) −2.52e−04 (2.82e−05)

0 – 3.78e−04 (2.37e−05)

1 –

1 –

Σ0,g

νh 1.9332 (2.0989)

0.0657 (0.0462) 0.6434 (0.0426) 0 –

δ0 −0.0011 (0.0002)

νhQ 2.6371 (0.2879)

Q

0.5358 (0.0293)

Σh 1.55e−05 (1.68e−06)

Σgh −0.8932

δ1,g √

diag (Ω ) × 1200 3m 0.2243

Φgh

0.0083 (0.0005)

Q

1 –

3.32e−05 (2.78e−05)

0.9943 (0.0061)

Φg

0.9950 (0.0007)

36 647.69

µg −1.34e−05 (5.14e−05)

Φh

µQg 0 –

LLF =

Σh νh 3.00e−05 –

Φg

(0.0184) −0.0366 (0.0111)

G = 2, H = 1

Σ0,g

Σ1,g

(0.0587) 0.0538 (0.0397)

3.61e−11 (1.11e−11) −1.89e−11 (8.52e−12)

δ1,h

δ1,g

1 –

1 –

1 –

2 yr 0.1248

3 yr 0.1236

0 –

−5.80e−12 (1.30e−12)

0.0063 (0.0004) −0.0035 (0.0003)

0 – 0.0046 (0.0003)

√

2 yr 0.1251

3 yr 0.1235

4 yr 0.1077

diag (Ω ) × 1200 3m 0.2245

4 yr 0.1078

Maximum likelihood estimates with quasi-maximum likelihood standard errors. Left: Gaussian A0 (3) model. Right: non-Gaussian A1 (3) model. The identifying restrictions

µQg = 0, δ1,g = ι, and δ1,h = 1 are imposed during estimation.

the level factor and δ1,h is positive. This is the global maximum in this case. In our sample, volatility is high during episodes where interest rates are high, so the level factor tends to explain the volatility best and δ1,h is positive. The next two columns present what happens when ht is the slope or curvature factor. Due to the nature of the data we are using, the likelihood function drops significantly from the global maximum to these alternative local maxima. In theory, there are six potentially different local maxima for each combination of eigenvalues and sign of δ1,h but in practice there are only three well-behaved local maxima. For the rest, we observe parameters hitting a boundary or eigenvalues of the autoregressive parameters being numerically 1. Each of these local maxima correspond to models where volatility declines in the 1970s, which contradicts the data. Those points can be locally unidentified, meaning that there exists a region of the parameter space where a subset of the parameters are unidentified. For models with unspanned volatility factors, the issues of multiple modes do not appear, because the likelihood surface of USV models is similar to spanned Gaussian models. 6. Model comparison In this section, we use our methodology to estimate a collection of popular and prominent ATSMs. The models that we estimate include three factor spanned models AH (3) with H = 0, 1, 2, 3 volatility factors and four factor spanned models AH (4) with H = 0, 1. We also estimate four USV models from Section 4.2.2: three versions of U1 (4)(φ, φ 2 , ψ) as well as the U1 (4)(φ, φ 2 , φ 4 ) model. We report the estimates as well as quasi-maximum likelihood standard errors as in White (1982) for eight models that have better empirical performance, leaving out the estimates for models with lower likelihoods for brevity.

6.1. Estimates and model fit Three factor spanned models Estimates with standard errors for the A0 (3) and A1 (3) models are included in Table 2. These two models have historically drawn most of the attention in the term structure literature. For the A0 (3) model, the concentrated likelihood drops the number of parameters entering the numerical optimizer from 22 to 10 (excluding Ω ), while this number drops from 24 to 16 for the A1 (3) model. The likelihood for the A0 (3) model at 37 080.94 is significantly higher than the A1 (3) model at 36 647.69. Surprisingly, among the models that we estimated, the A3 (3) model has the highest likelihood with a value of 37 385.28, substantially higher than the A0 (3) model. We report estimates for this model on the right panel of Table 3. In order to satisfy Q the admissibility restrictions, all values in Σh , Σh−1 Φh Σh and −1 15 Σh Φh Σh must be non-negative. We also impose the discreteQ time equivalent of the Feller condition, i.e. νh,i > 1 and νh,i > 1

for i = 1, . . . , H.16 In the A3 (3) model, we found that some of these parameters were near their boundaries. We fix them at the boundary when calculating the standard errors. For comparison purposes, we also report on the left panel of Table 3 estimates of Q an A3 (3) model where the matrices Φh and Φh are restricted to be

15 Multiple non-Gaussian factors tend to make the inversion problem more complicated. We solve this with randomized starting values, and use several iterations of ‘‘fminsearch’’ in MATLAB for each starting value to get an initial value that does not produce negative volatility. It takes roughly 20 trials (within a minute) to locate a starting value for the A3 (3) model. 16 No arbitrage requires that the stochastic discount factor M is strictly positive t +1

which is guaranteed by the Feller condition, see Appendix D for details.

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

69

Table 3 Maximum likelihood estimates for two A3 (3) models. Diagonal

LLF =

36 949.94

Flexible

Σh νh

LLF =

37 385.28

Σh νh

5.48e−05 –

2.75e−04 –

3.27e−04 –

5.35e−07 –

4.45e−05 –

9.35e−06 –

7.2325 (5.2324)

10.7472 (4.9148)

5.3456 (1.2905)

1 –

3.1176 (2.7588)

1 –

0 – 0.9511 (0.0136) 0 –

0 – 0 – 0.8133 (0.0298)

0.9587 (0.0102) 0.3109 (0.2807) 0

0.0010 (0.0003) 0.9133 (0.0259) 0.0226 (0.0018)

0 – 0.5975 (0.1841) 0.7876 (0.0408)

3.10e−04 –

6.63e−04 –

8.10e−07 –

1.43e−05 –

9.35e−06 –

12.0937 (5.2723)

10.8396 (1.9412)

1 –

1 –

0.0002 (0.0000) 0.8850 (0.0333) 0.0206 (0.0047)

0 – 1.3287 (0.2776) 0.6122 (0.0490)

0 – 1.43e−05 (3.05e−06) 0 –

0 – 0 – 9.35e−06 (2.43e−06)

1 –

−1

2 yr 0.1274

3 yr 0.1236

νh

νh

Φh

Φh

0.9918 (0.0056) 0 – 0 –

Σh νh

Σh νh

Q

2.01e−05 –

Q

νhQ 2.6464 (0.7721)

δ0

νhQ

-0.0070 (0.0016)

1.5151 (1.7822)

Φh

Φh

Q

0.9996 (0.0005) 0 – 0 –

Q

0 – 0.9479 (0.0018) 0 –

0 – 0 – 0.4837 (0.0239)

0.9890 (0.0014) 0 – 0.4968 (0.0964)

0 – 2.56e−05 (6.60e−06) 0 –

0 – 0 – 6.12e−05 (5.95e−06)

5.35e−07 (5.84e−08) 0 – 0 –

1 –

−1

1 –

Σh 7.58e−06 (1.27e−06) 0 – 0 –

Σh

δ1,h 1 –

δ1,h –

√

diag (Ω ) × 1200 3m 0.2247

δ0 −9.12e−04 (4.82e−04)

–

√ 2 yr 0.1324

3 yr 0.1286

4 yr 0.1112

diag (Ω ) × 1200 3m 0.2230

4 yr 0.1077

Maximum likelihood estimates with asymptotic quasi-maximum likelihood standard errors. Left: A3 (3) model with diagonal matrices Φh and Φh . Right: A3 (3) model. The identifying restrictions µh = 0, and δ1,h = (1, 1, −1) are imposed during estimation. Q

diagonal in which case the estimated parameters are not close to the boundaries. Finally, we note that the A2 (3) model (estimates not reported for brevity) had the lowest likelihood among the models we estimated with a value of 36 120.34. The bond loadings for the A0 (3), A1 (3), and diagonal A3 (3) models are almost the same, see Fig. 1. This happens for two reasons. First, the functional form of the bond loading recursions are the same up to Jensen’s inequality. Second, the size of the measurement errors in the cross-section of yields are small relative to the magnitude of the time series shocks causing an efficient estimator (like maximum likelihood) to emphasize the fit of the cross section. The estimated factors of these models have high correlations (ranging from 0.95 to 1). Interestingly, the bond loadings for the unrestricted A3 (3) model appear to be nonQ stationary, even though the eigenvalues of Φh for this model are inside the unit circle and are nearly identical to the other models.17 The average pricing errors are similar in these models, with the more restricted diagonal A3 (3) model having larger values. The measurement errors are about 22 basis points for 3 months, 12 for

17 Unlike Gaussian models, the bond loading recursions b¯ for multi-factor nonn, h Gaussian models in (6) are non-linear difference equations. The stability conditions Q for these equations appear to be determined by more than the eigenvalues of Φh . ¯ We leave the conditions for stability of bn,h for future research.

2 years, 12 for 3 years and 11 for 4 years. These errors have the same magnitude as those reported in Ang and Piazzesi (2003). Differences between these models are largely driven by their ability to fit the time series component of the likelihood. For example, in the A1 (3) model, the conditional mean under P is more restricted than the A0 (3) model as the Gaussian factors cannot enter the conditional mean of the non-Gaussian factors. The economic implication of this restriction for the A1 (3) model is that the level factor does not depend on the past values of slope and curvature factors. This is apparently counterintuitive. If the slope is high in the last period, i.e., the long rate is much higher than the short rate (more than explained by compensating for risk), then it means the market expects the short rate will increase in the future. On average, the next periods’ short rate or level will increase. Four factor USV models The two best fitting USV models are the U1 (4)(φ > φ 2 > ψ) and U1 (4)(ψ > φ > φ 2 ) models, whose estimates are in Table 4. See Section 4.2.2 for the definition of the models. We calculate the likelihood for USV models by the particle filter, see Creal (2012). The likelihood for the best fitting U1 (4)(φ > φ 2 > ψ) model is 37 331.25, which is lower than the unrestricted A3 (3) model but substantially higher than other spanned models. The bond loadings for USV models are not as flexible compared to spanned models due to the USV restrictions explained in

70

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

Fig. 1. Bond loadings for six different three factor models. Bond loadings B as a function of maturity n. Top row from left to right: A0 (3) model, A1 (3) model, A3 (3) model. Q Second row: A3 (3) model with diagonal Φh , U1 (4)(φ > φ 2 > ψ) model, U1 (4)(ψ > φ > φ 2 ) model. Reported in each graph are the eigenvalues of the risk-neutral feedback matrix Φ Q . To make the restricted A3 (3) model easily comparable, we report the absolute value of the loadings for this model.

Section 4.2.2. The first USV model constrains the two largest eigenvalues of Φ Q for the level and slope factors to be related as φ and φ 2 . Without any restrictions on the eigenvalues of Φ Q as in the A0 (3) and A1 (3) models, these are estimated to be 0.995 and 0.954. In the U1 (4)(φ > φ 2 > ψ) model, if the level factor has φ = 0.995, the USV restriction would require the slope factor to be more persistent φ 2 = 0.990 than it would be without the restriction (0.954). Conversely, if the model tried to fit the slope factor first, the USV restriction would not allow the level factor to be persistent enough. As a compromise, the two eigenvalues in the U1 (4)(φ > φ 2 > ψ) model are closer to each other than what the data would like with estimated values φ = 0.9868 and φ 2 = 0.9738. Consequently, the loadings on the level and slope factors in this model lie in between the loadings for spanned models, see Fig. 1. Due to the more restrictive loadings, it is not surprising that the average pricing errors for this model are larger than for spanned models. Moreover, the Gaussian factors whose eigenvalues share the relationship φ and φ 2 are highly correlated with a correlation of −0.928. This number is 0.43 in absolute value for the A0 (3) and A1 (3) models for example. Estimation of an A0 (3) that includes the same restrictions on Q Φg as the U1 (4)(φ > φ 2 > ψ) model but no stochastic volatility has a likelihood of 37 040.5. This indicates that this single USV restriction is rejected relative to the benchmark Gaussian A0 (3) model of Table 2 by a likelihood ratio test. On the other hand, the addition of unspanned stochastic volatility factors increases the likelihood by 37 331.25 − 37 040.5 = 290.75. The primary source is a significantly better fit of the dynamics of volatility, as we discuss further below. When the USV restriction is imposed on the second and third largest eigenvalues as in the U1 (4)(ψ > φ > φ 2 ) model of Table 4, it constrains the bond loadings of the slope and curvature factors Q because the eigenvalues in Φg associated with these factors now share the relationship φ and φ 2 . This restriction causes them to be closer than they would otherwise have been if left unrestricted. The U1 (4)(φ, φ 2 , φ 4 ) model imposes even stronger restrictions on

Φg because it only has a single free parameter. It has likelihood Q

36 889.44 (parameter estimates not reported), which is well below the benchmark A0 (3) model. Four factor spanned models Finally, we consider two four factor spanned models: the Gaussian A0 (4) model and the non-Gaussian A1 (4) model. There are a total of 35 parameters in the A0 (4) model and only 20 of these parameters enter the numerical optimizer, while there are 39 parameters in the A1 (4) model and 15 of these can be concentrated out. Parameter estimates and standard errors for both models are in Table 5. While adding another Gaussian factor increases the likelihood relative to the A0 (3) model to 37 194.34 for the A0 (4), the additional factor is not as important as adding an USV factor. Repeated eigenvalues When we estimated the A1 (4) model, we found that it had repeated eigenvalues and our estimates in Table 5 have imposed them using the Jordan decomposition. If we use a diagonal matrix for Φ Q , the matrix B1 will be singular and one element in Φ Q is unidentified.18 To illustrate this point, Table 6 reports different sets of parameter values all with the same likelihood. Across the four Q Q local maxima, the values of Φh , the last eigenvalue of Φg , and the log-likelihood function are almost identical but the first two Q eigenvalues of Φg vary across different optima. The last column of Table 6 shows the results when we impose repeated eigenvalues Q (from the model reported in Table 2). The estimates of Φh , the Q last eigenvalue of Φg and the likelihood function have the same Q values as before. However, the first two eigenvalues of Φg are identical by definition and are equal to the average of the first two eigenvalues in those local maxima. The log-likelihood value also does not change.

18 For a matrix with repeated eigenvalues, the Jordan decomposition imposes the additional restrictions necessary to obtain identification. With two repeated eigen Q ′ Q values in Φg , the Jordan decomposition is vec Φg = (λ1 1 0 ; 0 λ1 0 ; 0 0 λ2 ) where λ1 and λ2 are the unique eigenvalues.

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

71

Table 4 Maximum likelihood estimates for two U1 (4)(φ, φ 2 , ψ) models. G = 3, H = 1

LLF =

G = 3, H = 1

37 331.25

µg 2.11e−04 (5.04e−07)

LLF =

37 241.05

−1.17e−03 (8.42e−05)

8.49e−06 (5.76e−06)

0.3251 (0.1009) 0.4463 (0.0715) 0.0349 (0.0233)

0.0520 (0.0161) −0.0923 (0.0149) 0.9995 (0.0029)

0 –

0 —

0.8231 –

0.9967 (0.0004)

0 – 2.13e−05 (4.40e−05) 7.05e−05 (1.43e−04)

0 – 0 – 3.31e−04 (9.73e−05)

µg −1.88e−04 (3.66e−07)

−2.37e−04 (2.39e−06)

7.42e−04 (8.80e−05)

0.1181 (0.0008) 0.8736 (0.0015) −0.0372 (0.0001)

0.1107 (0.0009) −0.0555 (0.0005) 0.6347 (0.0070)

1.1316 (0.0427) −0.2881 (0.0213) 0.0348 (0.0162)

0 –

0 –

0.9738 –

0.4931 (0.0019)

0.9073 (0.0045)

0 – 2.44e−04 (1.14e−07) −4.27e−04 (2.10e−06)

0 – 0 – 1.66e−10 (3.90e−13)

1.21e−03 (7.24e−05) −1.30e−03 (4.63e−06) −3.67e−04 (1.06e−04)

0 – 0 – 0 —

0 – 0 – 0 –

8.35e−04 (4.03e−05) 0 – 0 –

0 – 0 – 0 –

0 – 0 – 0 –

νh

Φh

Σh

νh

Φh

Σh

1 –

0.9598 (0.0086)

0.0100 –

1.0658 (0.0319)

0.9539 (0.0101)

0.0100 –

1 –

1 –

1 –

1 –

1 –

2 yr 0.1298

3 yr 0.1316

Φg 1.0724 (0.0020) −0.0788 (0.0008) −0.0321 (0.0000)

Φg

µQg 0 –

δ0

µQg

0.0061 (0.0000)

0 –

Φg

Q

Σ0,g 1.16e−03 (2.41e−06) −1.23e−03 (2.83e−06) 5.04e−05 (9.23e−08)

Σ0,g

Σ1,g 7.05e−04 (1.44e−06) 0 – 0 –

Σ1,g

δ1,g 1 –

δ1,g √

√

diag (Ω ) × 1200 3m 0.2266

0.0115 (0.0007)

Φg

Q

0.9868 (0.0004)

δ0

2 yr 0.1298

3 yr 0.1316

4 yr 0.1101

diag (Ω ) × 1200 3m 0.2266

4 yr 0.1101

Maximum likelihood estimates with quasi-maximum likelihood standard errors for two unspanned models. Left: U1 (4)(φ > φ 2 > ψ) model. Right: U1 (4)(ψ > φ > φ 2 ) Q model. The USV and the identifying restrictions Σh = 0.01, µg = 0, and δ1,g = (1, 1, 1) are imposed during estimation.

6.2. Volatility The blue lines in Fig. 2 are the conditional volatilities from four different models; across the rows are the A1 (3), A3 (3), U1 (4)(φ > φ 2 > ψ) and U1 (4)(ψ > φ > ψ 2 ) models.19 The three columns represent maturities of 1 month, 1 year, and 5 years. The volatilities for the USV models are the filtered (one-sided) estimates calculated from the particle filter. To provide a point of comparison, we also plot in these graphs the unconditional volatility from the A0 (3) model in green and in red estimates of the conditional volatilities from the multivariate generalized autoregressive score model of Creal et al. (2011) and Creal et al. (2013).20 The estimated volatilities from spanned models (first two rows) are much less volatile than GAS volatility. A similar observation was made by Collin-Dufresne et al. (2009). With only a single volatility factor, the A1 (3) model does not fit yield volatility at any maturity. The A3 (3) model adds flexibility with two more factors,

19 The volatilities from the remaining models have the same qualitative features and are not shown. 20 The generalized autoregressive score model with time-varying covariance matrix is similar to a multivariate GARCH model. To make the volatilities of yields (1) comparable across models, we use a VAR(1) for the conditional mean of yields Yt and allow the errors to have time-varying volatilities and correlations.

which improves the fit for volatility at shorter horizons, especially for the recent episode of low volatility. At longer horizons, the fit to volatility appears to be the same as the A1 (3) model. This is because the estimated level (volatility) factor is similar across both the A1 (3) and A3 (3) models, and the long term volatility loads mostly on it. Ultimately, spanned models lack flexibility because the non-Gaussian state variables ht serve a dual role: they must simultaneously fit the conditional mean and variance. The maximum likelihood estimator chooses the parameter vector θ to fit the conditional mean first before fitting the conditional variance. USV models are designed to fit volatility by separating the role of Gaussian and non-Gaussian factors. The U1 (4)(φ > φ 2 > ψ) model performs much better than spanned models at fitting the volatility across different maturities, although it only has a single stochastic volatility factor. It does a particularly good job for short and medium maturities. At longer maturities, it fits the high volatility periods of the early 1980s well, although misses the period of low volatility early in the sample. Our results suggest, however, that not all USV models fit yield volatility equally well and researchers must be careful when choosing USV restrictions. We demonstrate this point by comparing the performance of the two best USV models U1 (4)(φ > φ 2 > ψ) and U1 (4)(ψ > φ > φ 2 ). Instead of having stochastic volatility on the level factor, the U1 (4)(ψ > φ > φ 2 ) model has stochastic

72

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

Table 5 Maximum likelihood estimates for the A0 (4) and A1 (4) models. G = 4, H = 0

LLF =

37 195.84

µg 3.89e−04 (1.73e−04)

G = 3, H = 1

LLF =

Σh νh

µg 3.07e−04 (1.76e−04)

−9.73e−04 (2.11e−04)

1.28e−03 (3.07e−04)

−8.19e−04 (3.68e−04)

5.64e−05 –

1.0336 (0.0352)

0.0874 (0.0536)

−0.0142

0.0908 (0.0646)

0.9901 (0.0073)

−0.0783

0.8338 (0.0893) 0.1541 (0.1312) −0.1405 (0.0901)

0.2044 (0.1298) 0.6935 (0.1731) −0.0398 (0.0924)

(0.1734) 0.1926 (0.2259) 0.5456 (0.1377)

0 –

0 –

0 –

0.9604 (0.0116)

0.8764 (0.0303)

0.6964 (0.0504)

6.94e−04 (2.33e−04)

0 –

0 –

0 –

2.11e−05 (3.75e−06)

Σgh

Σ0,g

−1.47e−03 (2.80e−04) 1.66e−03 (1.96e−04) −8.06e−04 (3.69e−04)

9.77e−04 (3.98e−04) −1.39e−03 (4.33e−04) 5.65e−04 (2.68e−04)

0 – 8.74e−04 (3.21e−04) −7.18e−04 (3.59e−04)

0 – 0 – 4.04e−04 (2.68e−05)

0.9248 (0.4737) −0.2171 (0.0333) −1.5635 (0.4198)

8.37e−04 (2.90e−04) −9.03e−05 (2.69e−05) −7.96e−04 (2.72e−04)

Φg

36 729.22

νh −3.76e−05 (1.84e−05)

−3.10e−04 (1.79e−04)

1.1695 (0.4487) 0.8156 (0.0456) −0.1235 (0.4361)

0.0911 (0.0974) 0.0028 (0.0097) 0.6508 (0.0974)

0 –

0 –

– –

0.7021 (0.0324)

0 – 7.43e−13 (4.22e−12) 9.41e−11 (1.75e−11)

0 – 0 – 4.03e−10 (6.48e−11)

0 – 8.15e−04 (7.96e−05) 1.12e−03 (4.10e−04)

0 – 0 – 4.68e−03 (3.04e−04)

1 –

2.6780 (1.8951)

Φh

(0.0626) 0.0889 (0.0766) −0.0851 (0.0528)

(0.0475)

−0.0788

µQg

Φgh

Φg

0.0037 (0.0260) 0.0016 (0.0031) −0.0352 (0.0257)

0.8668 (0.0644) 0.0153 (0.0062) −0.0095 (0.0555)

Σh νh

µQg

2.71E−05 –

0 –

Φh

Φg

0.9952 (0.0011)

0.9121 (0.0075)

Q

0 –

Φg

Q

Q

0.9922 (0.0024)

Σ0,g

νhQ 1.2839 (0.3516)

Q

Σh

Σ1,g 1.04e−02 (2.48e−03) −6.75e−04 (2.75e−04) −9.01e−03 (2.69e−03)

δ0 3.92e−03 (6.45e−04)

δ1,g 1 —

1 –

1 –

1 –

δ1,g

1 –

1 –

1 –

3 yr 0.1005

4 yr 0.0981

√

√

diag (Ω ) × 1200 3m 0.2390

δ0 −4.44e−04 (3.67e−04) δ1,h

3 yr 0.1013

diag (Ω ) × 1200 3m 0.2408

4 yr 0.0982

Maximum likelihood estimates with quasi-maximum likelihood standard errors. Left: Gaussian A0 (4) model. Right: non-Gaussian A1 (4) model. The identifying restrictions

µQg = 0, δ1,g = ι, and δ1,h = 1 are imposed during estimation. Table 6 Repeated eigenvalues.

Φ Φ

Q h Q g

LLF

tions needed to impose USV are not unique and the choice of which restriction to impose is not innocuous.

Local 1

Local 2

Local 3

Repeated

0.9952 0.9130 0.9112 0.7021

0.9952 0.9164 0.9075 0.7025

0.9952 0.9126 0.9115 0.7021

0.9952 0.9121 – 0.7021

36 729.22

36 729.20

36 729.22

36 729.22

Estimates from the A1 (4) model when repeated eigenvalues are not imposed compared to when they are. The table illustrates how this can create identification problems in affine models.

volatility on the slope factor. It fits the volatility at shorter maturities equally well but, at longer maturities, it exhibits almost no stochastic volatility. This is because the eigenvalue associated with the slope factor that has stochastic volatility is φ = 0.9073 as opposed to φ = 0.9868 in the first USV model. The bond loadings on this factor decay rapidly as maturity increases, meaning that long maturities have no stochastic volatility. To summarize, although USV models are designed to fit the volatility, the restric-

7. Conclusion We provide new estimation procedures for non-Gaussian affine term structure models with spanned or unspanned stochastic volatility. The new estimation approach for spanned models leverages the fact that many of the parameters can be concentrated out of the likelihood function. By optimizing the concentrated likelihood, it provides exactly the same solution as maximizing the original likelihood. But, it improves the estimation dramatically by reducing the number of parameters that need to be maximized numerically. We demonstrate the improvement in performance with our method. Using our procedure, we show what characteristics of spanned non-Gaussian models cause local maxima to exist and how alternative local maxima may have dramatically different economic implications. We apply a similar idea to the concentrated likelihood to estimate unspanned models. Estimating a wide range of popular models, we find that models with spanned volatility have similar cross sectional fit for yields.

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

73

Fig. 2. Estimated conditional volatility of yields from four different affine models. Estimated conditional volatility of yields from different affine models compared to the multivariate generalized autoregressive score model and the A0 (3) with constant volatility. Across the columns are the one month, one year, and 5 year maturities. Top row: A1 (3) model. Second row: A3 (3) model. Third row: U1 (4)(φ > φ 2 > ψ) model. Bottom row: U1 (4)(ψ > φ > φ 2 ) model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

They fit better than the USV models, because the latter impose restrictions on the cross section in order to introduce USV factors. Models with USV fit the volatility better by design. The choice of how to impose USV restrictions is not innocuous as some USV models can severely limit yield volatility at particular maturities.

USV models make an effort to fit the volatility of yields. Future work on term structure models aiming to fit both the conditional mean and volatility of yields simultaneously will likely require (1) multiple unspanned volatility factors, and (2) the ability to relax the restrictions that USV impose on the cross section.

74

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

Appendix A. Distributions

×

Using the definition of the modified Bessel function of the first kind,21 the p.d.f. can be expressed as

A.1. Gamma and multivariate gamma distributions A gamma r.v. wt +1 ∼ Gamma (νh , κ) has p.d.f. p (wt +1 |νh , κ) = w ν −1 w h κ −νh exp − tκ+1 and E (wt +1 ) = νh κ and V (wt +1 ) = Γ (ν ) t +1 1

h 1 νh νh κ 2 . The Laplace transform is E [exp (uwt +1 )] = 1−κ , which u exists only if κ u < 1. A multivariate gamma random vector ht +1 ∼ Mult. Gamma(νh , Σh , µh ) can be obtained by shifting and rotating a vector of uncorrelated gamma r.v.’s. Let ht +1 = µh + Σh wt +1 where wt +1 is an H × 1 vector with elements wi,t +1 ∼ Gamma νh,i , 1 for i = 1, . . . , H. The H × 1 vector of (non-negative) location parameters is µh , Σh is a full rank H × H matrix of (non-negative) scale parameters, and νh > 0 is a H × 1 vector of shape parameters. The

p.d.f. of ht +1 can be determined by a standard change-of-variables p (ht +1 |νh , Σh , µh ) H

ν −1 1 e′i Σh−1 [ht +1 − µh ] h,i Γ ν h ,i i =1 ′ −1 × exp −ei Σh [ht +1 − µh ]

where ei is an H × 1 unit vector. The mean and variance are E [ht +1 ] = µh + Σh νh and V [ht +1 ] = Σh diag (νh ) Σh′ . The Laplace transform is E exp u′ ht +1

=

∞

exp u′ ht +1 p (ht +1 |νh , Σh , µh ) dht +1

0

= exp u′ µh

∞

exp u′ Σh wt +1

0

H

1

ν −1

Γ νh,i

i=1

wi,ht ,+i 1

× exp (−wt +1 ) dwt +1 νh,i H 1 = exp u′ µh 1 − e′i Σh′ u i =1 H ′ ′ ′ = exp u µh − νh,i log 1 − ei Σh u .

p (ht +1 |νh , Φh ht , Σh , µh )

= |Σh | exp − −1

H

ei Σh [ht +1 − µh ] + ei Σh Φh [ht − µh ] ′

−1

′

−1

i=1

×

νh,2i −1

H

e′i Σh−1 [ht +1 − µh ]

− νh,2i −1

e′i Σh−1 Φh [ht − µh ]

i =1

× Iνh,i −1 2 e′i Σh−1 [ht +1 − µh ] e′i Σh−1 Φh [ht − µh ] . The Laplace transform can be derived from the law of iterated expectations E exp u′ ht +1

= Ez Eh|z exp u′ ht +1 νh,i +zi H ′ 1 = Ez exp u µh 1 − e′i Σh′ u i=1 νh,i zi H H ′ 1 1 Ez = exp u µh 1 − e′i Σh′ u 1 − e′i Σh′ u i=1 i=1 νh,i H ′ 1 = exp u µh 1 − e′i Σh′ u i=1 H e′i Σh−1 Φh [ht − µh ] e′i Σh′ u × exp 1 − e′i Σh′ u i =1 H e′i Σh′ u ′ −1 = exp u′ µh + ′ ′ ei Σh Φh (ht − µh ) 1 − e Σ u i h i=1 H ′ ′ − νh,i log 1 − ei Σh u

i=1

where e′i Σh′ u denotes the ith element of the H × 1 vector Σh′ u. The Laplace transform exists only if e′i Σh′ u < 1 for i = 1, . . . , H. A.3. Mixture of Gaussian and mult. NCG distributions

i=1

The Laplace transform exists only if ei Σh u < 1 for i = 1, . . . , H. ′

1

1

e′i Σh−1 [ht +1 − µh ]

Γ νh,i + zi,t +1 zi,t +1 ! i =1 zi,t +1 =0 z × e′i Σh−1 [ht +1 − µh ] e′i Σh−1 Φh [ht − µh ] i,t +1 .

We define the distributions found in the paper. The notation for them is local to the appendix.

= |Σh−1 |

∞ νh,i −1

H

′

A.2. Multivariate non-central gamma distributions A H × 1 non-central gamma (NCG) random vector ht +1 ∼ Mult.-N.C.G. (νh , Φh ht , Σh , µh ) is a Poisson mixture of multivariate gamma r.v.’s as in (12)–(14). The process ht remains positive and well-defined as long as µh ≥ 0, Σh−1 Φh Σh ≥ 0, νh > 1, and elements of Σh cannot be negative. A standard multivariate NCG random variable (the discrete-time CIR process) is obtained by setting µh = 0 and letting Σh be a diagonal matrix. Further properties of the univariate NCG process are described in Gouriéroux and Jasiak (2006). As long as Σh has full rank, the p.d.f. can be found by integrating out the Poisson r.v.’s

From standard results in statistics, the multivariate (G × 1) Gaussian r.v. gt +1 ∼ N gt +1 |µg , Σg Σg′ has Laplace transform E exp u′ gt +1 = exp µ′g u + 12 u′ Σg Σg′ u for any real (G × 1) vector u. Consider a (G + H ) × 1 vector xt +1 = (h′t +1 , gt′+1 )′ where ht +1 is an H × 1 vector having a multivariate NCG distribution p (ht +1 |νh , Φh ht , Σh , µh ) and gt +1 is a G × 1 vector of conditionally Gaussian r.v. gt +1 ∼ N µg + Σgh ht +1 , Σg Σg′ . Let u = (u′h , u′g )′ where uh and ug are H × 1 and G × 1 vectors, respectively. Using the law of iterated expectations, the Laplace transform is E exp u′ xt +1

= E exp u′g gt +1 exp u′h ht +1 ′ ′

= Eh Eg |h exp ug gt +1 exp uh ht +1 1 = Eh exp (µg + Σgh ht +1 )′ ug + u′g Σg Σg′ ug exp u′h ht +1 2

p (ht +1 |νh , Φh ht , Σh , µh )

= |Σh | exp − −1

H i =1

ei Σh [ht +1 − µh ] + ei Σh Φh [ht − µh ] ′

−1

′

−1

21 This is defined as I (x) = λ Stegun (1964).

x λ ∞ 2

1 z =0 Γ (λ+z +1)z !

2 z x 4

, see Abramowitz and

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

1 = exp u′g µg + u′g Σg Σg′ ug Eh exp u′g Σgh + u′h ht +1 2 1 = exp u′g µg + u′g Σg Σg′ ug + u′gh µh

+

H

+

i=1

1 − e′i Σh′ ugh

ei Σh Φh (ht − µh )

Appendix B. Bond pricing Bond pricescan be solved by induction. Guess that bond prices are Ptn = exp a¯ n + b¯ ′n,h ht + b¯ ′n,g gt for some coefficients a¯ n , b¯ n,h , and b¯ n,g . At maturity n = 1 when the payoff is Pt0+1 = 1, we find

Pt1 = Et

exp (−rt ) Pt0+1

exp (−rt ) Ptn+−11 = Et

Q

gt +1

Q

Q

Q

Q ∼ N µQg + ΦgQ gt + Φgh ht + Σgh × ht +1 − (IH − ΦhQ )µh + Σh νhQ + ΦhQ ht , Σg ,t Σg′ ,t .

1 ′ bn−1,g Σg ,t Σg′ ,t b¯ n−1,g

′ ˜ 0,g Σ ˜ 0′ ,g = Cgg Σ0,g Σ0′ ,g Cgg Σ −

+

˜ i ,g Σ ˜ i′,g = Σ

i=1

1 − ei Σh bn−1,gh

Mt +1

e′i Σh−1 Φh [ht − µh ]

exp(−rt )p(gt +1 |It , ht +1 , zt +1 ; θ , Q)p(ht +1 |It , zt +1 ; θ , Q)p(zt +1 |It ; θ , Q) p(gt +1 |It , ht +1 , zt +1 ; θ , P)p (ht +1 |It , zt +1 ; θ , P) p(zt +1 |It ; θ , P)

where the distributions are Gaussian, gamma, and Poisson. The exact SDF is

νhQ,i log 1 − e′i Σh′ b¯ n−1,gh

exp −rt − 12 εg ,t +1 εg ,t +1 −

i =1

Mt +1 =

= exp −δ0 + a¯ n−1 + µ ¯

Q′ g bn−1,g

+ µh bn−1,h + Φ Σgh b¯ n−1,g ′

Q′ h

′

1

2

−

νhQ,i log 1 − e′i Σh′ b¯ n−1,gh

i =1

−

H

e′i Σh′ b¯ n−1,gh

i =1

1 − ei Σh bn−1,gh ′

′¯

e′i Σh−1 Φh µh + b¯ ′n−1,g ΦgQ − δ1′ ,g gt Q

Q′

exp − ε

Q

ε

1 ′ 2 g ,t + 1 g ,t + 1

−

H

e′i Σh−1 Φh [ht − µh ] Q

i=1

H

ei Σh Φh [ht − µh ] ′

−1

i =1

zi,t +1 Q Γ νh,i + zi,t +1 e′i Σh−1 Φh [ht − µh ] Q × Γ νh,i + zi,t +1 e′i Σh−1 Φh [ht − µh ] i=1 where εg ,t +1 = Σg−,t1 gt +1 − µg − Φg gt − Φgh ht − Σgh εh,t +1 and wt +1 = Σh−1 [ht +1 − µh ]. No arbitrage requires that Mt +1 > 0 which implies that wt +1 > 0. The Feller conditions νh > 1 and νhQ > 1 guarantee that wt +1 > 0 under both probability measures. H

′ ¯ − νhQ′ Σh′ Σgh bn−1,g + b′n−1,g Σ0,g Σ0′ ,g b¯ n−1,g H

′ ′ −1 Cgg Σj,g Σj′,g Cgg ej Chh ei .

j =1

We define the stochastic discount factor as

H

−

H

Appendix D. Stochastic discount factor

Q

′¯

′

′ ′ −1 Cgg Σi,g Σi′,g Cgg ei Chh ch

i=1

=

e′i Σh′ b¯ n−1,gh

H

2

Q Q Q + µQg + ΦgQ gt + Φgh ht − Σgh (IH − Φh )µh + Σh νh ′ ′ ′ + ΦhQ ht b¯ n−1,g + Σgh b¯ n−1,g + bn−1,h µh H

only under a new

−1 cg + Cgh ([IH − Φh ] µh + Σh νh ) µ ˜ g = cg + Cgg µg − Cgg Φg Cgg −1 −1 − Cgh Φh − Cgg Φg Cgg Cgh + Cgg Φgh Chh ch −1 −1 ˜ Φgh = Cgh Φh − Cgg Φg Cgg Cgh + Cgg Φgh Chh −1 ˜ gh = Cgh + Cgg Σgh Chh Σ

This expectation has the same form as the Laplace transform provided in Appendix A. We find

is a member of

−1 ˜ h = Chh Φh Chh Φ −1 ˜ g = Cgg Φg Cgg Φ

Ptn = exp −δ0 − δ1′ ,g gt − δ1′ ,h ht + a¯ n−1 +

′ ′

′

parameters θ˜ . The proof is immediate by comparing the Laplace transform of these random variables before and after rotating them. The mapping between the parameters θ˜ and θ is

µ ˜ h = ch + Chh µh ˜ h = Chh Σh Σ

where the expectation is taken with respect to the distribution of ′ the random vector xt +1 = h′t +1 , gt′+1 ht +1 ∼ Mult-NCG νh,i , Φh ht , Σh , µh

The admissibility restrictions needed to keep ht positive are: (1.) Chg = 0; (2.) Chh is such that all elements Chh Σh are nonnegative; (3.) ch is such that ch + Chh µh is non-negative; (4.) Chh and Cgg are full rank. For some values of θ , these restrictions may allow ch and Chh to be negative. the same family of distributions as h′t , gt

× exp a¯ n−1 + b¯ ′n−1,h ht +1 + b¯ ′n−1,g gt +1 = exp −δ0 − δ1′ ,h ht − δ1′ ,g gt + a¯ n−1 EQt × exp b¯ ′n−1,h ht +1 + b¯ ′n−1,g gt +1

C.1. Proof of Proposition 1

exp −δ0 − δ1′ ,h ht − δ1′ ,g gt

Q

Q

Appendix C. Factor rotations

an n-period bond whose price in the next period is Ptn+−11 . We find Q

Q

Under these restrictions, the process h˜ ′t , g˜t′

= exp −δ0 − δ1′ ,h ht − δ1′ ,g gt

such that a¯ 1 = −δ0 , b¯ 1,g = −δ1,g and b¯ 1,h = −δ1,h . Next, consider Ptn = Et

′ ¯ for i = 1, . . . , H and b¯ n−1,gh = Σgh bn−1,g + b¯ n−1,h . The Laplace ′ ′¯ transform exists only if ei Σh bn−1,gh < 1 for i = 1, . . . , H.

−1

′ where ugh = Σgh ug + uh is an H × 1 vector. The Laplace transform exists only if e′i Σh′ ugh < 1 for i = 1, . . . , H. This is the key expression for solving for closed-form zero-coupon bond prices.

Q

Q

where Σg Σg′ is a GH × GH matrix with diagonal elements Σi,g Σi′,g

′

e′i Σh−1 Φh + b¯ ′n−1,g Φgh − Σgh Φh

′ 1 − δ1,h + IH ⊗ b¯ n−1,g Σg Σg′ ιH ⊗ b¯ n−1,g ht 2

i =1

e′i Σh′ ugh

1 − ei Σh bn−1,gh ′¯

′

νh,i log 1 − e′i Σh′ ugh

H

75

e′i Σh′ b¯ n−1,gh ′

i=1

2

−

H

Q

ν −νh,i wi,ht ,+i 1

76

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

For intuition, break the log-SDF mt +1 = log (Mt +1 ) into three terms; one for each of the shocks

=

H

−e i Σ h Λ h Σ h w t −1

′

zi,t +1 − e′i Σh−1 Φh Σh wt e′i Σh−1 Φh Σh wt

i=1

mt +1 = −rt + mg ,t +1 + mh,t +1 + mz ,t +1

=

where mi,t +1 is the compensation for risk i. Let λg = µg − µg , Q λh = νh −νhQ , Λg = Φg − ΦgQ , Λh = Φh − ΦhQ , and Λgh = Φgh − Φgh . Q

For the Gaussian part, we find mg ,t +1 = − 21 λ′gt λgt − λ′gt ϵg ,t +1

Σg−,t1 εg ,t +1 is a standard, zero mean Gaussian shock.

where ϵg ,t +1 = The price of Gaussian risk is

λgt = Σg−,t1

λg + Λg gt + Λgh ht − Σgh [Σh λh + Λh (ht − µh )] .

This generalizes the expression for Gaussian ATSMs, with a timevarying quantity of risk Σg ,t . Recall that wt +1 = Σh−1 (ht +1 − µh ). We will write risk compensation in terms of wt +1 . mh,t +1 =

H i=1

log Γ νh,i + zi,t +1 − νh,i + zi,t +1 − 1

is Poisson r.v. standardized to have

−1 e′i Σh Φh Σh wt

−1 e′ Σh Λh Σh wt

mean 0 and variance 1 and λzt ,i = i

Appendix E. Log-likelihood for spanned models

T 1

2 t =2 1

−

2 t =2

tr

T H

−

T 1

tr Ω −1 ηt ηt′ −

T

2 t =2

Σg ,t −1 Σg′ ,t −1

−1

e′i Σh−1 (ht − µh ) −

t =2 i=1

=

wi,t +1 − νh,i − zi,t +1 log 1 + νh,i + zi,t +1

−λh,i

i=1

≈

2

t =2 i=1

i=1

where ϵw,t +1,i =

wi,t +1 −νh,i −zi,t +1 √ νi +zi,t +1

εgt εgt′ − (T − 1) log |Σh | T H

log e′i Σh−1 [ht − µh ] −

is a gamma r.v. stan-

νh,i +zi,t +1

t =2 i =1

Q

H i =1

−1

−1

Q h

−zi,t +1

e′i Σh−1 Λh Σh wt ei Σh Φh Σh wt ′

−1

′

−1

−1

dℓ θˆc (θm ),θm ′ dθm

Proof. where

=

dℓ θˆc (θm ),θm

∂ℓ θˆc ,θm

′ dθm ∂ℓ(θc ,θm ) |θc =θˆc ∂θc′

=

, where θˆc (θm ) = arg maxθc ℓ (θc , θm ).

′ ∂θm ∂ℓ θˆc ,θm ′ ∂θm

+

ˆ ∂ℓ(θc ,θm ) |θc =θˆc dθdc θ(θ′ m ) ∂θc′ m

= 0 by the definition of θˆc .

=

∂ℓ θˆc ,θm ′ ∂θm

,

F.1. Proof of Proposition 2 Proof. First, we show that optimizing the original likelihood and optimizing the concentrated likelihood lead to equivalent ∂ℓ(θc ,θm ) solutions. Solving maxθ ℓ (θ) requires = 0 and ∂θ ′

+ e′i Σh−1 Λh Σh wt

22 Our derivation uses the approximation that Γ (a+x) ∝ xa−b 1 + O Γ (b+x) x. We also use the fact that log(1 + x) = x for small x.

ei Σh Φh [ht −1 − µh ] ′

Lemma 1. The derivative of the concentrated log-likelihood function can be computed as the partial derivative of the log-likelihood func-

Φh (ht − µh )

− ei Σh Φ (ht − µh ) − zi,t +1 log ei Σh + log zi,t +1 ! + e′i Σh−1 Φh (ht − µh ) H e′i Σh−1 Λh Σh wt = zi,t +1 log 1 − + e′i Σh−1 Λh Σh wt e′i Σh−1 Φh Σh wt i =1 ≈

′

In this appendix, we prove the propositions in the paper. We begin with a preliminary lemma.

i =1

′

2

Appendix F. Proofs and analytical derivatives

tion:

zi,t +1 log e′i Σh−1 Φh (ht − µh ) − log zi,t +1 !

T H νh,i − 1

where Iλ (x) is the modified Bessel function of the first kind and ei is the H × 1 unit vector.

Consider the non-Gaussian part due to the Poisson distribution H

e′i Σh−1 Φh (ht −1 − µh )

ei Σh [ht − µh ]

2

dardized to have mean 0 and variance 1.22 The market price of risk λ is λwt ,i = √ h,i . We note that V (wit |zt ) = νh,i + zi,t .

m z ,t + 1 =

t =2 i =1

The compensation for gamma risks is approximately mh,t +1 ≈

−λ′wt ϵw,t +1

log |Ω |

T H log Iνh,i −1 × log e′i Σh−1 Φh [ht −1 − µh ] +

λh,i wi,t +1 − νh,i − zi,t +1 −√ . √ νh,i + zi,t +1 νh,i + zi,t +1

H

2

log Σg ,t −1 Σg′ ,t −1

i=1

T −1

t =2 i =1

T H νh,i − 1

+

−λh,i log wi,t +1 − log νh,i + zi,t +1

H

.

−1 e′i Σh Φh Σh wt

i=1 H

−1 zi,t +1 −e′i Σh Φh Σh wt

ϵ z ,t + 1 =

−

× log e′i Σh−1 (ht +1 − µh ) + e′i Σh−1 (ht +1 − µh ) H Γ νh,i + zi,t +1 Q − νh,i − νhQ,i = log Γ νh,i + zi,t +1 i=1 ′ −1 × log ei Σh (ht +1 − µh ) H νh,i −νhQ,i ≈ log νh,i + zi,t +1 − λh,i log wi,t +1 =

ϵz ,t +1,i .

The log SDF is mz ,t +1 ≈ − λ′zt ϵz ,t +1 where

i=1

e′i Σh−1 Φh Σh wt

i=1

ℓ(θ ) = CONST − (T − 1) log |B1 | −

× log e′i Σh−1 (ht +1 − µh ) − e′i Σh−1 (ht +1 − µh ) ×

e′ Σh−1 Λh Σh wt

− i

Dropping the initial condition, the conditional log-likelihood for the general affine model is given by

− log Γ νhQ,i + zi,t +1 + νhQ,i + zi,t +1 − 1

H

H

1 x

for large

∂ℓ(θm ,θc ) = 0. The ∂θc′ dℓ θˆc (θm ),θm

erty

′ dθm

m

solution to maxθm ℓ θˆc (θm ) , θm has the prop-

= 0. By Lemma 1,

∂ℓ θˆc ,θm ′ ∂θm

= 0. This with

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

∂ℓ(θc ,θm ) |θc =θˆc ∂θc′

= 0 by the definition of θˆc constitute the two F.O.C.’s

for the original problem. Second, showing θˆc can be solved by least squares given θm is straightforward. A value for θm maps into a value for the loadings A and B though (9)–(10), and therefore xt by (16). Given the factors, Eq. (11) is a linear regression with heteroskedasticity, hence GLS. Given yields, factors and bond loadings, Ω in Eq. (17) is the variance–covariance matrix in a linear regression with homoskedasticity, hence OLS. F.2. Proof of Proposition 3 We use the result in Lemma 1 to derive the gradients for the concentrated log likelihood based on the original log-likelihood, which makes the derivation easier.

Proof. A direct result of Lemma 1 is

′ dθm

∂ℓ θˆc ,θm ,A(θm ),B(θm ) ′ ∂θm

dℓ θˆc (θm ),θm ,A(θm ),B(θm )

=

. Applying the chain rule,

∂ℓ θˆc , θm , A (θm ) , B (θm ) ∂θm′ ∂ℓ θˆc , θm , A, B

dℓ

=

dvech(Σ0,g

dℓ

=

e i Σh ˆ − 1 Σh (ht − µh ), and hˆ t −1 = Σh Φh (ht −1 − µh ). The matrices Sg and Sh extract the Gaussian gt = Sg xt and non-Gaussian factors ht = Sh xt from xt . T −1 ∂ℓ = − εgt′ Σg ,t −1 Σg′ ,t −1 Σgh Σh ei ∂νh,i t =2 T T ′ˆ ∂ I h˜ it ν − 1 h, i 1 e ht 1 + + log i ˆ 2 t =2 ∂νh,i ˜ it t =2 Iνh,i −1 h e′i hˆ t −1

dvech Σi,g

ei Σh [ht − µh ]

′

−1

Φh [ht −1 − µh ] , hˆ t =

T ′ −1 ∂ℓ ′ ′ ′ = − vec Σ Σ Σ ε h g , t − 1 gt gh g , t − 1 t − 1 ∂ vec (Φh )′ t =2 T H ′ −1 − vec Σh′ ei (ht −1 − µh )′ t =2 i =1 T H ′ νh,i − 1 −1 vec Σh′ e i ( h t − 1 − µh ) ′ − ˆ t =2 i =1 2e′ hˆ i t −1 T H I h˜ it ν h,i ν −1 h ,i + + h˜ it I h˜ t =2 i =1

ν h, i − 1

′ˆ

×

2ei ht h˜ it

vec

′ −1

Σh

it

e i ( h t − 1 − µh ) ′

′

∂ℓ ∂ℓ ∂A + ′ ′ ∂ vech(Σ0,g ) ∂ A ∂ vech(Σ0,g )′

∂A ∂ℓ ∂ℓ + ′ ′ ∂ vec (Σh )′ ∂ A ∂ vec (Σh ) ′ ∂ vec B ∂ℓ + ′ ′ ∂ vec (B ) ∂ vec (Σh )′ ∂ℓ ∂ℓ ∂A ′ + ′ ∂ A ∂ vec Σgh ′ ∂ vec Σgh ′ ∂ vec B ∂ℓ + ′ ′ ∂ vec (B ) ∂ vec Σgh ′

−1

dℓ

=

)′

′ =

dvec Σgh

(x) + Iλ+1 (x), see, Abramowitz

∂ I (x)

∂ℓ ∂ℓ ∂ A + ′ ∂µ′h ∂ A ∂µ′h dℓ

We now provide the analytical gradient. Note that θˆc = (µ ˆ g, ˆ g, Φ ˆ gh , Ω ˆ ) are optimized by θˆc = argmaxθc ℓ (θc , θm ) detailed in Φ −1

I x λ

∂ vec B′ ∂ℓ ∂ℓ ∂ A ∂ℓ = + ∂δ1′ ,g ∂ A′ ∂δ1′ ,g ∂ vec (B′ )′ ∂δ1′ ,g ′ ∂ vec B ∂ℓ ∂ℓ ∂ A ∂ℓ = ′ ′ + ′ ′ ′ ∂δ1,h ∂ A ∂δ1,h ∂ vec (B ) ∂δ1′ ,h ∂ℓ ∂ℓ ∂ A ∂ℓ ∂ℓ ∂ A = ′ Q′ = ′ Q′ Q′ ∂ A ∂νh ∂ A ∂µQg ′ ∂νh ∂µg ∂ vec B′ ∂ℓ ∂ℓ ∂A ∂ℓ = + ′ ∂ A′ ∂ vec ΦgQ ′ ∂ vec (B′ )′ ∂ vec ΦgQ ′ ∂ vec ΦgQ ′ ∂ vec B ∂ℓ ∂ℓ ∂ℓ ∂A Q ′ = ′ Q ′ + Q ′ ′ ′ ∂ A ∂ vec Φgh ∂ vec (B ) ∂ vec Φgh ∂ vec Φgh ∂ vec B′ ∂ℓ ∂ℓ ∂A ∂ℓ ′ = ′ + . ∂ A ∂ vec ΦhQ ′ ∂ vec (B′ )′ ∂ vec ΦhQ ′ ∂ vec ΦhQ

dℓ

′

λ

∂ℓ ∂ A ∂ℓ = ′ ∂δ0 ∂ A ∂δ0

dvec (Σh )′

=

λ and Stegun (1964). The derivative ∂λ is a complicated expression that is easier to compute numerically. The derivatives for parameters entering the loadings are calculated by the chain rule. Take derivatives of ℓ w.r.t. the loadings A1 , A2 , B1 and B2 . Then, take derivatives of the loadings w.r.t. the parameters.

dµ′h

F.3. Gradients

Let h˜ it = 2

∂ Iλ (x) ∂x

where we have used

The derivatives of the remaining parameters that enter both the loadings and the P dynamics are

∂ℓ θˆc , θm , A, B ∂ A (θ ) m + = ∂θm′ ∂ A′ ∂θm′ ∂ℓ θˆc , θm , A, B ∂ vec B (θm )′ . + ∂ vec (B′ )′ ∂θm′

Proposition 2.

77

′ =

∂A ∂ℓ ∂ℓ ′ + ′ ∂ A ∂ vech Σi,g ′ ∂ vech Σi,g ′ ∂ vec B dℓ + ∂ vec (B′ )′ ∂ vech Σi,g ′

T ∂ℓ ′ vec(Σgh (Σg ,t −1 Σg′ ,t −1 )−1 εgt νh′ )′ ′ =− ∂ vec (Σh ) t =2

− (T − 1)vec

′ ′ Σh−1

+

T H

vec

t =2 i =1

+

T H

vec

−1 ′ ˆ ′ Σh ei hˆ t −1

′

t =2 i =1

T H ′ ′ νh,i − 1 vec Σh−1 ei hˆ ′t − 2e′i hˆ t t =2 i =1 ′ T H −1 ′ ˆ ′ νh,i − 1 ˆ + vec Σh ei h t − 1 ˆ t =2 i =1 2e′ hˆ i t −1

′ ′ Σh−1 ei hˆ ′t

78

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

T H I h˜ it ν h, i νh,i − 1 − + h˜ it Iνh,i −1 h˜ it t =2 i =1 ×

2e′i hˆ t h˜ it

′ ˆ ′ vec Σh−1 ei hˆ t −1

−

T H

t =2 i=1

′ −

T H t =2 i=1

T H Iνh,i h˜ it νh,i − 1 − + h˜ it Iνh,i −1 h˜ it t =2 i =1 ˆ

×

2e′i hˆ t −1 h˜ it

vec

∂ℓ ′ = ∂ vec Σgh

vec

Σg ,t −1 Σg′ ,t −1

−1

′ εgt εht

′

T −1 ∂ℓ εgt εgt′ − IG vec Σg ,t −1 Σg′ ,t −1 ′ = ∂ vech Σ0,g t =2 ′ −1 × Σg ,t −1 Σg′ ,t −1 Σ0,g DGL

−1 ∂ℓ Σg ,t −1 Σg′ ,t −1 vec εgt εgt′ − IG ′ = ∂ vech Σi,g t =2 ′ −1 × Σg ,t −1 Σg′ ,t −1 Σi,g hi,t −1 DGL T

−1 ∂ℓ εgt′ Σg ,t −1 Σg′ ,t −1 Σgh (IH − Φh ) ′ = − ∂µh t =2 T H T H + e′i Σh−1 + e′i Σh−1 Φh t =2 i =1 t =2 i=1 T H T H νh,i − 1 ′ −1 νh,i − 1 ′ −1 ei Σh Φh e i Σh + − ˆ 2e′i hˆ t t =2 i=1 t =2 i =1 2e′i hˆ t −1 T H Iνh,i h˜ it νh,i − 1 − + h˜ it Iνh,i −1 h˜ it t =2 i =1 2 ′ ˆ × ei hˆ t e′i Σh−1 Φh + e′i hˆ t −1 e′i Σh−1 h˜ it T

T ′ ∂ℓ vec xt ηt′ Ω −1 ′ ′ =

∂ vec B2

2

H

1 vec Σi,g Σi′,g Shi B− 1 +

i =1

T H

1 e′i Σh−1 Sh B− 1

t =2 i =1

νh,i − 1 1 e′i Σh−1 Sh B− − 1 ′ˆ 2e h t t =2 i =1 i T

+

H

T H t =2 i =1

ei Σh Φh Sh B1 ′

−1 − vec xt −1 εgt′ Σg ,t −1 Σg′ ,t −1 1 ′ × Φg Sg + Φgh − Σgh Φh Sh B− 1

−1

−1

1

+ vec 2

×

T H νh,i − 1 ′ −1 1 + ei Σh Φh Sh B− 1 ˆ ′ ˆ t =2 i =1 2e h i t −1

H

−1 −1 ′ ′ Σg ,t −1 Σg′ ,t −1 IG − εgt εgt Σg ,t −1 Σg′ ,t −1

vec Σi,g Σi′,g

1 ′ Shi B− 1 ⊗ x t −1

i=1

+

T H

1 vec xt e′i Σh−1 Sh B− 1

′

t =2 i=1

T H νh,i − 1 1 ′ − vec xt e′i Σh−1 Sh B− 1 ′ˆ 2ei ht t =2 i=1 +

T H

1 vec xt −1 e′i Σh−1 Φh Sh B− 1

′

t =2 i=1

T H νh,i − 1 1 ′ vec xt −1 e′i Σh−1 Φh Sh B− + 1 ˆ t =2 i=1 2e′ hˆ i t −1

−

T H

t =2 i=1

Iνh,i h˜ it

νh,i − 1 + h˜ it Iνh,i −1 h˜ it

ˆ

×

t =2

T −1 ∂ℓ 1 ′ ′ = −ηt′ Ω −1 B2 B− 1 + εgt Σg ,t −1 Σg ,t −1 ′ ∂ A1 t =2 1 × IG − Φg Sg − Φgh + Σgh − Σgh Φh Sh B− 1 −1 ′ 1 −1 ′ IG − εgt εgt Σg ,t −1 Σg′ ,t −1 + vec Σg ,t −1 Σg′ ,t −1

×

Iνh,i h˜ it νh,i − 1 2e′ hˆ t 1 i e′i Σh−1 Φh Sh B− + 1 h˜ it h˜ it Iνh,i −1 h˜ it

∂ vec B1 t =2 1 ′ −1 Sg − Σgh Sh B− + vec xt εgt′ Σg ,t −1 Σg′ ,t −1 1

t =2

T ∂ℓ = ηt′ Ω −1 ∂ A′2 t =2

Iνh,i h˜ it

ˆ νh,i − 1 2e′ hˆ t −1 ′ −1 1 i + ei Σh Sh B− 1 h˜ it h˜ it Iνh,i −1 h˜ it

T 1 ′ ∂ℓ 1 ′ + −vec xt ηt′ Ω −1 B2 B− ′ ′ = −(T − 1)vec B− 1 1

′ ′ Σh−1 ei hˆ ′t

T

−

2e′i hˆ t −1 h˜ it T H t =2 i=1

×

2e′i hˆ t h˜ it

1 vec xt e′i Σh−1 Sh B− 1

′

Iνh,i h˜ it νh,i − 1 + h˜ it Iνh,i −1 h˜ it

′

1 vec xt −1 e′i Σh−1 Φh Sh B− . 1

The derivatives of A and B w.r.t. the parameters can be computed recursively as a function of maturity along with a¯ n , b¯ n,g and b¯ n,h . The derivatives of Bg and Bh have separate recursions. Let b¯ n,g ,ψ denote the derivatives of the Gaussian loadings b¯ n,g at maturity n w.r.t. a parameter ψ . All recursions are written assuming that ψ is a full vector/matrix of parameters with no restrictions. If the matrix has fewer parameters than entries, the user will have to multiply the respective recursion by a selection matrix. Let d¯ n−1 = diag ιH − Σh′ b¯ n−1,gh be a diagonal H × H matrix. Let c¯n′ −1

=

′ Q′ Q′ 1 −1′ ¯ −2 ν Q′ d¯ − n−1 − µh Φh Σh dn−1 Σh be an 1 × H

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81 Q′ h

b¯ n,h,Σgh = Φ

vector.

79

2 ′¯ Σh−1′ d¯ − n−1 Σh bn−1,h,Σgh

H ×GH

= a¯ ′n−1,µQ + b¯ ′n−1,g

a¯ ′

Q n,µg

g

b¯ n,h,δ1g

1×G

′

+ b¯ ′n−1,h + b¯ ′n−1,g Σgh ΦhQ

′

a¯ n,µh = a¯ n−1,µh

H ×G

1×H 1 −1 Q − b¯ ′n−1,gh Σh d¯ − n−1 Σh Φh ′ ¯′ ¯′ = a¯ ′ Q − log ιH − bn−1,gh Σh − bn−1,g Σgh Σh

a¯ ′

Q n,νh

n−1,νh

b¯ n,h,Φ Q g

′ 2 ¯′ − ΦhQ′ Σh−1′ IH − d¯ − n−1 Σh ⊗ bn−1,g Q Q ′ ¯ 2 ′¯ bn−1,g ,δ1g = ΦhQ′ Σh−1′ d¯ − n−1 Σh bn−1,gh,δ1g + Φgh − Σgh Φh + IH ⊗ b¯ ′n−1,g Σg Σg′ ιH ⊗ b¯ n−1,g ,δ1g Q Q ′ ¯ 2 ′¯ Q = ΦhQ′ Σh−1′ d¯ − bn−1,g ,Φ Q n−1 Σh bn−1,gh,Φg + Φgh − Σgh Φh g

H × G2

+ IH ⊗ b¯ ′n−1,g Σg Σg′ ιH ⊗ b¯ n−1,g ,ΦgQ

1×H

′

= a¯ ′n−1,Σ0,g + vec b¯ n−1,g b¯ ′n−1,g Σ0,g DGL

a¯ ′n,Σ0,g

1×G(G+1)/2

2 ′¯ = ΦhQ′ Σh−1′ d¯ − n−1 Σh bn−1,h,Σh 1 −1 − ΦhQ′ Σh−1′ ⊗ b¯ ′n−1,gh Σh d¯ − n − 1 Σh 2 ¯′ + ΦhQ′ Σh−1′ d¯ − n−1 ⊗ bn−1,gh

¯

bn,h,Σh H ×H (H +1)/2

a¯ ′n,δ1h = a¯ ′n−1,δ1h + µ′h b¯ n−1,h,δ1h + c¯n′ −1 b¯ n−1,h,δ1h 1×H

a¯ ′

= a¯ ′n−1,Φ Q + µ′h b¯ n−1,h,Φ Q + c¯n′ −1 b¯ n−1,h,Φ Q

Q

n,Φgh

gh

gh

gh

1×GH

= a¯ ′n−1,Σi,g + µ′h b¯ n−1,h,Σi,g + c¯n′ −1 b¯ n−1,h,Σi,g

a¯ ′n,Σi,g

1×G(G+1)/2

i = 1, . . . , H Q n,Φh 2 1×H

h

h

′ ¯ b¯ n−1,gh,Φ Q = Σgh bn−1,g ,Φ Q + b¯ n−1,h,Φ Q . g

h

a¯ n,Σgh 1×GH

a¯ n,δ1g

Q′ h

′ ′ + µQg ′ + µ′h ΦhQ′ Σgh − νhQ′ Σh′ Σgh b¯ n−1,g ,δ1g + b¯ ′n−1,g Σ0,g Σ0′ ,g b¯ n−1,g ,δ1g = a¯ ′n−1,Φ Q + µ′h b¯ n−1,h,ΦgQ + c¯n′ −1 db¯ n−1,gh,ΦgQ

Q n,Φg 1×G2

g

+ µ + µh Φ Σgh − ν Σh Σgh Q + b¯ ′n−1,g Σ0,g Σ0′ ,g b¯ Q′ g

′

Q′ h

′

Q′ h

′

′

b¯

Q

n−1,g ,Φg

n−1,g ,Φg

a¯ ′n,Σh

1×H (H +1)/2

Many of the derivatives of the loadings are zero for all maturities including b¯ n,g ,µh , b¯ n,h,µQ , b¯ n,g ,δ1h , b¯ n,g ,Φ Q , b¯ n,g ,Σh , b¯ n,g ,Σ0,g , g G× H

gh

H ×G(G+1)/2

1×G

a¯ ′

g

h

G× H

H ×G

= a¯ ′n−1,Σh + µ′h b¯ n−1,h,Σh + c¯n′ −1 b¯ n−1,h,Σh ′ 1 + vec Σh−1′ d¯ − Σ ′ b¯ µ′ Φ Q′ Σ −1′ ′ n−1 ′h′ n−1,gh h h h ′ ′ − vec Σgh b¯ n−1,g νh + vec b¯ n−1,gh c¯n−1 Σh−1′

G×GH

H ×G(G+1)/2

and = ιN and with b¯ 1,g ,δ1g = −IG , b¯ 1,h,δ1h = −IH . All other initial conditions start at zero.

h

G×GH

G×H

G.1. Proof of Proposition 4 We provide proofs for the U1 (4)(φ, φ 2 , ψ) and U1 (4)(φ, φ 2 ,

φ 4 ) models.

Proof for U1 (4)(φ, φ 2 , ψ): Showing bn,h = 0 is equivalent to ¯bn,h = 0. We prove b¯ n,h = 0 by induction over n. At n = 1, we have b¯ 1,h = −δ1,h = 0. Next, suppose b¯ n,h = 0, then under the restriction that Σgh = 0, we find that b¯ n,gh = 0. Imposing the restrictions on Σ1,g Σ1′ ,g and from (6), the loading recursion reduces to

b¯ n+1,h = Φgh,1 b¯ n,g ,1 + Φgh,2 b¯ n,g ,2 + Φgh,3 b¯ n,g ,3 + Q

Q

G×G

= ΦgQ′ b¯ n−1,g ,ΦgQ + IG ⊗ b¯ ′n−1,g

b¯ n,h,δ1h = Φ

2 ′¯ Σh−1′ d¯ − n−1 Σh bn−1,h,δ1h

Q

− IH Q − Φgh ,3

H ×H Q′ h

b¯ n,h,Φ Q = Φ gh

2 ′¯ Q Σh−1′ d¯ − n−1 Σh bn−1,h,Φ gh

¯′

+ IH ⊗ bn−1,g

H ×GH

b¯ n,h,Σi,g

Q′ h

= Φ

H ×G(G+1)/2

2 ′¯ Σh−1′ d¯ − n−1 Σh bn−1,h,Σgi

′ + ei vec b¯ n−1,g b¯ ′n−1,g Σi,g DGL

1 − φn 1−φ

b¯ 2n,g ,1 Σ12,g ,11 .

2n

i = 1, . . . , H

Q′

h

1 −1 − IH ⊗ b¯ ′n−1,g Σgh + IH ⊗ b¯ ′n−1,gh Σh d¯ − n − 1 Σh

Q δ1,g ,1 − Φgh ,2

1 − ψn 1−ψ

δ1,g ,3 +

1 − φ 2n 1 − φ2

1 (1 − φ n )2 2 (1 − φ)2

δ1,g ,2 δ12,g ,1 Σ12,g ,11 .

Collect terms related to (φ n )0 , (φ n )1 , (φ n )2 , ψ n , we get the following three equations

2 ′¯ Q b¯ n,h,Φ Q = Φh Σh−1′ d¯ − n−1 Σh bn−1,h,Φ h H ×H 2

2

Q

b¯ n+1,h = −Φgh,1

G×G2 Q′ h

1

The restrictions on Φg together with (7) implies the solution n

Q

h

H ×H

1−φ 1−φ for b¯ n,g : b¯ n,g ,1 = − 1−φ δ1,g ,1 , b¯ n,g ,2 = − 1−φ 2 δ1,g ,2 , b¯ n,g ,3 = n − 11−ψ δ . Substituting these above gives −ψ 1,g ,3

b¯ n,g ,δ1g = ΦgQ′ b¯ n−1,g ,δ1g − IG n,g ,Φg

G×G(G+1)/2

Appendix G. USV restrictions

Q

∂ vec(A) ∂δ0

b¯

G× H 2

G× H 2

b¯ n,h,Σ0,g , b¯ n,g ,Φ Q , b¯ n,g ,Σi,g , b¯ n,g ,Σgh , b¯ n,g ,ν Q , b¯ n,h,ν Q .

′

+ vec b¯ n−1,g µ′h Φ − ν Σh′ + c¯n′ −1 = a¯ ′n−1,δ1g + µ′h b¯ n−1,h,δ1g + c¯n′ −1 b¯ n−1,gh,δ1g Q′ h

′

g

H × G2

′ ′ ′ 1 ′¯ + vec Σgh b¯ n−1,g − Σh−1′ d¯ − n−1 Σh bn−1,gh µh = a¯ ′n−1,Σgh + µ′h b¯ n−1,h,Σgh + c¯n′ −1 b¯ n−1,h,Σgh

′

′ ¯ bn−1,g ,δ1g + b¯ n−1,h,δ1g b¯ n−1,gh,δ1g = Σgh H ×G

= a¯ ′n−1,Φ Q + µ′h b¯ n−1,h,Φ Q + c¯n′ −1 b¯ n−1,h,Φ Q

a¯ ′

where we also need to account for the derivatives of b¯ n,gh = ′ ¯ Σgh bn,g + b¯ n,h as

b¯ n+1,h =

δ1,g ,1 δ1,g ,2 δ1,g ,3 Q Q − Φgh − Φgh ,2 ,3 2 1−φ 1−φ 1−ψ 2 1 δ1,g ,1 δ1,g ,3 Q + Σ12,g ,11 + Φgh,3 ψn 2 2 (1 − φ) 1−ψ Q − Φgh ,1

80

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

δ12,g ,1 δ1,g ,1 2 + Φ − Σ1,g ,11 φ n 1−φ (1 − φ)2 2 δ1,g ,2 1 δ1,g ,1 Q 2 Σ1,g ,11 φ 2n . + Φgh,2 + 1 − φ2 2 (1 − φ)2

quantity Q θ |θ (i) = Q1 θm,b , θc |θ (i) + Q2 θm,h |θ (i) are

Q gh,1

−

(φ n )0 , (φ n )1 , (φ n )2 , ψ n are all zero. Hence, b¯ n+1,h = 0. Proof for U1 (4)(φ, φ 2 , φ 4 ): We prove b¯ n,h = 0 by induction. First, at maturity n = 1, we have b¯ 1,h = −δ1,h = 0. Next, suppose b¯ n,h = 0, then b¯ n,gh = 0 because Σgh = 0 under the USV restrictions. From (6), the non-Gaussian loading simplifies to 1 + b¯ ′n,g Σ1,g Σ1′ ,g b¯ n,g . 2

−

T 1

2 t =2 T 1

2 t =2

−

T −1 H

+

T −1 H νh,i − 1

−

T −1 H νh,i − 1

+

H T −1

The parameter restrictions on Φg together with (7) implies the Q

δ1,g ,2 , b¯ n,g ,3

4n

δ . Substituting these in, we find = − 11−φ −φ 4 1,g ,3 b¯ n+1,h = −Φgh,1 Q

1−φ

Q − Φgh ,3

1 − φ 4n 1 − φ4

1 1 − φ 2n

+

Q δ1,g ,1 − Φgh ,2

2 1 − φ2

δ1,g ,3 +

1 − φ2

1 (1 − φ n )2 2 (1 − φ)2

Q −Φgh ,1

+

δ

2 1 ,g ,1

1

2 (1 − φ)

2

Φ

Q gh,1

+

Φgh,2 Q

+

−

δ1,g ,2 δ12,g ,1 Σ12,g ,11

+

Φ

δ

1

Q − Φgh ,3

2 1,g ,2

2 1 − φ2

δ1,g ,3 1 − φ4

2 Σ1,g ,22 2

2

t =1 i=1

E log e′i Σh ht −1

E log e′i Σh Φh ht −1

E log Iνh,i −1

−1

2

e i Σh h t ′

−1

ei Σh Φh ht −1

′

−1

.

(1)

and Pitt (2011) to maximize the log-likelihood log p Y1:T ; θ over θm,h .

For the results in the paper, we used the third option. H.2. Particle filter

2 2n 2 Σ1,g ,22 φ

particle filter for the E-step in (a.). Let We implement a basic (m)

q ht |ht −1 , gt +1 , gt , θ

δ1,g ,3 1 δ 2 4n + Σ1,g ,22 φ 4 1−φ 2 1 − φ2 2 2 1 ,g ,2

e′i Σh−1 E [ht ]

expectations of zt and ht . By introducing zt , the maximization of Φh is analytical. • Option #3: Instead of optimizing over Q2 θm,h |θ (i) at each iteration, a valid EM algorithm can optimize the log-likelihood over θm,h . When H = 1, we can use the particle filterof Malik

2 δ1,g ,2 1 δ1,g ,1 + Σ12,g ,11 2 1−φ 2 (1 − φ)2

1 − φ2 Q gh,3

Σ12,g ,11 +

δ1,g ,2 1 − φ2

T −1 H

• Option #1: Maximize the intermediate quantity Q2 θm,h |θ (i) numerically over θm,h as above. • Option #2: From the definition of ht +1 in (2)–(4), there is also the variable zt . When calculating Q2 θm,h |θ (i) , take the

δ12,g ,1 δ1,g ,1 2 − Σ1,g ,11 φ n 1−φ (1 − φ)2

δ12,g ,2

Q − Φgh ,2

2

2 2 2 δ1,g ,2 Σ1,g ,22 .

δ1,g ,1

Maximizing Q1 θm,b , θc |θ (i) is similar to estimation of a Gaussian ATSM. When maximizing Q1 θm,b , θc |θ (i) , the analytical gradient of the intermediate quantity can be derived from the gradients in Appendix E. There are several options for maximizing Q2 θm,h |θ (i) , all of which lead to different types of EM algorithms. Each EM algorithm leads to the same maximum but they converge at different rates.

2

1−φ

εgt εgt′

t =1 i =1

Collect terms related to (φ n )0 , (φ n )1 , (φ n )2 , (φ n )4 b¯ n+1,h =

−1

e′i Σh−1 Φh E [ht −1 ]

t =1 i =1

2

1 − φ 2n

Σg ,t −1 Σg′ ,t −1

t =1 i=1

2

+ b¯ 2n,g ,2 Σ12,g ,22 . 1−φ 2n 1−φ 2

2 t =2

E log Σg ,t −1 Σg′ ,t −1

1

1 − φn

tr E

T 1

log |Ω |

t =1 i=1

(G.1)

δ1,g ,1 , b¯ n,g ,2 = −

2

and

1 Q Q Q b¯ n+1,h = Φgh,1 b¯ n,g ,1 + Φgh,2 b¯ n,g ,2 + Φgh,3 b¯ n,g ,3 + b¯ 2n,g ,1 Σ12,g ,11

solution for b¯ n,g : b¯ n,g ,1 = −

Q2 θm,h |θ (i) = −(T − 1) log |Σh | −

Substituting the parameter restrictions on Σ1,g Σ1′ ,g into the equation, we find

1−φ n 1−φ

tr Ω −1 ηt ηt′ −

T −1

Q

b¯ n+1,h = Φ ¯

Q1 θm,b , θc |θ (i) = − (T − 1) log |B1 | −

The restrictions on Φgh guarantee that the coefficients in front of

Q′ gh bn,g

heavier than the target distribution.

The restrictions on Φgh guarantee that the coefficients in front of Q

(φ n )0 , (φ n )1 , (φ n )2 , (φ n )4 are all zero. Therefore, b¯ n+1,h = 0.

be an importance density whose tails are

Appendix H. EM algorithm H.1. Intermediate quantity Given the identifying restriction µh = 0, we drop this parameter for convenience. The two components of the intermediate

For t = 1, . . . , T , run: (m) • For m = 1, . . . , M, draw from a proposal distribution: ht ∼ (m) q ht |ht −1 , gt +1 , gt , θ . (m) • For m = 1, . . . , M, calculatethe importance weight: wt ∝ m) w ˆ t(− 1

(m) (m) (m) p gt +1 |gt ,ht ,θ p ht |ht −1 ,θ (m) (m) q ht |ht −1 ,gt +1 ,gt ,θ

.

• For m = 1, . . . , M, normalize the weights: w ˆ t(m) =

(m)

wt M

m=1

(m) .

wt

D.D. Creal, J.C. Wu / Journal of Econometrics 185 (2015) 60–81

• Calculate the effective sample size: ESSt =

1

M

m=1

. (m) 2 w ˆt

M • If ESSt < 0.5M resample h(t m) with probabilities m=1 M w ˆ t(m) and set w ˆ t = 1/M. m=1

At time t = 1, q (h1 ; θ) does not depend on any previous particles. Simple proposal distributions are to draw from the transition density p (ht +1 |ht ; θ ) of the model (2)–(4). To calculate the expectations within Q θ |θ (i) , we use the algorithm of Godsill et al. (2004) that draws samples from the pos-

(m)

terior. Store w ˆt

, h(t m)

M m=1

for t = 1, . . . , T during the particle (m)

filter. On a backwards pass, sample h˜ T = hT and then for t = T − 1, . . . , 1.

(m)

with probability w ˆT

m) (m) • For m = 1, . . . , M, calculate backwards weights wt(+ 1|t ∝ wt (m) p h˜ t +1 |ht ; θ .

• For m = 1, . . . , M, normalize the weights: m) w ˆ t(+ 1|t =

(m)

wt +1|t M (m) m=1 wt +1|t (m)

• Sample h˜ t = ht

. (m)

with probability w ˆ t +1|t .

To calculate the expectations Q θ |θ (i) , we repeat this backwards pass a large number of times. For the final climb after the EM algorithm, we optimize over all θ using the algorithm from Malik and Pitt (2011). To implement this particle filter, we resample at every time period and use the resampling algorithm described in their appendix.

References Abramowitz, M., Stegun, I.A., 1964. Handbook of Mathematical Functions. Dover Publications Inc., New York, NY. Adrian, T., Crump, R.K., Moench, E., 2012. Pricing the term structure with linear regressions. J. Financ. Econ. 110 (1), 110–138. Aït-Sahalia, Y., 2008. Closed-form likelihood expansions for multivariate diffusions. Ann. Statist. 36, 906–937. Aït-Sahalia, Y., Kimmel, R.L., 2010. Estimating affine multifactor term structure models using closed-form likelihood expansions. J. Financ. Econ. 98, 113–144. Andersen, T., Benzoni, L., 2010. Do bonds span volatility risk in the U.S. treasury market? A specification test for affine term structure models. J. Finance 65 (2), 603–653. Ang, A., Piazzesi, M., 2003. A no-arbitrage vector autoregression of term structure dynamics with macroeconomic and latent variables. J. Monetary Econ. 50, 745–787. Bauer, M.D., (2014) Bayesian Estimation of Dynamic Term Structure Models Under Restrictions on Risk Pricing. Federal Reserve Bank of San Francisco Working Paper. Bauer, M.D., Rudebusch, G.D., Wu, J.C., 2012. Correcting estimation bias in dynamic term structure models. J. Bus. Econom. Statist. 30 (3), 454–467. Cheridito, P., Filipovic, D., Kimmel, R.L., 2007. Market price of risk specifications for affine models: theory and evidence. J. Financ. Econ. 84 (1), 123–170. Christensen, J.H., Diebold, F.X., Rudebusch, G.D., 2011. The affine arbitrage-free class of Nelson–Siegel term structure models. J. Econometrics 164 (1), 4–20. Cochrane, J.H., Piazzesi, M., (2008) Decomposing the Yield Curve. Booth School of Business, University of Chicago, Unpublished Manuscript. Collin-Dufresne, P., Goldstein, R.S., 2002. Do bonds span the fixed income markets? Theory and evidence for unspanned stochastic volatility. J. Finance 57 (4), 1685–1730. Collin-Dufresne, P., Goldstein, R.S., Jones, C., 2008. Identification of maximal affine term structure models. J. Finance 63 (2), 743–795. Collin-Dufresne, P., Goldstein, R.S., Jones, C., 2009. Can the volatility of interest rates be extracted from the cross section of bond yields? An investigation of unspanned stochastic volatility. J. Financ. Econ. 94 (1), 47–66. Cox, J.C., Ingersoll, J.E., Ross, S.A., 1985. A theory of the term structure of interest rates. Econometrica 53 (2), 385–407.

81

Creal, D.D., 2012. A survey of sequential Monte Carlo methods for economics and finance. Econometric Rev. 31 (3), 245–296. Creal, D.D., Koopman, S.J., Lucas, A., 2011. A dynamic multivariate heavy-tailed model for time-varying volatilities and correlations. J. Bus. Econom. Statist. 29 (4), 552–563. Creal, D.D., Koopman, S.J., Lucas, A., 2013. Generalized autoregressive score models with applications. J. Appl. Econometrics 28 (5), 777–795. Dai, Q., Singleton, K.J., 2000. Specification analysis of affine term structure models. J. Finance 55 (5), 1943–1978. de Jong, P., 1991. The diffuse Kalman filter. Ann. Statist. 19 (2), 1073–1083. Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. J. Roy. Statist. Soc. Ser. B 39 (1), 1–38. Diebold, F.X., Rudebusch, G.D., 2013. Yield Curve Modeling and Forecasting. Princeton University Press, Princeton, NJ. Diez de Los Rios, A., (2013) A New Linear Estimator for Gaussian Dynamic Term Structure Models, Bank of Canada, Unpublished Manuscript. Duffee, G.R., 2002. Term premia and interest rate forecasts in affine models. J. Finance 57 (1), 405–443. Duffee, G.R., 2011. Information in (and not in) the term structure. Rev. Financ. Stud. 24 (9), 2895–2934. Duffee, G.R., (2012) Bond pricing and the macroeconomy. Working Paper, Johns Hopkins University. Duffie, D., Kan, R., 1996. A yield factor model of interest rates. Math. Finance 6 (4), 379–406. Durbin, J., Koopman, S.J., 2012. Time Series Analysis by State Space Methods, second ed. Oxford University Press, Oxford, UK. Fama, E.F., Bliss, R.R., 1987. The information in long maturity forward rates. Amer. Econ. Rev. 77, 680–692. Godsill, S.J., Doucet, A., West, M., 2004. Monte Carlo smoothing for nonlinear time series. J. Amer. Statist. Assoc. 99 (465), 156–168. Gouriéroux, C., Jasiak, J., 2006. Autoregressive gamma processes. J. Forecast. 25 (2), 129–152. Gouriéroux, C., Jasiak, J., Sufana, R., 2009. The Wishart autoregressive process of multivariate stochastic volatility. J. Econometrics 150 (2), 167–181. Gouriéroux, C., Monfort, A., Polimenis, V., (2002) Affine term structure models. Institut National de La Statistique et des Etudes Economiques. Gürkaynak, R.S., Wright, J.H., 2012. Macroeconomics and the term structure. J. Econ. Lit. 50 (2), 331–367. Hamilton, J.D., Wu, J.C., 2012. Identification and estimation of Gaussian affine term structure models. J. Econometrics 168 (2), 315–331. Hamilton, J.D., Wu, J.C., 2014. Testable implications of affine term structure models. J. Econometrics 178, 231–242. Heidari, M., Wu, L., 2003. Are interest rate derivatives spanned by the term structure of interest rates?. J. Fixed Income 13, 75–86. Jacquier, E., Johannes, M., Polson, N.G., 2007. MCMC maximum likelihood for latent state models. J. Econometrics 137 (2), 615–640. Joslin, S., (2010) Can Unspanned Stochastic Volatility Models Explain the Cross Section of Bond Volatilities? Working Paper, University of Southern California, Marshall School of Business. Joslin, S., Singleton, K.J., Zhu, H., 2011. A new perspective on Gaussian affine term structure models. Rev. Financ. Stud. 27, 926–970. Kim, D.H., Orphanides, A., (2005) Term Structure Estimation with Survey Data on Interest Rate Forecasts. Federal Reserve Board, Finance and Economics Discussion Series 2005–48. Kim, D.H., Wright, J.H., (2005) An Arbitrage-free Three-factor Term Structure Model and The Recent Behavior of Long-term Yields and Distant-horizon Forward Rates. Federal Reserve Board, Finance and Economics Discussion Series 2005–33. Le, A., Singleton, K.J., Dai, Q., 2010. Discrete-time affine term structure models with generalized market prices of risk. Rev. Financ. Stud. 23 (5), 2184–2227. Li, H., Zhao, F., 2006. Unspanned stochastic volatility: evidence from hedging interest rate derivatives. J. Finance 61 (1), 341–378. Malik, S., Pitt, M.K., 2011. Particle filters for continuous likelihood evaluation and maximisation. J. Econometrics 165, 190–209. Meng, X.-L., Rubin, D.B., 1993. Maximum likelihood estimation via the ecm algorithm: A general framework. Biometrika 80 (2), 267–278. Piazzesi, M., 2010. Affine term structure models. In: Ait-Sahalia, Y., Hansen, L.P. (Eds.), Handbook of Financial Econometrics. Elsevier, New York, pp. 691–766. Richard, J.F., Zhang, W., 2007. Efficient high-dimensional importance sampling. J. Econometrics 141 (2), 1385–1411. Trolle, A.B., Schwartz, E.S., 2009. A general stochastic volatility model for the pricing of interest rate derivatives. Rev. Financ. Stud. 22 (5), 2007–2057. Wei, G.C.G., Tanner, M.A., 1990. A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms. J. Amer. Statist. Assoc. 85 (411), 699–704. White, H.L., 1982. Maximum likelihood estimation of misspecified models. Econometrica 50 (1), 1–25.