2.1. The Nelson-Siegel term structure representation

Nelson & Siegel (Reference Nelson and Siegel1987) propose a parametric representation of the yield curve of the form

(2.1) \begin{equation}y(t,T)=X^{(1)}+X^{(2)}\left(\dfrac{1-e^{-\lambda\tau}}{\lambda\tau}\right)+X^{(3)}\left(\dfrac{1-e^{-\lambda\tau}}{\lambda\tau} - e^{-\lambda\tau }\right),\end{equation}

where y(t, T) is the time-t spot rate with tenor $\tau=T-t$ , and $ X^{(1)},X^{(2)},X^{(3)}$ and $\lambda>0$ are model parameters. The model leads to a convenient interpretation of the parameters, with $X^{(1)}$ , $X^{(2)}$ , and $X^{(3)}$ , respectively, characterizing the level, slope, and curvature of the yield curve, while $\lambda$ is a dilation parameter. Although this model has first been considered in a static setting, Diebold & Li (Reference Diebold and Li2006) embed it in a dynamic setting where time-varying factors $X^{(1)}_{t},X^{(2)}_{t},X^{(3)}_{t}$ are considered instead of $X^{(1)},X^{(2)},X^{(3)}$ . Several auto-regressive models are used for the dynamics of such parameters, thereby making it possible to forecast the term structure.

A pitfall of the dynamic Nelson-Siegel model of Diebold & Li (Reference Diebold and Li2006) is its incompatibility with no-arbitrage theory as discussed Filipović (Reference Filipović1999). Christensen et al. (Reference Christensen, Diebold and Rudebusch2011) therefore adapt the Nelson-Siegel framework and develop the AFNS model, an arbitrage-free term structure model whose spot rate formulas match exactly these of the Nelson-Siegel model up to a correction term. It is obtained by specifying the short rate risk-neutral dynamics as the following affine Gaussian Itô diffusion:

(2.2) \begin{align}r_t \equiv X^{(1)}_t + X^{(2)}_t, \end{align}

(2.3) \begin{align} \left[ \begin{array}{c} dX^{(1)}_t \\[4pt] dX^{(2)}_t \\[4pt] dX^{(3)}_t \end{array} \right] \equiv \underbrace{\left[ \begin{array}{c@{\quad}c@{\quad}c} 0 & 0 & 0 \\[4pt] 0 & \lambda & -\lambda \\[4pt] 0 & 0 & \lambda \end{array} \right] }_{ \kappa^\mathbb{Q} }\underbrace{\left[ \begin{array}{c} \theta^{\mathbb{Q}}_1-X^{(1)}_t \\[4pt] \theta^{\mathbb{Q}}_2-X^{(2)}_t \\[4pt] \theta^{\mathbb{Q}}_3-X^{(3)}_t \end{array} \right] }_{ \theta^{\mathbb{Q}}-\textbf{X}_t } dt + \Sigma \!\left( {\begin{array}{c} dW^{\mathbb{Q}}_{t,1} \\[6pt] dW^{\mathbb{Q}}_{t,2} \\[6pt] dW^{\mathbb{Q}}_{t,3} \\[4pt] \end{array} } \right) \end{align}

where $\big\{W^{\mathbb{Q}}_{t,1}\big\}_{t \geq 0}$ , $\big\{W^{\mathbb{Q}}_{t,2}\big\}_{t \geq 0}$ and $\big\{W^{\mathbb{Q}}_{t,3}\big\}_{t \geq 0}$ are independent Brownian motions under the risk-neutral probability measure $\mathbb{Q}$ , $\Sigma$ is a $3\times 3$ positive semi-definite matrix and $\theta^{\mathbb{Q}}_1,\theta^{\mathbb{Q}}_2,\theta^{\mathbb{Q}}_3,\lambda$ are real numbers. As shown in Christensen et al. (Reference Christensen, Diebold and Rudebusch2011), this leads to the following formulas for spot rates:

(2.4) \begin{align}y(t,T) = -\frac{1}{\tau}\sum^3_{j=1} \tilde{B}^{(j)}(t,T) X^{(j)}_{t} - \frac{\tilde{C}(t,T)}{\tau},\end{align}

(2.5) \begin{align} \tilde{C}(t,T) = \sum^3_{j=2} \big(\kappa^\mathbb{Q} \theta^\mathbb{Q}\big)_j \int^T_{t} \tilde{B}^{(j)}(s,T)ds + \frac{1}{2} \sum^3_{j=1} \int^T_{t} \!\left(\Sigma^\top \tilde{B}(s,T) \tilde{B}(s,T)^\top \Sigma \right)_{j,j} ds\end{align}

with

(2.6) \begin{eqnarray}\tilde{B}(s,T) &\equiv& \left[\tilde{B}^{(1)}(s,T), \,\, \tilde{B}^{(2)}(s,T), \,\, \tilde{B}^{(3)}(s,T)\right]^\top \notag\\[4pt] &\equiv& \left[{-}\tau, \,\, -\left(\dfrac{1-e^{-\lambda\tau}}{\lambda} \right), \,\, \left(\tau e^{-\lambda\tau } -\dfrac{1-e^{-\lambda\tau}}{\lambda} \right)\right]^\top \end{eqnarray}

and $()_j$ or $()_{j,j}$ denoting respectively the j $^{th}$ element of a vector or the element at row j and column j of a matrix, respectively.

Striking similarities can be observed when comparing (2.1) with (2.4) and (2.6). Indeed, identifying $X^{(1)}_{t}, X^{(2)}_{t}, X^{(3)}_{t}$ with $X^{(1)}, X^{(2)}, X^{(3)}$ , discrepancies between the spot rates from both models are only caused by the corrective term $- \frac{\tilde{C}(t,T)}{\tau}$ appearing in (2.4). This leads Christensen et al. (Reference Christensen, Diebold and Rudebusch2011) to describe the arbitrage-free Nelson-Siegel model (2.2)–(2.3) as the “closest match to the Nelson-Siegel yield function” within some class of exponential affine term structure models.

2.2. Discrete-time arbitrage-free Nelson-Siegel model

The purpose of this section is to define a discrete-time model analogous to model (2.2)–(2.3). In what follows, an arbitrage-free market with discrete time steps $t=0,1,...,T$ , each separated by a time elapse $\Delta$ , is considered. Market dynamics are defined on probability space $(\Omega, \mathcal{F}_T,\mathbb{P})$ where $\mathbb{P}$ is the real-world probability measure and $\mathcal{F}\,:\!=\,\left\lbrace \mathcal{F}_t \right\rbrace_{t=0}^{T}$ is a filtration characterizing the information flow in the market.

In the discrete-time context, the time-t short risk-free rate $r_t$ is the $\mathcal{F}_t$ -measurable rate effective for the period $[t,t+1)$ . The short rate dynamics considered are

(2.7) \begin{align}r_t = X^{(1)}_t + X^{(2)}_t,\end{align}

(2.8) \begin{align}X_{t+1} = X_t+\kappa^{\mathbb{P}} \big(\theta^{\mathbb{P}}-X_t\big)+\Sigma Z_{t+1}^{\mathbb{P}}, \end{align}

where the stochastic factor vectors are $X_{t} =\left[X^{(1)}_t, X^{(2)}_t, X^{(3)}_t\right]^\top$ , $t=0,\ldots,T$ , $\kappa^{\mathbb{P}}$ and $\Sigma$ are constant parameter matrices of dimension $3\times 3$ , and $\theta^{\mathbb{P}}$ is a constant column vector of dimension 3. The process $Z^{\mathbb{P}}= \big\{ Z_{t}^{\mathbb{P}} \big\}^T_{t=1}$ with $ Z_{t}^{\mathbb{P}}= \left[Z^{\mathbb{P}}_{t,1}, Z^{\mathbb{P}}_{t,2}, Z^{\mathbb{P}}_{t,3}\right]^\top$ is a multivariate standard Gaussian white noise under $\mathbb{P}$ with contemporaneous correlation parameters $\rho^{(Z)}_{i,j} \equiv \text{corr}\!\left(Z_{t,i}^{\mathbb{P}},Z_{t,j}^{\mathbb{P}}\right)$ , $t=1,\ldots,T$ and $i,j=1,2,3$ . The initial value of the factor process is fixed at $X_0= x_0\in \mathbb{R}^3$ . $\Sigma$ is assumed to be a diagonal matrix with strictly positive elements. Equations (2.7)–(2.8) are meant to mimic (2.2)–(2.3), except that the model is specified under the physical measure $\mathbb{P}$ instead of under the risk-neutral measure $\mathbb{Q}$ . A more general class of Gaussian discrete-time affine structure models is studied in Wüthrich & Merz (Reference Wüthrich and Merz2013), but without specializing to the DTAFNS model which is the object of this paper.

To obtain the specification of spot rates under the DTAFNS model, risk-neutral dynamics need to be specified. Indeed, spot rates are given by

\begin{equation*}y(t,T) = -\frac{\log P(t,T)}{(T-t)\Delta}\end{equation*}

where P(t,T) is the time-t of a risk-free zero coupon of unit face value with maturity on time point T. In discrete-time models of the short rate, the relationship between P(t,T) and the risk-free rate is

\begin{equation*}P(t,T) = \mathbb{E}^{\mathbb{Q}}\!\left[ \exp\!\left(-\Delta \sum^{T-1}_{j=t} r_{j}\right) \bigg\vert \mathcal{F}_t\right].\end{equation*}

To obtain risk-neutral dynamics, following the lines of Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) for instance, a discrete-time version of the Girsanov theorem can be invoked to apply a translation to the drift of short rate factors $X_t$ without altering their volatilities, while keeping a multivariate Gaussian distribution for innovations.Footnote 1 It follows from this theorem that for a given constant market price of risk matrix $\gamma \in \mathbb{R}^{3\times 3}$ , there exists a probability measure $\mathbb{Q}$ equivalent to $\mathbb{P}$ such that the process $ Z^{\mathbb{Q}}=\big\lbrace Z^{ \mathbb{Q} }_t\big\rbrace_{t=1}^{T} $ defined through $Z^{\mathbb{Q}}_{t+1}=Z^{\mathbb{P}}_{t+1}-\gamma X_t $ is also an $\mathcal{F}$ -adapted standard Gaussian white noise under $\mathbb{Q}$ , still with contemporaneous correlation matrix $\rho^{(Z)} = \left[ \rho^{(Z)}_{i,j} \right]^3_{i,j=1}$ . Substituting $ Z^{\mathbb{P}}_{t+1}=Z^{\mathbb{Q}}_{t+1}+\gamma X_t $ into (2.8) leads to

\begin{equation*} X_{t+1}=X_t+\kappa^{\mathbb{P}} \theta^{\mathbb{P}} - \left(\kappa^{\mathbb{P}} - \Sigma \gamma\right) X_t +\Sigma Z_{t+1}^{\mathbb{Q}}.\end{equation*}

Assuming there is a solution $(\kappa^{\mathbb{Q}},\theta^{\mathbb{Q}})$ with $\kappa^{\mathbb{Q}} \in \mathbb{R}^{3\times 3}$ and $\theta^{\mathbb{Q}} \in \mathbb{R}^{3\times 1}$ to the system of equations

(2.9) \begin{eqnarray}\begin{cases}\kappa^{\mathbb{Q}} \theta^{\mathbb{Q}} = \kappa^{\mathbb{P}} \theta^{\mathbb{P}},\\[4pt] \kappa^{\mathbb{Q}} = \kappa^{\mathbb{P}} - \Sigma \gamma,\end{cases}\end{eqnarray}

then

(2.10) \begin{align}X_{t+1} = X_t+\kappa^{\mathbb{Q}} (\theta^{\mathbb{Q}}-X_t)+\Sigma Z_{t+1}^{\mathbb{Q}} \end{align}

and thus, the factor process $\{X_t\}^T_{t=0}$ preserves its auto-regressive structure under the risk-neutral measure. Note that as a future research avenue, other pricing kernels can be considered. For instance, these associated with essentially affine models defined in Duffee (Reference Duffee2002) are shown in that work to enhance the forecasting ability of the model by generalizing the physical dynamics of term structure factors while retaining the same risk-neutral dynamics.

The DTAFNS model relies on specific choices of parameters mimicking the AFNS model. Zero-coupon bond prices under such specific choices are hereby derived. The proofs for the next results are provided in Appendix A.

Lemma 2.1. For some integer $\tau>0$ and a real number $r \neq 1$ ,

(2.11) \begin{align}\zeta_0(r,\tau) \equiv \sum_{u=1}^{\tau-1}r^{u} =\dfrac{r-r^{\tau}}{1-r},\end{align}

(2.12) \begin{align}\zeta_1(r,\tau) \equiv \sum_{u=1}^{\tau-1}u r^{u} = \dfrac{r- \tau r^{\tau}+(\tau-1) r^{\tau+1}}{(1-r)^2},\end{align}

(2.13) \begin{align}\zeta_2(r,\tau) \equiv \sum_{u=1}^{\tau-1}u^2 r^{u} = \dfrac{ -(\tau-1)^2 r^{\tau+2} + (2\tau^2-2\tau-1)r^{\tau+1} - \tau^2 r^{\tau} + r^2+r}{(1-r)^3}.\end{align}

Proposition 2.1. Consider $\lambda\in (0,1)$ and assume

(2.14) \begin{equation} \kappa^{\mathbb{Q}} = \left( {\begin{array}{c@{\quad}c@{\quad}c} 0 & 0 & 0 \\[4pt] 0 & \lambda & -\lambda \\[4pt] 0& 0 & \lambda \\[4pt] \end{array}} \right). \end{equation}

In other words, suppose that the state variables $X_t$ have the following risk-neutral dynamics:

\begin{align*} \left( {\begin{array}{c} X^{(1)}_{t+1}-X^{(1)}_t \\[4pt] X^{(2)}_{t+1}-X^{(2)}_t \\[4pt] X^{(3)}_{t+1}-X^{(3)}_t \\[4pt] \end{array} } \right) = \left( {\begin{array}{c@{\quad}c@{\quad}c} 0 & 0 & 0 \\[4pt] 0 & \lambda & -\lambda \\[4pt] 0& 0 & \lambda \\[4pt] \end{array}} \right) \left[ {\begin{array}{c} \theta^{\mathbb{Q}}_1- X^{(1)}_t \\[4pt] \theta^{\mathbb{Q}}_2- X^{(2)}_t \\[4pt] \theta^{\mathbb{Q}}_3 - X^{(3)}_t \end{array} } \right] + \left( {\begin{array}{c@{\quad}c@{\quad}c} \Sigma_{1,1} & 0 & 0 \\[4pt] 0 & \Sigma_{2,2} & 0 \\[4pt] 0& 0 & \Sigma_{3,3} \\[4pt] \end{array} } \right) \left( {\begin{array}{c} Z^{\mathbb{Q}}_{t+1,1} \\[4pt] Z^{\mathbb{Q}}_{t+1,2} \\[4pt] Z^{\mathbb{Q}}_{t+1,3} \\[4pt] \end{array} } \right).\end{align*}

Then, the time-t arbitrage-free price of a zero-coupon maturing at T is, for $t=0,\ldots,T-1$ ,

(2.15) \begin{equation}P(t,T)= A_\tau\exp \big[{-}\Delta B_\tau^\top X_t\big], \end{equation}

where $\tau=T-t$ , $B_\tau=\left[B^{(1)}_\tau, \,\, B^{(2)}_\tau, \, \, B^{(3)}_\tau\right]^\top$ and

\begin{align*}B^{(1)}_\tau&= \tau, \quad B^{(2)}_\tau= \dfrac{1-(1-\lambda )^\tau}{\lambda }, \quad B^{(3)}_\tau= \frac{1-(1-\lambda )^{\tau-1}}{\lambda } - (\tau-1) (1-\lambda )^{\tau-1}, \\[4pt]\log A_\tau &= -\Delta\theta^{\mathbb{Q}}_2 \!\left(B^{(1)}_\tau - B^{(2)}_\tau\right) + \Delta \theta^{\mathbb{Q}}_3 B^{(3)}_\tau + \frac{1}{2} \Delta^2\upsilon_\tau\end{align*}

with

\begin{align*} \upsilon_\tau &= \sum^3_{i=1} \sum^3_{j=1}\upsilon^{(i,j)}_\tau\\[4pt]\upsilon^{(1,1)}_\tau &= \Sigma^2_{1,1}\dfrac{\tau(\tau-1)(2\tau-1)}{6},\\[4pt]\upsilon^{(2,2)}_\tau &= \frac{\Sigma^2_{2,2}}{\lambda ^2} \left( \tau - 2 \!\left[\frac{1-(1-\lambda )^\tau }{\lambda }\right] + \frac{1-(1-\lambda )^{2 \tau}}{1-(1-\lambda )^2}\right),\\[4pt]\upsilon^{(3,3)}_\tau &= \unicode{x1D7D9}_{ \{ \tau > 1\} } \frac{\Sigma^2_{3,3}}{\lambda ^2} \Bigg[ \tau-2 + \zeta_0\!\left( (1-\lambda )^2,\tau-1\right)+\lambda ^2 \zeta_2\!\left( (1-\lambda )^2,\tau-1\right)\\[4pt]& \quad - 2\zeta_0\!\left( (1-\lambda ),\tau-1\right) - 2\lambda \zeta_1\!\left( (1-\lambda ),\tau-1\right) + 2\lambda \zeta_1\!\left( (1-\lambda )^2,\tau-1\right) \bigg],\\[4pt]\upsilon^{(1,2)}_\tau &= \upsilon^{(2,1)}_\tau = \rho_{1,2}\Sigma_{1,1} \Sigma_{2,2} \frac{1}{\lambda }\left( \frac{\tau(\tau-1)}{2} - \zeta_1 ((1-\lambda ),\tau)\right),\\[4pt]\upsilon^{(1,3)}_\tau &= \upsilon^{(3,1)}_\tau = \unicode{x1D7D9}_{ \{ \tau > 1\} } \rho_{1,3}\Sigma_{1,1} \Sigma_{3,3} \frac{1}{\lambda } \bigg[\frac{\tau(\tau-1)}{2}-1 -\zeta_0\!\left( (1-\lambda ),\tau-1\right) \\[4pt] & \quad -(\lambda +1)\zeta_1\!\left( (1-\lambda ),\tau-1\right)-\lambda \zeta_2\!\left( (1-\lambda ),\tau-1\right) \bigg],\\[4pt]\upsilon^{(2,3)}_\tau &= \upsilon^{(3,2)}_\tau = \unicode{x1D7D9}_{ \{ \tau > 1\} } \rho_{2,3}\Sigma_{2,2} \Sigma_{3,3} \bigg( \frac{\tau-2- (2-\lambda )\zeta_0\!\left( (1-\lambda ),\tau-1\right) + (1-\lambda )\zeta_0\!\left( (1-\lambda )^2,\tau-1\right) }{\lambda ^2} \\[4pt] &\quad \quad \quad \quad \quad \quad + \frac{- \zeta_1\!\left( (1-\lambda ),\tau-1\right) + (1-\lambda )\zeta_1\!\left( (1-\lambda )^2,\tau-1\right) }{\lambda } \bigg).\end{align*}

This entails that spot rates have the form

(2.16) \begin{equation} y(t,T) = X_t^{(1)} + \dfrac{ (1-(1-\lambda )^\tau)}{\tau \lambda }X_t^{(2)} + \bigg [ \dfrac{ (1- (1-\lambda )^{\tau-1} )}{\tau \lambda } - \dfrac{\tau-1}{\tau} (1-\lambda )^{\tau-1} \bigg ] X_t^{(3)} - \dfrac{1}{\Delta\tau} \log A_\tau.\end{equation}

To obtain risk-neutral dynamics of the form proposed in Proposition 2.1, the following physical dynamics is considered:

\begin{equation*} \kappa^{\mathbb{P}} = \left( {\begin{array}{c@{\quad}c@{\quad}c} \kappa^{\mathbb{P}}_{1,1} & 0 & 0 \\[4pt] 0 & \kappa^{\mathbb{P}}_{2,2} & -\lambda \\[4pt] 0& 0 & \kappa^{\mathbb{P}}_{3,3} \\[4pt] \end{array}} \right), \quad \theta^{\mathbb{P}} = \left( {\begin{array}{c} 0 \\[4pt] \theta^{\mathbb{P}}_2 \\[4pt] \theta_3^{\mathbb{P}} \end{array} } \right)\end{equation*}

which allows meeting the second condition in (2.9) with

(2.17) \begin{eqnarray}\gamma &=& \left( {\begin{array}{c@{\quad}c@{\quad}c} \gamma_1 & 0 & 0 \\[4pt] 0 & \gamma_2 & 0 \\[4pt] 0& 0 & \gamma_3 \\[4pt] \end{array}} \right), \quad \gamma_i = \left( \kappa^{\mathbb{P}}_{i,i} -\unicode{x1D7D9}_{ \{i>1\} } \lambda \right)/\Sigma_{i,i}, \, i=1,2,3.\end{eqnarray}

Furthermore, substituting (2.17) in the first equation of (2.9) leads to

\begin{equation*} \theta_2^{\mathbb{Q}} = \lambda^{-1}\!\left(\kappa_{2,2}^{\mathbb{P}} \theta_2^{\mathbb{P}} + \kappa_{3,3}^{\mathbb{P}} \theta_3^{\mathbb{P}} - \lambda \theta_3^{\mathbb{P}}\right), \quad \theta_3^{\mathbb{Q}} = \lambda^{-1} \kappa_{3,3}^{\mathbb{P}} \theta_3^{\mathbb{P}}.\end{equation*}

Conversely, physical parameters can be retrieved from risk-neutral parameters through

(2.18) \begin{equation} \kappa^{\mathbb{P}}_{i,i} = \unicode{x1D7D9}_{ \{i>1\} } \lambda+\Sigma_{i,i}\gamma_i, \, i=1,2,3,\quad \theta_3^{\mathbb{P}}=\lambda \theta_3^{\mathbb{Q}}/\kappa^{\mathbb{P}}_{3,3}, \quad \theta_2^{\mathbb{P}}=\frac{\lambda}{\kappa^{\mathbb{P}}_{2,2}} \left[\theta_2^{\mathbb{Q}} -\frac{\theta_3^{\mathbb{Q}}}{\kappa^{\mathbb{P}}_{3,3}}\left(\kappa^{\mathbb{P}}_{3,3}-\lambda\right) \right].\end{equation}

3.1. Estimation framework

Several methods have been considered in the literature to estimate factor models of interest rates, such as conventional maximum likelihood estimation (Chen & Scott, Reference Chen and Scott1993), maximum likelihood with Kalman filters (Duan & Simonato, Reference Duan and Simonato1999; Lemke, Reference Lemke2006; Park, Reference Park2014; Augustyniak et al., Reference Augustyniak, Godin and Hamel2021), the method of simulated moments (Dai & Singleton, Reference Dai and Singleton2000), or Bayesian approaches (Hong et al., Reference Hong, Niu and Zeng2019). In this paper, the Kalman filter approach is applied. Observed spot rates (observable data) are considered to be a noisy version of the true spot rates (the signal). Setting up the Kalman filter requires to first derive the recursive Gaussian linear relationship between the observable quantities and factors driving the signal, which is hereby provided.

Assume that on each time t a set of M annualized continuously compounded spot rates with fixed times-to-maturity $n_1, ..., n_M$ are observed, which are denoted by the column vector

\begin{equation*}\hat{\textsf{y}}(t) \equiv (\hat{y}(t,t+n_1), ..., \hat{y}(t,t+n_M))^\top.\end{equation*}

Define also the corresponding model-implied spot rate vector

\begin{equation*} \textsf{y}(t) \equiv \left(y(t,t+n_1),...,y(t,t+n_M)\right)^\top\end{equation*}

where (2.16) leads to

(3.1) \begin{align}y(t,t+n) = \dfrac{-1}{n\Delta}\log A_n+\dfrac{1}{n}B_n^\top X_t.\end{align}

The relationship between observed and model-implied spot rate is assumed to be

\begin{equation*} \hat{\textsf{y}}(t)=\textsf{y}(t)+\eta_{\textsf{t}}\end{equation*}

where $\{\eta_{\textsf{t}}\}^T_{t=1}$ is a multivariate Gaussian white noise with diagonal variance-covariance matrix H. Furthermore, for simplicity, all values on the diagonal of matrix H are assumed to be identical and equal to some parameter h. This leads to

(3.2) \begin{align}\hat{\textsf{y}}(t) &=\textsf{a}+\textsf{B}X_t+\eta_{\textsf{t}}, \end{align}

where

(3.3) \begin{align}\textsf{a} & \equiv\left(\dfrac{-1}{n_1\Delta} \log A_{n_1},\dots,\dfrac{-1}{n_M\Delta} \log A_{n_M}\right)^\top,\\\beta_{n} & \equiv\left(\dfrac{B^{(1)}_n}{n},\dfrac{B^{(2)}_n}{n},\dfrac{B^{(3)}_n}{n}\right)^\top, \nonumber\end{align}

(3.4) \begin{align}\textsf{B} \equiv \left({\beta}_{n_1},\dots, {\beta}_{n_M}\right)^\top .\qquad\qquad\qquad\qquad\end{align}

Furthermore, the transition equation describing the dynamics of the latent factors $X_t$ is obtained from (2.8):

(3.5) \begin{equation}X_{t+1}=\textbf{b}+\textbf{D}X_t+{\xi}_{t+1},\end{equation}

where

(3.6) \begin{align}\textbf{b}\equiv\kappa^{\mathbb{P}} \theta^{\mathbb{P}},\quad \textbf{D}\equiv I-\kappa^{\mathbb{P}}, \quad {\xi}_{t+1}&\equiv \Sigma Z_{t+1}^{\mathbb{P}}.\end{align}

The sequence $\{\xi_{\textsf{t}}\}^T_{t=1}$ is therefore a multivariate Gaussian white noise with covariance matrix

(3.7) \begin{equation} \textbf{Q} =\Sigma \rho \Sigma.\end{equation}

The representation (3.2)–(3.5) of state space variables allows estimating the DTAFNS model using a Kalman filter. The general Kalman filter algorithm is presented in Appendix D of the Online Appendix. The Kalman filter allows deriving the log-likelihood of observed spot rate curves for a candidate set of parameters. Algorithm 1 indicates steps to calculate the log-likelihood function in the context of the DTAFNS model. In this algorithm, p denotes joint density functions. In the present work, the implementation of the Kalman filter is performed with the R package FKF. Moreover, to find maximum likelihood estimates of the parameters, consistently with Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021), an optimization of the likelihood is conducted with the R package Rsolnp (Ghalanos & Theussl, Reference Ghalanos and Theussl2015) applying the nonlinear augmented Lagrange multiplier optimization method of Ye (Reference Ye1987). The optimization is performed under the following constraints:

\begin{equation*} \Sigma_{i,i}\geq 0, \quad h \geq 0, \quad -1 \leq \rho_{i,j} \leq 1, \, i,j=1,2,3, \quad \lambda \in (0,1).\end{equation*}

Additional constraints could be applied to enforce the positive semi-definiteness of matrices $\Sigma$ and $\rho$ , but optimization results obtained with our dataset satisfy such constraints without explicitly specifying them to the optimizer.

Algorithm 1 Kalman filter algorithm for the calculation of likelihood function and smoothed state densities

The smoothed state density is the density of time-t factors conditional on all observations from the sample, namely $p(X_{t}|\hat{y}(1),\dots,\hat{y}(T))$ . Such distribution is multivariate Gaussian with mean $\bar{x}_{t|T}$ and covariance matrix $P_{t|T}$ . Such mean and variance are obtained with the Kalman smoother algorithm (Shumway & Stoffer, Reference Shumway and Stoffer2017) which is also shown in Algorithm 1.

3.2. Model estimation and analysis

This section presents the DTAFNS model estimation using historical Canadian spot rate data.

3.2.2. Estimated parameters and performance

Table 1 gives estimates of DTAFNS model parameters obtained with Algorithm 1 which applies the Kalman filter. Factor 1 innovations are negatively correlated with these of the two other factors, which are positively correlated between themselves. The negative correlation between the level factor $X^{(1)}$ and the slope and curvature factors $X^{(2)}$ and $X^{(3)}$ tends to reduce the frequency of very large or low values of the short rate $r_t$ . Furthermore, it implies that when the long-term rates increase through higher values of $X^{(1)}$ , the slope given by $-X^{(2)}$ also tends to go up. Thus, movements in the long-end of the curve are not necessarily matched by movements of a similar magnitude and direction in the short-end. The hump in the curve reflected by $X^{(3)}$ also tends to decrease during periods of increasing of long-term spot rates. The speed of reversion $\kappa^{\mathbb{P}}_{i,i}$ is much lower for the first factor ( $i=1$ ) than for the other two, which is not surprising due to the overall rates level having had a clear lasting downward trend in the historical sample, representing close-to-non-stationary behavior. Finally, the estimation implies a long-term average of short rates of $\theta^{\mathbb{P}}_{2}= 0.0301$ under physical dynamics.

Table 1. Maximum likelihood estimates of the DTAFNS model parameters

$\rho$ = $\displaystyle\left[\begin{array}{c@{\quad}c@{\quad}c} 1 & -0.6303 & -0.4097 \\[4pt] -0.6303 & 1 & 0.2993 \\[4pt] -0.4097 & 0.2993 & 1 \\[4pt] \end{array}\right] $ , $P_{1|0}$ = $\left[\begin{array}{c@{\quad}c@{\quad}c} 4.45\!\times\! 10^{-6} & 0 & 0 \\[4pt] 0 & 4.45\!\times\! 10^{-6} & 0 \\[4pt] 0 & 0 & 4.45\!\times\! 10^{-6} \\[4pt] \end{array}\right] $

Notes: Maximum likelihood parameter estimates for the DTAFNS model presented in Section 2.2 obtained with the Kalman filter (see Algorithm 1). The data sample includes Canadian end-of-month yield curves from January 1986 to January 2022. $\bar{x}_{i,1|0}$ refers to element i of the vector $\bar{x}_{1|0}$ .

Fig. 1 shows the smoothed value (the smoothed distribution expected value) of latent factors provided by the Kalman smoother procedure described in Algorithm 1. The black curve is the short rate implied by the model which is obtained by summing the smoothed values of the first two factors, consistently with (2.7). The decreasing trend of the short rate throughout the data sample is recovered by the downward trend of Factor 1. Conversely, Factor 2 and Factor 3, driving the slope and humps in the curve respectively, exhibit more stationary fluctuations.

Figure 1 Model-implied factors and short rate time series. Notes: Time series of DTAFNS model-implied factors which correspond to the smoothed state inferences $ E^{ \mathbb{P}}\!\left[ x_t^{(i)}|\hat{y}(1),\hat{y}(2),..., \hat{y}(T)\right]$ , $i =1, 2, 3$ provided by Algorithm 1, and implied short rates obtained by summing smoothed values of the first two factors. The model is estimated on the end-of-month Canadian spot rate curves extending from January 1986 to January 2022.

To assess the ability of the DTAFNS model to match observed short- and long-term spot rates of the term structure, time series of observed spot rates of three-month and 10-year tenors are compared to the model-implied spot rates for the same tenors. Such time series are provided in Fig. 2. The model can be seen to adequately match short- and long-term rates throughout the entire sample, as differences between the model-implied and observed spot rates are very small.

Figure 2 Model-implied and observed spot rate time series for 3-month and 10-year tenors. Notes: Time series of observed 3-month (short-term) and 10-year (long-term) maturity spot rates (dotted curves), with corresponding spot rates implied by the fitted DTAFNS model. The dataset considered is the end-of-month Canadian spot rate curves extending from January 1986 to January 2022.

3.3. Performance assessment and benchmarking

This section analyzes the adequacy of the fit and the predictive performance of the DTAFNS model through benchmarking.

3.3.1. Benchmarks

The discrete-time Gaussian three-factor model of Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021), subsequently referred to as the DG3 model, and a version of the dynamic Nelson-Siegel (DNS) model introduced by Diebold & Li (Reference Diebold and Li2006) are benchmarks considered in this study for comparison of performance results. See Appendix B for a formal specification of these two models. The main differences between DG3 and the DTAFNS model are that (i) the latent factors in the former only regress on themselves instead of allowing cross-factor auto-regressions, (ii) the DTAFNS has a short rate equal to the sum of two of the factors instead of the three as in the benchmark, and (iii) the DTAFNS has a fixed structure for the matrix $\kappa^{\mathbb{Q}}$ driving the auto-regression intensity, unlike the Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) model allowing optimizing such quantities. Regarding the second benchmark, unlike in the DTAFNS model, the DNS model does not include a deterministic adjustment term in spot rates to rule out arbitrage opportunities. Additionally, a nested version of the DTAFNS model which sets all correlation parameters to zero, i.e. $\rho_{i,j}=0$ , $i,j=1,2,3$ , is also considered for benchmarking purposes. This benchmark, which is denoted by the acronym DTAFNS-U, is used to evaluate the impact of the innovation correlations on performance. Table 2 shows parameter estimates for the DTAFNS-U model, which are slightly different but qualitatively similar to these of the DTAFNS model presented in Table 1.

Table 2. Maximum likelihood estimates of the DTAFNS-U model parameters

$\rho$ = $\displaystyle\left[\begin{array}{c@{\quad}c@{\quad}c} 1 & 0 & 0 \\[4pt] 0 & 1 & 0 \\[4pt] 0 & 0 & 1 \\[4pt] \end{array}\right] $ , $P_{1|0}$ = $\left[\begin{array}{c@{\quad}c@{\quad}c} 4\!\times\! 10^{-6} & 0 & 0 \\[4pt] 0 & 4\!\times\! 10^{-6} & 0 \\[4pt] 0 & 0 & 4\!\times\! 10^{-6} \\[4pt] \end{array}\right] $

Notes: Maximum likelihood parameter estimates for the DTAFNS-U model obtained with the Kalman filter. The DTAFNS-U is identical to the DTAFNS model except that the term structure innovations are uncorrelated, i.e. the matrix $\rho$ is set to the identity matrix. The data sample includes Canadian end-of-month yield curves from January 1986 to January 2022. $\bar{x}_{i,1|0}$ refers to element i of the vector $\bar{x}_{1|0}$ .

3.3.2. Goodness-of-fit of the model and forecasting ability

To assess the model’s ability to adequately fit the observed data, we compare the log-likelihood of the DTAFNS, the DTAFNS-U, the DG3 and the DNS models. Two types of log-likelihoods are computed: an in-sample version obtained by fitting the models to the entire dataset, and an out-of-sample version which is obtained by considering an expanding window sequential approach. For the latter, in the first iteration, data from January 1986 up to the end of 2016 is considered as a first training set, and data from 2017 is considered as a test set. In any following iteration, the training sets are expanded by one year and the test year moves to the year following that of the previous iteration. The last iteration has a test set including data for both 2021 and 2022 since the dataset only includes data for January in 2022. The aggregated out-of-sample log-likelihood is obtained by summing these associated with all test sets. Table 3 provides the results. The best-performing model in terms of aggregate log-likelihood over the full sample is the DTAFNS-U. However, when looking at year-by-year performance, the DTAFNS outperforms the DTAFNS-U for the three first years, while only underperforming the latter for the last year. Moreover, the DTAFNS beats both the DG3 and DNS benchmarks out-of-sample for all years except for the DG3 in years 2020 and 2021–2022. While the in-sample log-likelihood of the DG3 is higher than that of the three other competing models, such better fit does not translate to better (out-of-sample) predictive ability. These tests imply that either the DTAFNS or the DTAFNS-U should be favored over the DG3 and DNS benchmarks.

Table 3. Log-likelihood of the DTAFNS model and its benchmarks

Notes: Comparison of the log-likelihood of the DTAFNS model of Section 2.2 and of the DG3 and DNS benchmarks described in Appendix B. The data sample is the Canadian end-of-month yield curve. The in-sample dataset starts in January 1986 and ends in January 2022. The out-of-sample estimation procedure uses an expanding window approach described in Section 3.3.2. The aggregated out-of-sample log-likelihood is ultimately obtained by summing the log-likelihood for all test years.

To further analyze the predictive ability of the competing interest rate models, Appendix C provides descriptive statistics on out-of-sample point forecasts of spot rates for various tenors. The difference between the performance of the DTAFNS model and that of benchmarks seems overall marginal, with slightly higher (lower) performance exhibited by the DTAFNS for long (medium) tenors. Therefore, our analyses indicate that point forecasting performance of the various models is somewhat similar, but the distributional predictive ability of the DTAFNS model is slightly superior.

To visualize the ability of the DTAFNS model and benchmarks to replicate observed yield curves, Fig. 3 presents the realized spot curves and model-implied counterparts for four selected dates. These four days are the same that are considered in Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) to display fitting performance, namely December 29, 2006, December 31, 2008, June 30, 2016 and October 31, 2018. Such dates exhibit several spot curve shapes: flat, upward sloping or humped. The figure shows that the DTAFNS and all three benchmarks are able to reproduce observed yield curves reasonably well.

Figure 3 Model-implied and observed yield curves. Notes: Realized and model-implied spot rate curves on the four following dates: December 29, 2006, December 31, 2008, June 30, 2016 and October 31, 2018. Dotted black line: observed spot rates. Full green line: DTAFNS model-implied curve. Dashed blue line: DG3 benchmark implied curve. Red dotted-dashed line: DNS benchmark implied curve. Pink dotted-dashed line: DTAFNS-U benchmark implied curve. The observed data are end-of-month Canadian spot rates provided by the Bank of Canada.

4.1. Rolling bond fund returns

The dynamics of the mixed fund model of Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) involves a regression on movements of interest rate factors, among other quantities. The inclusion of such regressors is justified by an analysis of returns of a rolling bond fund, which are shown to depend on such drivers. Therefore, to adapt the Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021), the rolling bond fund dynamics is derived under the DTAFNS dynamics.

Consider a fund containing a single zero-coupon of time-to-maturity $\tau$ that is rolled-over on each period. Denoting its time- $t+1$ log-return by $R^{(\tau)}_{t+1}$ , (2.15) leads to

(4.1) \begin{align} R^{(\tau)}_{t+1} &=\log \!\left(\dfrac{P(t+1,t+\tau)}{P(t,t+\tau)}\right) \nonumber\\[4pt] &= \log A_{\tau-1} - \log A_{\tau} - \Delta B_{\tau-1}^\top X_{t+1} +\Delta B_\tau^\top X_t\nonumber\\[4pt] &=\log A_{\tau-1} - \log A_{\tau} - \Delta \sum_{i=1}^3 \left( B_{\tau-1}^{(i)} X_{t+1}^{(i)}-B_{\tau}^{(i)}X_{t}^{(i)}\right).\end{align}

Moreover, recalling that $r_t = X_t^{(1)}+X_t^{(2)}$ and using (2.14), (2.18), Lemma A.2 and (4.1), the excess return can be expressed as

(4.2) \begin{align} R^{(\tau)}_{t+1} - \Delta r_t &= \log A_{\tau-1} - \log A_{\tau} - \Delta \sum_{i=1}^3 \left( B_{\tau-1}^{(i)} X_{t+1}^{(i)}-(B_{\tau}^{(i)}-1) X_t^{(i)}\right) \nonumber\\& \quad +\Delta \sum_{i=1}^3 X_t^{(i)}- \Delta \!\left(X_t^{(1)}+X_t^{(2)}\right) \nonumber \\[4pt] &= \log A_{\tau-1} - \log A_{\tau} - \Delta \sum_{i=1}^3 \left(B_{\tau-1}^{(i)} X_{t+1}^{(i)}-\left(B_{\tau-1}^{(i)}\left(1-\kappa^{\mathbb{Q}}_{i,i}\right)\right.\right.\nonumber\\& \quad\left.\left.-\unicode{x1D7D9}_{ \{i=3\} } (1-\lambda)^{\tau-1}\right)X_t^{(i)}\right) + \Delta X_t^{(3)}\nonumber \\[4pt] &=\log A_{\tau-1} - \log A_{\tau} - \Delta\! \sum_{i=1}^3 B_{\tau-1}^{(i)}\!\left(X_{t+1}^{(i)}-\left(1-\kappa^{\mathbb{Q}}_{i,i}\right)X_t^{(i)}\right) + \Delta \!\left(1-(1-\lambda)^{\tau - 1}\right)X_t^{(3)}.\end{align}

Mimicking Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) by considering a portfolio consisting of M positions with weights $\omega_1, \dots, \omega_M$ (where $\sum_{j=1}^M\omega_j=1$ ) that are rolled-over on bonds with fixed time-to-maturities $\tau_1, \dots, \tau_M$ , the log-return of the rolling bond fund between times t and $t+1$ , denoted by $R_{t+1}^{(B)}$ , is approximately characterized by

(4.3) \begin{align} R^{(B)}_{t+1} - \Delta r_t & \approx \omega_1 R^{(\tau_1)}_{t+1}+\dots+ \omega_M R^{(\tau_M)}_{t+1}-\Delta r_t\nonumber\\[4pt] &=\omega_1\!\left( R^{(\tau_1)}_{t+1}-\Delta r_t\right)+\dots+ \omega_M \!\left(R^{(\tau_M)}_{t+1}-\Delta r_t\right)\nonumber\\[4pt] &=\sum_{j=1}^M\omega_j\bigg(\log A_{\tau_j-1} - \log A_{\tau_j} - \Delta \sum_{i=1}^3 B_{\tau_j-1}^{(i)}\left(X_{t+1}^{(i)}-\left(1-\kappa^{\mathbb{Q}}_{i,i}\right)X_t^{(i)}\right)\nonumber \\[4pt] &\quad + \Delta \!\left(1-(1-\lambda)^{\tau_j - 1}\right)X_t^{(3)}\bigg) \nonumber \\[4pt] &=\sum_{j=1}^M\omega_j\!\left(\log A_{\tau_j-1} - \log A_{\tau_j}\right) - \Delta \sum_{i=1}^3\sum_{j=1}^M\omega_j B_{\tau_j-1}^{(i)}\!\left(X_{t+1}^{(i)}-\left(1-\kappa^{\mathbb{Q}}_{i,i}\right)X_t^{(i)}\right) \nonumber \\[4pt] &\quad + \Delta \sum_{j=1}^M \omega_j\!\left(1-(1-\lambda)^{\tau_j - 1}\right)X_t^{(3)}.\end{align}

Relationship (4.3) will serve as the basis to model the fixed income part of the mixed fund in the proposed model. The key difference with the Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) model is the presence of the last term proportional to $X_t^{(3)}$ in (4.3), which is not present in their work.

4.2. Risky assets returns

The adaptation of the mixed fund model proposed in this paper keeps the same model for excess equity returns as in Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021). Returns of q equity assets are represented with the exponential generalized auto-regressive conditional heteroskedastic (EGARCH) model of Nelson (1991):

(4.4) \begin{align} R_{t+1,j}^{(S)}-\Delta r_t &= \lambda_j^{(S)}\sqrt{h_{t,j}^{(S)}}-\dfrac{1}{2}h_{t,j}^{(S)} + \sqrt{h_{t,j}^{(S)}}Z_{t+1,j}^{(S)},\end{align}

(4.5) \begin{align} \log h_{t,j}^{(S)}=\omega^{(S)}_j+\alpha_j^{(S)}Z_{t,j}^{(S)}+\gamma^{(S)}_j\bigg (|Z_{t,j}^{(S)}|-\dfrac{2}{\sqrt{2\pi}}\bigg )+\beta_j^{(S)}\log h_{t-1,j}^{(S)}, \end{align}

where $\big(\omega_j^{(S)}, \alpha_j^{(S)}, \gamma_j^{(S)}, \beta_j^{(S)}\big)$ are parameters associated with the conditional volatility. The process $\big\{h_{t,j}^{(S)}\big\}_{t=1}^{T-1}$ is $\mathcal{F}$ -adapted. Parameters $\lambda_j^{(S)}$ represent equity risk premia. For $j=1,\ldots,q$ , $Z_{j}^{(S)} \equiv \big\{Z_{t,j}^{(S)}\big\}_{t=1}^{T}$ are standard Gaussian white noises under $\mathbb{P}$ , independent of $Z_{t}$ , where the correlation between $Z_{t,i}^{(S)}$ and $Z_{t,j}^{(S)}$ is denoted $\Gamma_{i,j}$ .

As in Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021), the following risk-neutral dynamics for the q equity assets is assumed:

(4.6) \begin{align} R_{t+1,j}^{(S)}-\Delta r_t = -\dfrac{1}{2}h_{t,j}^{(S)} + \sqrt{h_{t,j}^{(S)}}Z_{t+1,j}^{\mathbb{Q}(S)}, \end{align}

(4.7) \begin{align} \log h_{t,j}^{(S)} =\omega^{(S)}_j+\alpha_j^{(S)}\!\left(Z_{t,j}^{\mathbb{Q}(S)}-\lambda_j^{(S)}\right)+\gamma^{(S)}_j\bigg (|Z_{t,j}^{\mathbb{Q}(S)}-\lambda_j^{(S)}|-\dfrac{2}{\sqrt{2\pi}}\bigg )+\beta_j^{(S)} \log h_{t-1,j}^{(S)}, \end{align}

where, under the risk-neutral measure $\mathbb{Q}$ , $\left\{Z_{t,j}^{\mathbb{Q}(S)} \right\}^T_{t=1} = \left\{Z_{t,j}^{(S)} + \lambda_j^{(S)} \right\}^T_{t=1}$ are also standard Gaussian white noises for each $j=1,\dots,q$ , still with a contemporaneous correlation matrix $\Gamma =[ \Gamma_{i,j}]_{i,j=1}^q$ . The independence between the equity innovations $\left\{Z_{t,j}^{\mathbb{Q}(S)} \right\}^T_{t=1}$ and interest rate innovations $\left\{Z_{t}^{\mathbb{Q}} \right\}^T_{t=1}$ is assumed to also hold under $\mathbb{Q}$ .

4.3. Mixed fund dynamics

To motivate the structure of the proposed mixed fund model, consider a portfolio invested in the rolling bond fund with weight W and in equity indices with respective weights $\tilde\psi^{(S)}_1,\ldots,\tilde\psi^{(S)}_q$ , where $W+ \sum^q_{j=1} \tilde\psi^{(S)}_j =1 $ . Its excess return would be represented as

(4.8) \begin{eqnarray}\tilde{R}_{t+1}^{(F)}- \Delta r_t &\approx& W\!\left(R_{t+1}^{(B)} - \Delta r_t\right) + \sum_{j=1}^q \tilde\psi^{(S)}_j \!\left(R_{t+1,j}^{(S)} - \Delta r_t\right)\nonumber\\[4pt] &=& \tilde{\psi}_0+ \sum_{i=1}^3 \tilde{\psi}_i\!\left(X_{t+1}^{(i)}-\left(1\!-\!\kappa^{\mathbb{Q}}_{i,i}\right)X_t^{(i)}\right)\!+ \tilde{\psi}^{\prime}_3 X_t^{(3)}+ \sum_{j=1}^q \tilde\psi^{(S)}_j \!\left(R_{t+1,j}^{(S)} \!-\! \Delta r_t\right)\end{eqnarray}

where

(4.9) \begin{align}\tilde{\psi}_0 = W \sum_{j=1}^M\omega_j\left(\log A_{\tau_j-1} - \log A_{\tau_j}\right), \end{align}

(4.10) \begin{align} \tilde{\psi}_i = - W \Delta \sum_{j=1}^M\omega_j B_{\tau_j-1}^{(i)},\quad i=1,2,3,\end{align}

(4.11) \begin{align} \tilde{\psi}^{\prime}_3 = W \Delta \sum_{j=1}^M \omega_j\!\left(1-(1-\lambda)^{\tau_j - 1}\right).\end{align}

Equation (4.8) serves as the conceptual basis for the proposed mixed fund model. However, instead of relying on (4.9)–(4.11), model parameters can be directly estimated from the data in an econometric fashion. Indeed, allocations to fixed income and equity assets from the mixed funds are not exactly identical to those implied by previous assumptions – that is an investment in a rolling horizon bond fund and in equity indices – and basis risk needs to be taken into account. Moreover, as in Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021), basis risk is further represented through a noise component also following EGARCH dynamics. Mixed fund excess returns are thus assumed to be of the form

(4.12) \begin{align} R_{t+1}^{(F)} - \Delta r_t = \psi_0+ \sum_{i=1}^3 \psi_i\!\left(X_{t+1}^{(i)}\!-\!\left(1\!-\!\kappa^{\mathbb{Q}}_{i,i}\right)X_t^{(i)}\right)\!+\! \psi^{\prime}_3 X_t^{(3)}\!+\! \sum_{j=1}^q \psi^{(S)}_j \left(R_{t+1,j}^{(S)} - \Delta r_t\right) \!+\!\sqrt{h^{(F)}_{t}}Z^{(F)}_{t+1}, \end{align}

(4.13) \begin{align} \log h_{t}^{(F)} =\omega^{(F)}+\alpha^{(F)}Z_{t}^{(F)}+\gamma^{(F)}\bigg (|Z_{t}^{(F)}|-\dfrac{2}{\sqrt{2\pi}}\bigg)+\beta^{(F)} \log h_{t-1}^{(F)},\end{align}

where the basis risk innovation process $Z^{(F)}=\big\{Z_t^{(F)}\big\}_{t=1}^T$ is a standard Gaussian white noise under $\mathbb{P}$ , independent of interest rate innovations $Z^{\mathbb{P}}$ and equity innovations $Z_{j}^{(S)}, j=1,\dots,q$ . Basis risk volatility parameters $\big(\omega^{(F)},\alpha^{(F)}, \gamma^{(F)}, \beta^{(F)}\big)$ and linear coefficients $\big(\psi_0,\psi_1,\psi_2, \psi_3,\psi^{\prime}_3,\psi^{(S)}_1,\ldots,\psi^{(S)}_q\big)$ are model parameters to be estimated.

The risk-neutral dynamics can also be obtained along the lines of Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) by translating basis risk innovations $\big\{Z_{t}^{(F)} \big\}^T_{t=1}$ by some adapted process $\big\{\lambda_{t}^{(F)} \big\}^{T-1}_{t=0}$ : under $\mathbb{Q}$ , innovations $\big\{Z_{t}^{\mathbb{Q}(F)} \big\}^T_{t=1} \equiv \big\{Z_{t}^{(F)} + \lambda^{(F)}_{t-1} \big\}^T_{t=1}$ form a Gaussian white noise independent of risk-neutral term structure and equity innovations $\big\{Z_{t}^{\mathbb{Q}} \big\}^T_{t=1}$ and $\big\{Z_{t,j}^{\mathbb{Q}(S)} \big\}^T_{t=1}$ , respectively. The determination of the process $\big\{\lambda_{t}^{(F)} \big\}^{T-1}_{t=0}$ stems from the martingale property that must be satisfied by the mixed fund value process. Its value, which characterizes the risk-neutral dynamics of the mixed fund, is provided in the following proposition whose proof is in Appendix A.

Proposition 4.1. The risk-neutral dynamics of mixed fund excess returns are given by

(4.14) \begin{align} R_{t+1}^{(F)}-\Delta r_t = -\dfrac{1}{2}\bigg(\sigma_t^{(F)}\bigg)^2+\sigma_t^{(F)}\epsilon^{\mathbb{Q}(F)}_{t+1} \end{align}

where

(4.15) \begin{align} \bigg(\sigma_t^{(F)}\bigg)^2 &\equiv \sum_{i=1}^3\sum_{l=1}^3\psi_i \psi_l \Sigma_{i,i} \Sigma_{l,l}\rho_{i,l}+\sum_{j=1}^{q}\sum_{k=1}^q\psi_j^{(S)}\psi_k^{(S)}\Gamma_{j,k}\sqrt{h_{t,j}^{(S)}h_{t,k}^{(S)}}+h_t^{(F)}, \\[4pt] \phi_t &\equiv \psi_0 + \sum_{i=1}^3 \psi_i \!\left( \kappa^{\mathbb{Q}}_{i,i} \theta^{\mathbb{Q}}_i -\lambda \theta^{\mathbb{Q}}_3 \unicode{x1D7D9}_{ \{ i=2\} } \right) + (\psi_2 \lambda +\psi^{\prime}_3) X_t^{(3)} - \sum_{j=1}^q \psi_j^{(S)} \left( \dfrac{1}{2}h_{t,j}^{(S)}\right), \nonumber \\[4pt] \epsilon^{\mathbb{Q}(F)}_{t+1}& \equiv \dfrac{\sum_{i=1}^p\psi_i\Sigma_{i,i} Z_{t+1}^{\mathbb{Q}(i)}+\sum_{j=1}^q\psi_j^{(S)}\sqrt{h_{t,j}^{(S)}}Z_{t+1,j}^{\mathbb{Q}(S)}+\sqrt{h_{t}^{(F)}}Z_{t+1}^{\mathbb{Q}(F)}}{\sigma_t^{(F)}}, \nonumber \\[4pt] \log h_{t}^{(F)}&=\omega^{(F)}+\alpha^{(F)}\!\left(Z_{t}^{\mathbb{Q}(F)}-\lambda_{t-1}^{(F)}\right)+\gamma^{(F)}\bigg (|Z_{t}^{\mathbb{Q}(F)}-\lambda_{t-1}^{(F)}|-\dfrac{2}{\sqrt{2\pi}}\bigg )+\beta^{(F)} \log h_{t-1}^{(F)}, \nonumber \\[4pt] \lambda_{t}^{(F)} \equiv& \dfrac{1}{\sqrt{h_t^{(F)}}}\Bigg[\phi_t+\dfrac{1}{2}\bigg(\sigma_t^{(F)}\bigg)^2\bigg]. \nonumber\end{align}

and the process $\left\{\epsilon^{\mathbb{Q}(F)}_{t+1}\right\}^T_{t=1}$ is a sequence of independent standardized Gaussian variables under $\mathbb{Q}$ .

4.4. Estimation of the mixed fund model

The mixed fund model developed above is estimated on real data. The equity model (4.4) and (4.5) is estimated for two equity assets ( $q=2$ ), namely the S&P/TSX Composite and S&P 500 stock indices. For the mixed fund model (4.12)–(4.13), we consider the Assumption/Louisbourg Balanced Fund A.Footnote 4 This is a mixed bond and equity fund composed approximately of 39% of Canadian fixed income, 36% of Canadian equity, 15% of US equity, and 10% of other products. Monthly price return data (NAV returns for the mixed fund) in Canadian currency from February 1986 (from February 1996 for the mixed fund) to January 2022 are considered, a span that matches the term structure data used in earlier sections. Note that for the mixed fund, 14 out of the 326 monthly returns are missing values.

Tables 5 and 6 respectively show maximum likelihood estimated parameter and corresponding standard errors for the equity model (4.4)–(4.5) and the mixed fund model (4.12)–(4.13).Footnote 5 Note that initial variances $h^{(S)}_{0,j}$ , $j=1,2$ and $h^{(F)}_0$ are also considered as parameters to optimize during the estimation. Equity indices are highly correlated with $\Gamma_{1,2}=0.753$ . Furthermore, the persistence of the volatility for both indices is high and similar, with $\beta_1^{(S)}=0.628$ for the S&P TSX Composite and $\beta_2^{(S)}=0.697$ the S&P 500. Furthermore, the presence of a leverage effect is implied by negative values of $\alpha_j^{(S)}$ , $j=1,2$ . For the mixed fund model, all parameter estimates are significant at the 99% confidence level. Negative values for estimates of coefficient $\psi_i$ , $i=1,2$ reflect negative correlation between variations of the short rate and returns of the mixed fund. This is expected due to the fund being partially invested in fixed income, where bond prices are inversely related to interest rates. Interestingly, the estimate of coefficient $\psi^{\prime}_3$ , a parameter not present in the Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) model, is also significant and negative. It implies that mixed fund returns are negatively related to the humpiness of the term structure (and not only to changes in humpiness).

Table 5. Bivariate EGARCH equity model parameter estimates

Notes: Maximum likelihood estimates and standard errors (in parentheses) of the equity model (4.4)–(4.5). Indices $j=1$ and $j=2$ are respectively the S&P/TSX Composite and the S&P 500. The estimation is performed on monthly price return time series extending from February 1986 to January 2022 (432 returns by index).

The performance of our mixed fund model (4.12)–(4.13) is compared to that of Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021), which basically imposes the constraint $\psi^{\prime}_3=0$ and uses the DG3 instead of the DTAFNS model for the computation of smoothed values for term structure factors X. Our model has an in-sample log-likelihood of $1038.719$ , whereas that of the Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021) benchmark is $1035.714$ . Such reported figures are conditional log-likelihood for the mixed fund returns specification given smoothed values of the term structure factors for each respective models and equity index returns. Since the use of our mixed fund model does not lead to a deterioration in fitting performance, and since such model is compatible with the DTAFNS model which has arguably more favorable properties than the DG3 counterpart, the use of our mixed fund model shall be preferred over that of Augustyniak et al. (Reference Augustyniak, Godin and Hamel2021).

Table 6. Mixed fund model parameter estimates

Notes: Maximum likelihood estimates and standard errors (in parentheses) of the mixed fund model (4.12)–(4.13). The estimation is performed on a monthly NAV return time series for the Assumption/Louisbourg Balanced Fund A extending from February 1996 to January 2022 (312 returns).