Affine Term Structure Models with Unspanned Stochastic Volatility in Illiquid Derivative Markets

1. Introduction

Interest rate volatility is increasingly recognised as a distinguished feature in ATSMs and also an independent driver of time-varying risk premia in fixed income markets [1]. Typically, volatility is unobservable and potentially unspanned by the prices of traded bonds, raising significant challenges for both modelling and empirical inference. In particular, the inability to observe volatility directly—especially in markets with sparse derivatives trades due to lack of liquidity—necessitates models where volatility is treated as a latent factor. Such models must be capable of distinguishing between spanned and unspanned components of risk, a subtle but crucial distinction for pricing, hedging, and risk management.

Unspanned Stochastic Volatility (USV) refers to volatility risk that influences the dynamics of yields and risk premia but cannot be perfectly hedged or inferred from the cross-section of bond prices alone; see [1]. It results from a tension that exists between fitting a cross-section of yields and matching the time series dynamics of bond prices; see [2,3,4] and [5] among others. This introduces fundamental incompleteness into the fixed income market, making it impossible to replicate all sources of risk through portfolios made up solely on bonds—suggesting the presence of portfolio volatility risk [2]. The presence of USV can explain persistent empirical puzzles such as the failure of (ATSMs) to match observed option-implied volatilities or the dynamic behaviour of the market price of risk (MPR) [1,6] and [7].

One of our central objectives is to evaluate the role of model restrictions in enabling or suppressing the presence of USV. In doing so, we aim to test affine models both with USV, leveraging off the maximal models from [8] and [3] which in turn, are characterised by a greater number of identifiable parameters. The representation that we follow leads to specifications that are more flexible than the canonical form identified by [8]. These are referred to as a maximal $A_M\left(N\right)$ form of ATSMs which in some cases have more identifiable parameters. Our focus is on the impact of parameter restrictions on these maximal ATSMs, the three-factor $A_1$(3) and four-factor $A_1\left(4\right)$ models [3]. Subsequently, we compare and contrast both the $A_1\left(3\right)$USV and $A_1\left(4\right)$USV in terms of their time series and yield curve dynamics.

As in most emerging markets, we anticipate that market data constraints are significant due to non-liquid derivative trading. In such markets, data on derivatives may be sparse and irregular. [5] points out that his approach to test the presence of USV does not rely on derivatives market data to estimate the underlying risk factors under the physical measure. We consider the yield curve itself to reveal information about latent volatility. While the cross-section of yields provides insight about the first moment, higher-order dynamics — such as yield curvature, and term premium variation, — changes in the covariance structure of yields may be indirectly embedded within the latent volatility. Given our choice of models, we test whether the structures permit the identification of volatility as a latent but influential factor, even in the absence of explicit derivative prices.

There is a challenge of understanding the ’dual role’ of matching the interest rate volatility within the ATSMs. [3] reported significant discrepancies when comparing the stochastic volatility from non-Gaussian, square-root ATSMs to the observed state variable proxy of volatility. We explore the dual role of interest rate variance by considering the following aspects. First, we establish whether the three-factor $A_1\left(3\right)$ produces a time-varying volatility and whether the variance is also a factor for the yield curve. The $A_1\left(3\right)$ model was a preferred specification choice in [8] for its ability to characterise the unconditional volatility and a flexible correlation structure. Second, there is a need to investigate the joint dynamics of level, slope, curvature and volatility.

The market price of risk specification should be chosen carefully to ensure that the tension in simultaneous fitting of certain cross-sectional and time-series properties in the ATSMs is circumvented. The three-factor ATSMs with completely affine specification by [8] was found to be poor in preserving the affine time-series and cross-sectional properties of bond prices. On the other hand, the essentially affine specification allows greater flexibility in fitting the time-varying price of interest rate risk over time, while preserving the affine time-series and cross-sectional properties of bond prices. This is despite the need for a trade-off between time-varying conditional variances and the flexibility in fitting time-varying market price of risk [1].

Estimation methods will have a problem dealing with explanation of bond yields under USV. Whereas USV can match some aspects of time series and cross-section of yields, it is unclear whether USV restrictions will affect the model’s ability to capture the cross-section and time series of yields [7]. Econometrics methods, including generalised method of moments (GMM), simulated methods of estimation (SME), among others, are unsuitable for investigating models with USV [7] and references therein. Following their approach, we consider the Bayesian Markov Chain Monte Carlo (MCMC), with a Kalman filter integrated into it. The Kalman Filter is preferred for its ability to address the non-linear processes. It is only used here as a computational devise to evaluate the likelihood given a path of the latent part $V_t$ of the state variable.

Ultimately, we test the time series and term structure dynamics using multiple volatility representations: stochastic volatility (SV) from the models, GARCH-type dynamics, and macro variables such as oil shocks and FX volatility. This comparative approach helps us understand the robustness of results across different volatility specifications and provides evidence for or against the empirical importance of USV in fixed income markets. To test the model performance and goodness-of-fit, we apply both the root-mean-square (RMSE) and biasness tests, including the pairwise comparison of [9] (DM) to the models.

The rest of this paper is structured as follows: In section 2 we review the relevant literature. Section 3 covers the model establishment. Section 4 discusses a data collection method. In Section 5 - 7 we discuss the model implementation and analysis of results. In Section 8 we conclude.

2. Literature Review

There are ongoing attempts aimed at successfully addressing the dynamic term structure models (DTSM) that exhibit USV. The presence of USV render the fixed income markets to be incomplete, therefore, present some challenges to pricing and hedging of bond prices and their derivatives. In recent studies, linear-rational square root (LRSQ) models are being evaluated for the ability to capture a cross-sectional and time series dynamics simultaneously. [10] specifies a five-factor model exhibiting USV, where two of the factors are not spanned by the yield curve. Their specification confirms the presence of USV in conditional yield variance and bond risk premia, which in turn are also linked to macro-economic uncertainty. [11] test the predictive power of the spread between short and long-term U.S government yields using both the quadratic term structure and the autoregressive gamma-zero models in comparison to LRSQ. Their results provide evidence that LRSQ allows USV, accurately fitting the cross-section of yields and capturing the time series properties of the yield variance.

Specifically, the focus on ATSMs suggest other alternatives to the trade-off between fitting the cross-section of yields and capturing their time-series behaviour. [5] considers the process where a quadratic variation of yields is tied to the cross-section of average yields. He implements arbitrary portfolios of bonds to test the presence of USV using jump-robust estimators of the diffusive variance of Nelson-Siegel factors. His approach follows [12] on quadratic diffusive and affine jump-diffusive models, although they found them to be incapable of accommodating the observed yield volatility dynamics. Their model exhibit evidence for USV across all factors, which includes oil shocks and tax policy shocks also driving part of the USV, although a greater part of variation still requires some explanation.

[7] study the "dual role" of the short rate variance as predicted by the stochastic volatility models. They provide evidence that the three-factor ATSMs explain the shape of the yield curve but fails to explain both the quadratic time-variation in spot rate as estimated from GARCH and the implied variance from options. Evidence from the four-factor model with USV suggests that the model generates both realistic short rate volatility estimates and a good cross-sectional fit. The study suggest further that short rate volatility cannot be extracted from the cross-section of bond prices. They propose an alternative representation of ATSM by [8] which is essentially characterised by — (1) physical interpretations of state variables such as level, slope and curvature, (2) affine dynamics and tractability are preserved, (3) models are econometrically identifiable, and (4) necessary and sufficient conditions lead to parameters restrictions under which the model exhibits USV [13]. In their companion paper, [2] they derive these conditions under which restriction exhibits USV. The same paper provides evidence on additional explanatory power possessed by latent variables for explaining the time-series and cross-sectional behaviour of bond prices. The presence of USV enforce the requirement that bond and derivative prices are required for estimation of parameters. Finally, USV presents some challenges on how fixed income derivatives should be hedged.

[6] explore the differences among Gaussian, stochastic volatility and USV models in the context of the conditional volatility, which is an inherent feature in ATSMs. They follow [2] approach and restrictions that ensure that models exhibit USV. They estimate the models using both Eurodollar Futures and options data. Gaussian and stochastic volatility models match the conditional mean and volatility of the term structure well. They require both the yield curve dynamics and options data for differences to be distinguished. USV models fail to resolve the tension between futures and options fits. Additional factors such as among others, jumps should be considered for the model performance to be improved. [14] maintains USV by allowing a clear separation between factors dedicated to term structure fitting and those for options fitting; see [6] and references therein. Similarly, they also reported the tension in pricing and parallel tension in matching lower and higher order moments, also suggesting the direction of adding jump components. [15] use the term structure model with inflation and economic growth factors, together with latent variables, to investigate the effect of macro variables on bond prices and the yield curve dynamics. PCA and variance decomposition show that macro factors explain up to 85% of the variation in bond yields. They explain movements at the short end and middle of the yield curve while latent factors still account for most of the movement at the long end of the yield curve. They conclude that, incorporating macro factors in a term structure model improves forecasts.

Another contributing factor to the tension between the time series dynamics and yield curve properties, hence the presence of USV, is the misspecification of the MPR. [1] proposes the essentially affine model for MPR as they are found to allow greater flexibility in fitting variations in the price of interest rate risk over time, while retaining the affine time-series and cross-sectional properties of bond prices. [16] studies a no-arbitrage specification for MPR, with their improved fit coming from the time-series rather than cross-sectional features of the yield curve. A separate set of parameters with one being risk-neutral measure and physical measure for the other set, adds flexibility to their specification. The specification relieves the tension between matching the time-series behaviour of the interest rate process and matching the cross-sectional shape of the yield curve.

[13] describes the concepts of identification, identifiability, and maximality. These are the centre for local and global identifiability of models and parameters, whose presence determine the existence of economic interpretation. Identification deals with how and whether the state vector and parameter vector can be inferred from a particular data set. Identifiability deals with the issue of whether the state vector and parameter vector can be inferred from observing all conceivable financial data, as frequently as necessary. Maximality, is about the most general model within $A_M\left(N\right)$ as discussed by [8], that is identifiable given sufficiently informative data. Affine models are written in terms of latent variables with no clear economic meaning independently from the model and parameters, leading to representations that are locally but not globally identifiable and often meaningless. [3] proposes a representation in which the state vector is written in terms of theoretically observable state variables that have unambiguous economic interpretations. To circumvent these challenges, a process is required whereby latent state variables are rotated to observable state variable so that they have economic meaning. [3] proposes a new representation of affine models in which the state vector comprises infinitesimal maturity yields and their quadratic covariations. In contrast to the invariant transformation rotations by [17] and [8], their representation accommodates the USV [18].

Finally, the issue of parameter uncertainty and model flexibility compels us to consider a probabilistic method of estimation, the MCMC. We are also concerned about the global identifiability of our models and parameters, and therefor economic meaning. [19] performs a comparative study between MCMC, Particle Filter and the Kalman Filter. Key challenges are the discretisation bias and the presence of a latent factors. While all procedures are computationally demanding, they found the Kalman Filter a computational efficient estimation method to use in goodness-of-fit testing procedures. Incorporating their process to stochastic volatility, and application to observed data suggest that the volatility depends on an additional factor that varies independently of the short rate level.

[20] proposes an efficient and flexible method to compute the maximum likelihood estimators of continuous-time models when part of the state vector is latent, considering stochastic volatility and term structure models. In contrast to MCMC, their approach relies on closed-form approximations to estimate parameters and simultaneously infer the distribution of filters—of the latent states conditioning on observations. Computational complexity of these methods may present some challenges to our intended purpose. MCMC enable us to estimate the parameters under assumptions of prior parameter distributions and conditional likelihood functions to determine posterior distributions of parameter samples and latent variables with reasonable precision.

In state-space models, the use of MCMC methods introduces sampling noise due to the discreteness of the generated samples and the inherent gaps between each pair of measurements. [21] proposes the estimation framework which relies on the introduction of latent auxiliary data to complete the missing diffusion between each pair of measurements. This data augmentation is synonymously referred to as smoothing. [22] propose a simulation-based smoothing technique and the auxiliary mixture model. The auxiliary mixture model consists of state-dependent weights and efficient block sampling algorithms to jointly update all unobserved states given latent mixture indicators [7] and references therein. The same blocking approach is followed by [7] with the parameters vector broken into three blocks— those affecting the dynamics under risk-neutral measure, risk-premium parameters and the measurement errorstandard deviations.

3. Model Establishment

ATSMs are described in terms of an $N-$dimensional Markov process of a state variable X and its dynamics under the risk-neutral measure $\mathbb{Q}$ are written as:

$$d{X}_t={\kappa }^Q\left({\theta }^Q-X_t\right)dt+\Sigma \sqrt{S_t}d{W}_t^Q$$

(1)

where ${\theta }^Q$ is the long-term mean of the process $X_t$ under the risk-neutral measure, $W_t$ is an independent Brownian motion, $\kappa$ and $\Sigma$ are $N\times N$ matrices, and $S_t$ is a diagonal matrix with the ith element given by

$$S_{ii,t}={\alpha }_i+{\beta }_i^T·X_t$$

(2)

The spot rate is an affine function of $X_t$ and written as

$$r_t={\delta }_0+{\delta }_1{\dot{X}}_t$$

(3)

where $X_t$ is a Markov N-dimensional state vector at time t, ${\delta }_0$ is a scalar and ${\delta }_i$ is an $N-$dimensional vector for $i=1,...,N$.

Provided that parameters are admissible, it is known from [17] that a zero-coupon bond has a solution

$$P_t\left(\tau \right)=e^{A\left(\tau \right)-B{\left(\tau \right)}^T·X_t}$$

(4)

where $T\equiv (T-\tau )$ and $A\left(\tau \right)$ and $B\left(\tau \right)$ satisfy the following ordinary differential equations (ODEs) also known as (Ricatti equations)

$${\displaystyle \frac{dA\left(\tau \right)}{d\tau }}=-{\theta }^{\prime Q}{\mathcal{K}}^{Q^{\prime }}B\left(\tau \right)+{\displaystyle \frac{1}{2}}\sum _{i=1}^N{\left[{\Sigma }^{\prime }B\left(\tau \right)\right]}_i^2{\alpha }_i-{\delta }_0$$

(5)

$${\displaystyle \frac{dA\left(\tau \right)}{d\tau }}=-{\theta }^{\prime }{\mathcal{K}}^{Q^{\prime }}B\left(\tau \right)-{\displaystyle \frac{1}{2}}\sum _{i=1}^N{\left[{\Sigma }^{\prime }B\left(\tau \right)\right]}_i^2{\beta }_i+{\delta }_y$$

(6)

A solution to these ODEs is found through numerical integration, starting from the initial conditions $A\left(0\right)$=$B{\left(0\right)}_{N\times 1}$.

By inverting (4) a related yield $y_t$ is computed as

$$y_t\left(\tau \right)=-{\displaystyle \frac{log{P}_t\left(\tau \right)}{\tau }}={\displaystyle \frac{A\left(\tau \right)}{\tau }}+{\displaystyle \frac{B{\left(\tau \right)}^{\prime }{X}_t}{\tau }}$$

(7)

We adopt a framework of [7] in which they use a classification scheme of [8] for $A_1\left(3\right)$ and $A_1\left(4\right)$ models, each with $M=1$ conditional volatility vector. They choose a representation where for a state vector X, each state variable has a clear economic interpretation. They derive the state variables below using the Taylor series expansions with respect to maturity $\tau$

$$Y_t\left(\tau \right)=Y_t^0\left(\tau \right)+\tau Y_t^1\left(\tau \right)+{\displaystyle \frac{1}{2}}{\tau }^2{Y}_t^2+\dots$$

(8)

They assign the state vectors $X=[r,{\mu }^Q,V]$ and $X=[r,{\mu }^Q,{\theta }^Q,V]$ for three-factor and four-factor models, respectively, where state variables are defined as follows:

$$\begin{array}{cc}\hfill r_t& =Y_t\left(0\right)\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\mu }_t^Q& =2{\displaystyle \frac{\partial Y_t\left(\tau \right)}{\partial \tau }}{|}_{\tau =0}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\theta }_t^Q& =3{\displaystyle \frac{{\partial }^2{Y}_t\left(\tau \right)}{\partial {\tau }^2}}{|}_{\tau =0}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill V_t& ={\displaystyle \frac{1}{dt}}d{r}_t^2\hfill \end{array}$$

(12)

Following [3], alternative ways of representing (1) were evaluated to identify the most general and identifiable form than the $\Sigma \sqrt{S_t}$ for the diffusion matrix. Here the stochastic components of the model are expressed in terms of a covariance matrix rather than as I$\widehat{t}o$ diffusions and is found to be more suitable for introducing the parameters that have clear economic interpretations.

3.1. The $A_1\left(3\right)$ Model

Model $A_1\left(3\right)$ is a three-factor ATSM which allows the presence of stochastic volatility of interest rates. Its canonical form is expressed by a state vector $S=[x,r,{\mu }_1]$ can be rotated into three state variables short rate, risk-neutral drift of the short rate, and the variance of the short rate, resulting into an observable state vector $X=[r,{\mu }^Q,V]$. This is regardless of whether there is or no USV; see [18]. The risk-neutral dynamics for the state vector X can be written as

Instantaneous mean:

$${\displaystyle \frac{1}{dt}}{\mathbb{E}}^{\mathbb{Q}}\left[d{X}_t\right]=\left[\begin{array}{c}{m}_0+m_r{r}_t+m_{\mu }{\mu }_t^{\mathbb{Q}}+m_V{V}_t\\ {\mu }_t^{\mathbb{Q}}\\ g_V-k_V{V}_t\end{array}\right]$$

(13)

where:

• $m_0$ = is a constant term
• $r_t={\delta }_0+{\delta }_1^{\top }{X}_t$ = feedback from the short rate
• ${\mu }_t^{\mathbb{Q}}$ = slope of the yield curve
• $m_V$ = volatility feedback

and Covariance matrix:

$${\displaystyle \frac{1}{dt}}cov\left(d{X}_t,d{X}_t^T\right)\equiv {\Omega }_t={\Omega }_0+{\Omega }_V\left(V_t-\underline{V}\right)$$

(14)

where $\underline{V}$ is set to be a lower bound for $V_t$. Admissibility requirements are also met provided that ${\Omega }_0$ and ${\Omega }_V$ are positive semidefinite and positive definite, respectively.

$${\Omega }_0=\left[\begin{array}{ccc}\underline{V}& c_{r\mu }& 0\\ c_{r\mu }& {\sigma }_{\mu }& 0\\ 0& 0& 0\end{array}\right]and{\Omega }_V=\left[\begin{array}{ccc}1& c_{r\mu }& c_{rV}\\ c_{r\mu }& {\sigma }_{\mu }& c_{\mu V}\\ c_{rV}& c_{\mu V}& {\sigma }_V\end{array}\right]$$

(15)

There are six parameters in the drift and eight in the covariance matrix, resulting into 14 risk-neutral parameters in total. This implies that the model is a Q-maximal model because of the greater number of parameters, according to [8] and [3].

3.2. Market Price of Risk

The pricing kernel which is also known as stochastic discount factor (SDF) is defined as

$${\displaystyle \frac{d{\Lambda }_t}{{\Lambda }_t}}=-r_{t}dt-{\lambda }_t^{\top }d{W}_t$$

(16)

The essentially affine form of [1] specifies the market price of risk (MPR) ${\Lambda }_t$ as a nonlinear function of volatility

$${\Lambda }_t=-\sqrt{S_t}{\lambda }_1+\sqrt{S_t^{-1}}{\lambda }_2$$

(17)

where $S_t$ is a non-negative function of the volatility state $V_t$. This structure introduces both direct and inverse sensitivity to volatility in the risk premia. It captures effects such as increased compensation for high volatility environments, and asymmetric or leverage effects through $\sqrt{S_t^{-1}}$.

3.3. Model $A_1\left(3\right)$ Under Physical Measure

No-arbitrage conditions require that both risk-neutral measure $\mathbb{Q}$ and physical or historical measure $\mathbb{P}$ be equivalent. It is implied that volatilty or diffusion part of the SDE does not change in either measure and that only drift will be influenced by the risk premia. The drift for $A_1\left(3\right)$ state vector under physical measure $\mathbb{P}$ is specified as

$${\displaystyle \frac{1}{dt}}{E}^P\left[d{X}_t\right]=\left[\begin{array}{c}{\lambda }_{r0}+{\lambda }_{rr}{r}_t+(1+{\lambda }_{r\mu }{\mu }_t^Q+{\lambda }_{rV}{V}_t)\\ (m_0+{\lambda }_{\mu 0})+(m_r+{\lambda }_{\mu r})r_t+(m_{\mu }+{\lambda }_{\mu \mu }){\mu }_t^Q+(m_V+{\lambda }_{\mu V})V_t\\ ({\gamma }_V+{\lambda }_{V0})-({\kappa }_V-{\lambda }_{VV})V_t\end{array}\right]$$

(18)

See appendix A for a detailed derivation:

3.4. The $A_1\left(4\right)$ Model

In () we defined the state variable ${\theta }_t^Q$, which represents the curvature of the yield at short maturities. Adding this to the state vector results in $X=[r,{\mu }_t^Q,{\theta }_t^Q,V_t]$. The result is a four-factor model with risk-neutral dynamics represented by the instantaneous mean and covariance matrix as

$${\displaystyle \frac{1}{dt}}{E}^Q\left[d{X}_t\right]=\left[\begin{array}{c}{\mu }_t^Q\\ {\theta }_t^Q+V_t\\ a_0+a_r{r}_t+a_{\mu }{\mu }_t^Q+a_V{V}_t\\ {\gamma }_V-{\kappa }_V{V}_t\end{array}\right]$$

(19)

and

$${\displaystyle \frac{1}{dt}}cov(d{X}_t,d{X}_t^T)\equiv {\Omega }_t={\Omega }_0+{\Omega }_V(V_t-\underline{V})$$

(20)

where

$${\Omega }_0=\left[\begin{array}{cccc}\underline{V}& c_{r\mu }& c_{r\theta }& 0\\ c_{r\mu }& {\sigma }_{\mu }& c_{\mu \theta }& 0\\ c_{r\mu }& c_{\mu \theta }& {\sigma }_{\theta }& 0\\ 0& 0& 0& 0\end{array}\right]and{\Omega }_V=\left[\begin{array}{cccc}1& c_{r\mu }& c_{r\theta }& c_{rV}\\ c_{r\mu }& {\sigma }_{\mu }& c_{\mu \theta }& c_{uV}\\ c_{r\theta }& c_{\mu \theta }& {\sigma }_{\theta }& c_{\theta V}\\ c_{rV}& c_{\mu V}& c_{\theta V}& {\sigma }_V\end{array}\right]$$

(21)

Similar to $A_1\left(3\right)$ model, admissibility requirements are met provided that ${\Omega }_0$ and ${\Omega }_V$ are positive semidefinite and positive definite, respectively. We note that the model has a total of 22 free risk-neutral parameters, qualifying the model to be a maximal $A_1\left(4\right)$.

3.5. Parameter Restrictions

[2] describe the necessary and sufficient conditions for parameter restrictions under which an $A_1\left(3\right)$ model can display USV.

3.5.1. $A_1\left(3\right)$

Implementing the following restricted parameters characterise the $A_1\left(3\right)$ as a model with USV $m_r=-2c_V^2$

$m_{\mu }=3c_V$

$m_V=1$

${\sigma }_V^2=c_V^2$

3.5.2. $A_1\left(4\right)$

[7] proposes an extension to the $A_1\left(3\right)$ restictions and derive the following for the $A_1\left(4\right)$ model to display USV [7].

$a_r=-2c_{r\mu }^2(3c_{r\mu }-a_{\theta })$

$a_{\mu }=7c_{r\mu }^2-3c_{r\mu }{a}_{\theta }$

$a_V=3c_{r\mu }$

${\sigma }_{\mu }=c_{r\mu }^2$

${\sigma }_{\theta }=c_{r\mu }^4$

${\rho }_{r\theta }=c_{r\mu }^2$

${\rho }_{\mu \theta }=c_{r\mu }^3$

3.6. Estimation STRATEGY

Bayesian estimation procedures for continuous-time finance are well documented in many sources including [7,22] and references therein, and [23]. We aim to infer the posterior distribution of unknown parameters given observed data. This process can be computationally challenging, especially when the posterior has a complex or high-dimensional structure. Two powerful tools that can simplify these challenges are data augmentation and Gibbs-like posterior sampler. Data augmentation is a technique by which a latent auxiliary variable is introduced to fill the missing diffusion between each observed pair of measurments; see [21]. The Gibbs sampler is the most simplest algorithm that is possible to use directly to sample iteratively from all complete conditionals $(\Theta ,X)$, where $\Theta$ is a parameter vector and X a latent variable [23].

[7] considers two data augmentation approaches. First, augmenting with unobservable high frequency data, which enables the use of the Euler approximation, and provide a Gaussian density that is easy to work with. This introduces discretisation bias, especially when time steps h are not sufficiently small, and the true continuous-time model deviates from the linear or Gaussian assumptions underlying Euler methods. Second, augmenenting with the observed yield data with the theoretically observable term structure factors — state variables X. Its advantage over the former is that it treats the main sources of variation in the yield curve as observable proxies for the latent state process. This approach assumes that the latent factors driving the term structure are extracted directly from cross-sectional yield data. It results in a much simpler posterior structure and avoids discretisation altogether, enabling more stable and efficient inference.

We let $\mathcal{P}={\mathcal{P}}_1,{\mathcal{P}}_2,...,{\mathcal{P}}_T$ represent a time series of principal components of the yields $Y_t$. The posterior $p(\Theta |\mathcal{P})$ is approximated using PCA-augmented data as

                                            $p(\Theta |\mathcal{P})\propto p(\mathcal{P}|\Theta )p(\Theta )$

where $\Theta$ represents a vector of parameters.

The first term on the right hand side is the likelihood function and the second is the prior distribution of parameters. The likelihood function is intractable, making it difficult to evaluate the posterior. This would require the use of several techniques from the MCMC. The observable data $\mathcal{P}$ are augmented with the term structure factor data $X=\left\{X_1,X_2,...,X_T\right\}$. These factors are theoretically observable as they can be interpreted independently from the model being considered. In practice they may not be directly observed as there are only a finite maturity yields available. Uncertainity in these state variables is therefore integrated out using a Gibb-like posterior that altrernates between drawing $p\left(\Theta |\mathcal{P},X\right)$ and $p\left(X|\mathcal{P},\Theta \right)$ [7].

True dynamics from SDE (1) are then approximated as1

$$\left(X_{t+h}\sim X_t\right)\sim \mathcal{N}\left(h(a+b{X}_t),h{\Omega }_t\right)$$

(22)

where ${\Omega }_t$ represents the instantaneous covariance matrix, and h is a discretisation time step, for which a choice of a smaller number is suitable for approximating to a normal distribution, in avoidance of a using a non-central chi-square distribution of the square root factors.

The likelihood function is then specified in terms of the relationship between data and the state vector as

                                            ${\mathcal{P}}_t=$ PC loadings $\times Y_t$

It is clear from (7) that there is a linear relation between the principal components and state variables. Adding a Gaussian error vector ${ϵ}_t\sim N(0,{\Omega }_T)$ , we obtain

$${\mathcal{P}}_t=K+L{X}_t+{ϵ}_t$$

(23)

(23) is further defined as a "measurement equation" in a state-space model, which we shall use in Kalman filtering to estimating latent factor $V_t$.

When full conditionals $(\Theta ,X)$ are not tractable for a Gibbs-sample, Metropolis-Hastings (MH) steps within the Gibbs sampler become an alternative technique. The MH algorithm generates a Markov chain with a target posterior density $\pi \left(X\right)$. At each iteration t:

• It proposes $X^{*}\sim p\left(X^{*}\right|X^{(}t))$
• Accept $X^{*}$ with the following probability:

$$\alpha =min\left(1,{\displaystyle \frac{\pi \left(X^{*}\right)p\left(X^{\left(t\right)}\right)\left|X^{*}\right)}{\pi \left(X^{\left(t\right)}\right)p\left(X^{*}\right)p\left(X^{*}\right|X^{\left(t\right)})}}\right)$$

(24)

where:

$X^{*}$ represents a proposed state or candidate in the MH algorithm

$p\left(X^{*}\right|X^{\left(t\right)})$ is the proposal distribution that defines how the candidate state $X^{*}$ is chosen given the current state $X^{\left(t\right)}$.

$\alpha$ is the acceptance probability in the MH algorithm. It determines whether the proposed state $X^{*}$ is accepted as the new state $X^{\left(t\right)}$. If $X^{*}$ is accepted, the chain moves to $X^{*}$; otherwise, it stays at $X^{\left(t\right)}$.

There are several challenges and limitations for the MCMC such as computational burden, weak parameter identification, and poor mixing. The latter is caused by the presence of high corrlation between the parameters and latent states resulting in slow convergence. [23] recommends the use of artificial data to check the efficiency and convergence of the algorithm. [22] conducted a comparative analysis on three alternative MCMC algorithms in terms of their autocorrelation functions. The option with fast decay in autocorrelation was found to have a faster convergence to the other alternatives.

Appendix C discusses briefly the estimation methodology and algorithm as adapted from [7] approach; where parameters are broken into three blocks, ${\varphi }^Q$, ${\varphi }^{\Lambda }$ and ${\varphi }^{\lambda }$ for risk-neutral measure, measurement error standard deviations and risk-premium parameters, respectively.

4. Data Collection

We use a sample of weekly SA government treasury bond prices spanning the period October 2013 to September 2024 with maturities of 3 months, 5, 10, 12, 20, 25, and 30 years. For out-sample analysis we use data for the periods October 2023 - Septermber 2024 with the same maturities. The data were retrieved from the Thomson Reuters database. We extract the approximated zero-coupon yields by inverting the bond pricing equation. Specifically, for a bond with maturity $\tau$ and price P we compute the implied yield using the simplified inversion formula (7).

Traditionally, zero-coupon yields are typically bootstrapped from swap rates, LIBOR or any of its replacements. Swap-based bootstrapping depends on credit assumptions, interpolation, and curve construction techniques that may not reflect the true market dynamics. These instruments are derived, not directly traded, and do not represent actual cash instruments. In contrast, coupon bond prices that are observed in the market, contain rich maturity and liquidity information, and are directly priced by supply-demand dynamics. We therefore use coupon bond data as the basis for our term structure modeling.

To ensure internal consistency and eliminate arbitrage opportunities, we employ a no-arbitrage ATSMs where zero-coupon bond prices satisfy

$$P(t,T)=E_t\left[exp(-{\int }_t^T{r}_{s}ds)\right]$$

(25)

Assuming an affine structure for the short rate and state variables $X_t$, the price takes the form of (4). Here, $X_t$ is a vector of latent state variables, and $A\left(\tau \right)$, and $B\left(\tau \right)$ are functions satisfying Riccati differential equations, which ensure closed-form and arbitrage-free bond pricing.

We simulate the evolution of interest rates using the SDE (1). At time t, given each simulation paths, bond prices are computed, followed by zero yields. These zero yields produce a time series of term structures $y(t,\tau )$, each curve reflecting the simulated latent state at time t for a maturity $\tau$.

To extract the structure from the evolving yield curves, we apply the PCA. The zero yield data are then decomposed into orthogonal components to capture the majority of variance with a few statistical factors level, slope and curvature. Table 1 presents the PCA for our zero yields data. It is noted that the first three components have a cumulative variance of 99.97% which agrees with the empirical evidence of [24]; see Figure 1 below.

We also compute the FX volatility and oil shock using brent crude and USDZAR prices, respectively, for the period October 2013 to September 2024.

5. Scenario Determination

We apply the ATSM to extract the interest rate volatility from a cross-section of bond prices. ATSM represent both the yield curve term structure and the time series of bond prices over several trading periods. This multi-period and cross-sectional nature of the interest rates present several challenges including the volatility and trading risk, risk premia and market completeness. At the heart of term structure modelling, is the need to study both the term structure and time series simultaneously. ATSMs have been successfully used in this area as they are characterised by tractability, closed-form solutions, efficient approximation and closed-form moment conditions for empirical analysis [13] and references therein. Our study is based on the following choice of scenarios, some of which are already supported by the empirical evidence:

• The test whether state vectors and parameters vectors have a unique economic interpretation. This requires a suitable representation among ATSMs where latent state vectors are translated into observable factors. Both the state vector and model parameter vector should be globally identifiable so that their values can be compared directly across different countries, periods and even models.
• Among the three and four factor stochastic volatility models, evaluate their capability to break the dual role of predicting the variance of the short rate and simultaneously a linear combination of yields and the quadratic variation of the spot rate. There is empirical evidence that $A_1\left(3\right)$ model is unable to play the dual role [13]. We compare the $A_1\left(4\right)$ USV with $A_1\left(3\right)$USV.
• To determine whether estimation based on bond price only or a combination of both bond price and options data produce best results. In the absence of option price data, we test as an alternative the simulated at-the-money bond futures implied volatility, and the macro variables as sources of variation and the substitute for options data when estimating USV models.

Among several estimation approaches, we proceed with the Bayesian approach which combines the MCMC and Kalman filter. This probabilistic approach is more suitable where there are parameter uncertainty and model flexibility issues.

6. Model Implementation

We use (4) to invert zero yields from the SA government treasury bond. The zero yields are further transformed into nonlinear state space form. The purpose is to implement the three-factor and four-factor affine models to extract the short rate variance . The final outcome is to determine whether the affine models can explain simultaneously the time series and cross-sectional properties of bond prices. It should be noted that the ATSMs that we selected are based on the maximal forms of both $A_{(}3)$ and $A_1\left(4\right)$.

[3] propose an alternative approach to those of [17](DK) and [8] to rotate from latent to observable state vectors. Whereas DK’s technique involves inverting the term structure with respect to latent factors and cannot be implemented, [3] write these term structure dynamics in a nonlinear state space. The model parameters are then estimated by Bayesian and MCMC.

It is noted that the volatility state variable $V_t$ does not enter the bond pricing equation for those models exhibiting USV. A four-factor model with a state vector $X=[r,{\mu }_t,{\theta }_t,V_t]$ effectively becomes a three-factor model by excluding the state variable $V_t$. $V_t$ is therefore an additional volatility factor that is free to explain the time series patterns [13]. USV can also be regarded as a latent component $V_t$ that influences the risk premia and yield variance but is not directly priced by the cross-section of bond yields.

Despite the exclusion of $V_t$ from the bond price equation, it continues to influence the following, thus fulfilling the conditions for USV:

• The conditional variance of state transitions $\Sigma \left(X_{t+1}\right|X_t)$
• The market price of risk via ${\Lambda }_t=-\sqrt{S_t}{\lambda }_1+\sqrt{S_t^{-1}}{\lambda }_2$

Since $V_t$ is not spanned by the cross-section of bond yields, it must be identified through time-series variation and external instruments.

In our empirical implementation, we filter $V_t$ using the observed term structure which, by way of robustness check we compare against the variation of macro signals such as oil shocks, and GARCH-based FX volatility. This we believe, also assist to ascertain the identifiability and economic relevance of the USV factor in markets with limited market options data due to sparse trades.

We consider a state-space model:

$$\begin{array}{cc}\hfill y_t& =K+L{X}_t+{ϵ}_t,{ϵ}_t\sim \mathcal{N}(0,{\Sigma }_{ϵ})\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill X_{t+1}& =\mu +\Phi X_t+{\eta }_t,{\eta }_t\sim \mathcal{N}(0,Q\left(V_t\right))\hfill \end{array}$$

(27)

$X_t$ is an unknown state vector,and $y_t$ is the observed data vector. Both equations represent the Kalman filter. We apply the Kalman filter within the MCMC procedure to extract the posterior distribution of $V_t$.

7. Analysis of Results

This section presents the results from the two competing models $A_1\left(3\right)$ and $A_1\left(4\right)$ both with USV, focusing on both in-sample and out-of-sample performance. We begin with the analysis of posterior distribution for parameters to assess the point estimates and associated intervals for estimation consistency. This is followed by an assessment of the models’ ability to fit the yield curve in-sample, thereafter, an evaluation of their forecasting performance out-of-sample. Subsequently, we provide a closer examination of the time series behavior of key the latent factor, volatility, term structure implications, various sensitivity and robustness diagnostics.

7.1. Posterior Distributions of Key Parameters

Table 2 and Table 3 report the parameter point estimates for risk-neutral and risk premia, respectively. For each point estimate, based on mean values there are a corresponding credible interval bounds, which appear consistent across the two models. Point estimates lie well within the interval bounds. The variation between models is expected due to differences in model structure, but the overall consistency suggests stability in estimation. We applied the discretisation value $h={\displaystyle \frac{1}{52}}$, which produced reasonable estimates.

In addition, for a sample of risk-neutral drift parameters, which includes the physical ${\gamma }_V^P$ and $K_V^P$ we analysed the posterior plots for consistency.2 Figure 2 and Figure 3 plot the histograms and kernel densities for drift parameters for models $A_1\left(3\right)$ USV and $A_1\left(4\right)$ USV, respectively. As expected, the shapes of the histograms do not exhibit a multi-model but instead, are all unimodal. This is a sign of convergence and good mixing of the MCMC algorithm.

To gain more insight on the effectiveness of the Bayesian inference obtained through MCMC sampling, we used the python software function arviz.summary() of [25] to generate the Table 4 and Table 5. These tables present posterior summaries for key model parameters obtained via MCMC sampling for models $A_1\left(3\right)$USV and $A_1\left(4\right)$USV, respectively. For each parameter, we report the posterior mean, standard deviation, and 95% Highest Density Interval (HDI), along with diagnostic statistics including Monte Carlo standard errors (MCSE), effective sample sizes (ESS), and the Gelman–Rubin convergence statistic $(\widehat{R}$). All $\widehat{R}$ values are equal to 1, indicating satisfactory convergence. Effective sample sizes (ESS) were used to assess sampling efficiency and convergence. In Table 4, bulk ESS values ranged from 975 to 2550 across parameters, with all but one parameter exceeding 2000. The lowest bulk ESS (975) is still well above common diagnostic thresholds $\ge 400$, which is indicative of adequate sampling.

Tail ESS values ranged from 674 to 1500, demonstrating reliable estimation of posterior distribution tails and credible intervals. Similar trend for bulk ESS and tail ESS is also observed in Table 5. These values suggest that the Markov chains mixed well and that the posterior estimates are statistically stable. The effective sample sizes are sufficiently large, suggesting reliable estimation. Notably, all point estimates lie well within their respective HDI bounds, consistent with stable posterior inference.

7.2. Yield Curve Fit

In this section we evaluate the in-sample and out-sample model performance by way of the root mean square error (RMSE) and bias for each maturity. Model-fitted yields are derived from the point estimates in Table 2 and Table 3. Similar to [7] we avoid the process of integrating over the posterior distribution for computational burden. Posterior samples are therefore rerun only for sampled state variables and parameter values held fixed.

Using the short rate observed during short-term while assuming the constant yield to maturity, we compute the discount rate for each maturity which we then apply to convert the observed yields into par bonds. From the par bonds we compute the zero yields using (7). This exercise is necessary so that we may compare both the model-fitted yields and the observed ones as zero yields.

Table 6 presents the yield curve fit analysis results in terms of RMSE and bias, both in and out-sample3. The yield curve fitting performance of two models is assessed using both RMSE and forecast bias, evaluated in-sample and out-of-sample across maturities 0.25, 5, 10, 12, 20, 25, 30 years. We applied the DM test to compare forecast errors ${ϵ}_t=Y_t-\widehat{Y_t}$ and detect statistically significant differences; where $Y_t$ and $\widehat{Y_t}$ represent the observed and forecast yields, respectively.

In-sample, model $A_1\left(4\right)$ USV provides a significantly better fit to the yield curve in terms of RMSE, with a DM test statistic of 6.1233 and p-value of 0.0009. This suggest that $A_1\left(4\right)$ USV captures the cross-sectional variation in yields more accurately. Model $A_1\left(3\right)$ USV exhibits significantly a lower forecast bias, particularly at longer maturities, with a DM statistic of -8.00 and p-value of 0.0002, suggesting that the model aligns more closely with observed yields on average, despite being less precise.

Out-sample, the results are more nuanced. Model $A_1\left(4\right)$ USV exhibit lower RMSE at shorter maturities, with a statistically significant difference observed at the first maturity. Across the full maturity spectrum, the difference in RMSE between the models is not statistically significant with a DM statistic of 0.9999 and p-value of 0.3559. A similar pattern is seen in forecast bias with model $A_1\left(3\right)$ USV consistently showing lower bias, but the overall difference is again not significant with a DM statistic of 0.9571 and p-value of 0.3755.

Finally, the results suggest that $A_1\left(4\right)$ USV fits the in-sample yield curve more accurately, whereas $A_1\left(3\right)$ USV exhibits less bias and potentially more stable performance out-of-sample. There is a trade-off between precision and biaseness whose decision depend on the application, as we shall also discuss under the volatility regression below. Although previous studies have associated model misspecification with high autocorrelation in residuals [7], we did not perform an autocorrelation analysis in this study. Future research could incorporate this test to further assess model adequacy.

7.3. Time Series Dynamics

In this section we evaluate the properties of model-implied state variables. To estimate them, we run a posterior sampler by holding the parameters in Table 2 and Table 3 fixed. Smoothed estimates are obtained by averaging the results of the draws of state variables [7]. Both models produce a time series of state variables $X_t$, including a filtered volatility $V_t$. From these model-implied state variables, we compute various correlations against the actual yield time series variables. A 26-week rolling standard deviation of the log returns for the 3-month actual yields is fitted. As a comparative means, we also repeated the computation of the 26-week rolling standard deviation for the 12-year and 30-year maturities. Yield curve variables, slope was computed as a function of the 30-year and 3-month yield as $\left(Y_{30y}-Y_{0.25y}\right)$, and the curvature as $\left(Y_{30y}-2Y_{12y}+Y_{0.25y}\right)$. Other variables include the GARCH(1,1) volatility of the same 3-month log returns data used for computing the 26-week rolling standard deviation, FX volatility for SA Rand Dollar currency pair, and Brent crude volatility.

Table 7 compares models $A_1\left(3\right)$USV and $A_1\left(4\right)$USV in terms of the various correlations between model-implied variables and the actual time series variables for 3-month, 12-year and 30-year maturities. The $A_1\left(3\right)$ model shows strong alignment with yield curve decomposition literature of [24], by capturing over 99% of the variance in the yield data. All three components level, slope, and curvature show very high correlations with their empirical counterparts with over 96%, —especially curvature with 98.9%, consistent with theoretical expectations. Figure 1 illustrates how the first three components contribute 99.97% to total variation. As result, 0.03% variance explained is well below any reasonable threshold and negligible.

The $A_1\left(4\right)$USV model introduces a marginal component explaining only 0.03% of variance. Including this factor raises the average yield fit to $r_t=1.0$ due to overfitting but significantly worsens the curvature correlation —down to 0.733, likely due to numerical instability introduced by a non-meaningful fourth factor. Whereas, the $A_1\left(4\right)$USV model slightly improves average yield, it does so at the cost of degrading the structural interpretation of key term structure components.

The $A_1\left(4\right)$USV model exhibits a larger correlation of 0.550 between the model-implied and 26-week rolling for the short end, even though it decreases for both mid and long end. The $A_1\left(3\right)$USV model on the other hand exhibits very low or even negative correlations of 0.054, 0.002 and -0.037 for the short, mid and long end, respectively. The GARCH(1,1) displays the results that are almost similar except that suprisingly some slight increase in correlations is observed with the $A_1\left(4\right)$USV model.

The left panel of Figure 4 plots the 26-rolling volatility and the model-implied volatility $V_t$ against the period 2016 - 2024. Figure (a) displays the inability of $A_1\left(3\right)$USV model volatility to track the variation of short rate volatility with 26-week rolling, instead the model exhibits a flat pattern. The $A_1\left(4\right)$USV volatility on the other hand, tracks the variation of the short rate volatility closely with the 26-week rolling volatility.

As a robustness check, we compare the correlations between the model-implied volatility $V_t$ and FX volatility, and Brent crude. Correlations for the 3-month displays a positive relationship with the model-implied volatility $V_t$ of 0.408 for SA Rand Dollar, and 0.216 for the Brent crude volatilities. A poor picture is exhibited for the mid and short end of the yield curve regarding the correlations. For the yield curve factors level, slope and curvature there is generally a positive relationship with the model-implied volatility $V_t$.

The curvature factor exhbits negative correlation with model-implied volatility for both models, with $A_1\left(3\right)$USV showing slightly higher correlation of -0.119 when compared to more negative -0.421 for the $A_1\left(4\right)$USV models. The negative correlations between the model-implied volatility and the curvature factors— a yield curve component— confirm the existense of the tension between time series dynamics and yield curve dynamics.

7.4. Volatility Forecasting and Regression

In this section we analyse the in-sample and out-sample performance using two volatility proxies, the absolute one-week of changes in yield returns $E\left[\parallel \Delta Y\parallel \right]$ and realised volatility $\widehat{\sigma }$. We use the same data for parameter point estimates from Table 2 and Table 3 to extract our weekly forecasts.

Following the methodology adopted in [7], we implement a rolling-window approach using a two-year or approximately 104-week estimation window to construct one-week-ahead forecasts of the state variable. At each step, model parameters obtained from a posterior sampler are held fixed. The rolling process is used to re-estimate the latent state variables based on incoming data. This recursive structure allows for dynamic updating of state variables while maintaining stable parameter estimates. Forecasts are generated for both a 104-week in-sample period and a 45-week out-of-sample period. This design enables robust comparison of forecast accuracy across different subsamples and simulates a realistic forecasting scenario.

To generate each forecast, we first estimate the current values of the latent state variables using only information available up to time t This is done in a manner consistent with the previous section, except that in the forecasting context, no data beyond time t are used in the state filtering step. Importantly, the parameter estimates used in this process are fixed at their posterior modes or means, obtained from the full sample estimation. Given these current-state estimates, we simulate 10000 paths of the model forward one week to generate the predictive distribution of the state variables at time $t+1$. From this simulated distribution, we then construct a forecast distribution for each yield and volatility proxy of interest per maturity.

The second proxy, realised volatility is given as

$$\widehat{{\sigma }_{t,\tau }}=\sqrt{{\sum }_{i=1}^N\Delta Y_{t,i,\tau }^2}$$

(28)

where $Y_{t,i,\tau }$ is the $\tau$-maturity yield on the ith week following the observation at t.

Traditionally, realised volatility is better estimated using high-frequency return data such as daily to increase the certainity over reliability and accuracy of the forecast. In this study we use weekly yield data to construct both proxies for volatility. Our choice is motivated both by data availability and by the structure of the domestic bond market, where SA government treasury bonds are traded during weekly auctions. It is our view that these auctions anchor the pricing of yields, meaning that weekly data capture economically meaningful shifts in investor expectations, interest rate risk, and macroeconomic signals. Intermediate level frequencies such as weekly or monthly are better than both high-frequency and low-frequency levels in producing sign dependence in volatility, which is the source of forecastability in asset returns [26].

7.4.1. Forecasting and Model Performance

We report both in-sample and out-of-sample RMSEs in Table 8. In each case, the DM test is employed to assess the statistical significance of differences in forecast accuracy between models. The in-sample RMSE of $\parallel \Delta Y\parallel$ exhibits a mixed picture. For the first three maturities, model $A_1\left(4\right)$USV show a more reduction in RMSEs than $A_1\left(3\right)$USV, indicating an improved fit. For the remaining four maturities, $A_1\left(3\right)$USV performs better. Despite this maturity-dependent performance, the overall DM test yields a test-statistic of 0.053 with a p-value of 0.959, indicating that the observed differences are not statistically significant at conventional levels.

Turning to the in-sample RMSE of the estimated volatility series $\widehat{\sigma }$, $A_1\left(4\right)$USV consistently outperforms $A_1\left(3\right)$USV across all maturities. The DM statistic of 3.976 and a p-value of 0.007 support this result with significance at the 1% level. This suggests that $A_1\left(4\right)$USV more accurately captures the in-sample volatility dynamics of yields.

Out-sample, $A_1\left(4\right)$USV again shows improved performance in forecasting $\parallel \Delta Y\parallel$, with uniformly lower RMSEs. The DM statistic of 2.404 and a p-value of 0.053 suggest a statistical significance at the 5% level. This result supports $A_1\left(4\right)$USV’s superior ability to forecast near-term volatility in yield changes. Out-sample performance of the volatility forecasts is even more decisive. $A_1\left(4\right)$USV demonstrates markedly lower RMSEs in forecasting $\widehat{\sigma }$ with a DM statistic of 6.791 and a p-value of 0.001. This result is highly significant and highlights the robustness of$A_1\left(4\right)$ USV in capturing the underlying volatility dynamics beyond the estimation window.

It is worth noting, a tension between time series volatility performance and cross-sectional yield curve dynamics. While $A_1\left(4\right)$USV consistently outperforms $A_1\left(3\right)$ USV in forecasting volatility, both in-sample and out-sample, the initial inconsistency across maturities in the RMSE of $\parallel \Delta Y\parallel$ reflects a possible trade-off. Specifically, models that are optimised for volatility dynamics may not always align perfectly with those tailored for fitting the cross-section of the yield curve.

This tension is not surprising. Yield curve models often prioritise cross-sectional fit at a point in time, while volatility models emphasise dynamic consistency and predictive power over time. The results here suggest that $A_1\left(4\right)$USV, potentially incorporating richer dynamics or additional latent volatility structure, sacrifices some cross-sectional fit in the near term, as seen in maturities 12, 20,25 and 30-year but gains substantially in time series forecasting accuracy.

Ultimately, the decision between models depends on the application such as pricing and hedging for which cross-sectional accuracy may dominate. For risk management and forecasting, superior volatility dynamics as achieved by $A_1\left(4\right)$USV are likely to be more valuable.

We also notice that both models match the unconditional volatility of yield changes, with $A_2\left(4\right)$USV exhibiting a more refined pattern. The right panel of Figure 4 displays the unconditional volatility plotted against maturities, standard deviation of model-fitted yield changes $\Delta Y_t$ in blue and lying within the light-grey bound distribution of model-implied yield volatilities as determined by the point estimates. The model-fitted unconditional volatility is clearly fitting along the mean and median of the distribution. In the bottom plot (d) we note that $A_1\left(4\right)$ USV appears to be more of an improvement from $A_1\left(3\right)$ USV in the top plot (b). In $A_1\left(4\right)$ USV the model-fitted unconditional volatility ( blue line) matches the pattern of the distribution average (grey line) better when compared to $A_1\left(3\right)$ USV. Both plots exhibit a snake-shaped curve as discussed by [4] for US Treasury yields and swaps. From an improved plot (d), $A_1\left(4\right)$ USV the pattern commencing with the back of the snake, followed by a hump towards 5-year maturities, a hump and drops towards the 10-year maturity, thereafter decreasing in a stable manner towards longer maturities. High volatility over the period between 3-month to 5-year maturities, also termed by [4] the back-of -the-snake, is due to reactions from monetary policy events, changes in liquidity premia and macro economic dynamics, among others. It is also a key to the factor correlations.

7.4.2. Regression

To evaluate the extent to which model-derived volatility $V_t$ can be explained by both term structure dynamics and exogenous volatility sources, we estimate the following linear regression model:

$$V_t=\alpha +{\beta }^T{X}_t+{\gamma }_i{\Sigma }_{m_{i,t}}+{ϵ}_t$$

(29)

where $X_t$ denotes the yield curve factors derived from either the three-factor with level, slope, and curvature or a four-factor including latent volatility $V_t$ in term structure model. ${\Sigma }_{m_{i,t}}$ represent the M observable macro variable proxies at time t for $i=1,...,M$, and ${\gamma }_i$ is the coefficient for the i-th macro variable.

Table 9 presents the results of volatility regressions for the 3-months and 30-year volatilities. The GARCH model exhibits moderate explanatory power, with $R^2$ values between 0.221 and 0.477, suggesting economically meaningful relationships in both short- and long-end volatility regressions, regardless of the inclusion of PCA factors. Applying PCA on both models $A_1\left(3\right)$USV and $A_1\left(4\right)$USV significantly improves the model performance, increasing $R^2$ values from approximately 0.022 to 0.428 for both models, with minimal differences observed between them. Since both models achieve similar explanatory power after PCA with $R^{2}0.43$, $R^2$ alone cannot be used to select the better model.

Stochastic volatility models, $A_1\left(3\right)$USV and $A_1\left(4\right)$USV exhibit weak explanatory power in their original specification with $R^2$ values of 0.022 and 0.027, but performance improves markedly after applying PCA to $R^2$ values of 0.428 and 0.427, indicating that latent yield curve factors significantly enhance the models’ ability to capture volatility dynamics. Despite similar fit post-PCA, $A_1\left(3\right)$USV demonstrates a much stronger negative relationship between volatility and the yield curve factors with $\beta$ coefficient of –6.876, while $A_1\left(4\right)$USV exhibits a modest positive relationship with $\beta$ of 0.810, suggesting fundamentally different volatility structures. The estimated loadings on the level, slope, and curvature components further reveal that $A_1\left(4\right)$USV is more responsive to changes in yield curve shape, especially in terms of slope –0.173 and curvature –0.364. These findings imply that while both models satisfy the USV condition by showing significant latent factor influence, $A_1\left(4\right)$USV may better align with USV theory by isolating volatility innovations not captured by the term structure, though $A_1\left(3\right)$USV captures a stronger overall volatility response. In the second regression, $A_1\left(3\right)$USV and $A_1\left(4\right)$USV exhibits a beta of -2.144 and 0.302, respectively, suggesting that $A_1\left(4\right)$USV is still positively related to realised volatility and aligns better. A drop in $\beta$ coefficient from 0.810 to 0.302 only suggests that the model $A_1\left(4\right)$USV performs better in the short- than longer end.

The purpose of this regressions is to evaluate the extent to which the stochastic volatility models capture information embedded in the yield curve, and whether they span the term structure or USV properties. Two key insights emerge from both the 3-month and 30-year maturity regressions. First, time-series information, extracted through PCA, proves substantially more informative than raw cross-sectional data. In their unaugmented forms, both models exhibit weak explanatory power with $R^2$ approximately 0.022–0.027, indicating that volatility is poorly explained by contemporaneous yield levels alone. However, once latent factors level, slope, and curvature are introduced, explanatory power improves substaintially with $R^2$ increasing to 0.428 for 3-month volatility forecast and up to 0.221 for the 30-year forecasts. This highlights the importance of dynamic yield curve behavior in driving interest rate volatility.

Second, although both models are influenced by latent yield curve factors, they diverge in structure and responsiveness. $A_1\left(3\right)$ USV displays a strong contemporaneous volatility response, with large negative coefficients on latent factors. Yet, it shows limited sensitivity to yield curve shape in the long end regression. This implies that $A_1\left(3\right)$ USV primarily reacts to broad level shifts, lacking a nuanced view of the term structure, particularly over longer horizons.

7.4.3. Market Price of Risk

Figure 5 plots a comparison of MPR components from the $A_1\left(3\right)$ and $A_1\left(4\right)$ USV models. The $A_1\left(3\right)$ model yields relatively stable and low-magnitude MPRs, while the $A_1\left(4\right)$ USV model exhibits consistently higher values and greater variability, suggesting that the additional factor introduces a more pronounced risk compensation mechanism. This difference, largely a result of model specification and scale, highlights how additional factors can amplify perceived risk premia. Given that both sets of MPRs are derived from Bayesian posterior point estimates under the risk-neutral measure, the results reflect not just data fit but also the influence of prior beliefs. Overall, the 4-factor model implies a richer, possibly more realistic, structure for pricing risk.

8. Conclusion

We assume that restricted ATSMs with USV can be effectively estimated through Bayesian methods in emerging markets using only time-series data. As a result, we do not follow a traditional joint approach of estimation using both bond price and options data. We evaluate the sensitivity of oil shock and exchange rate volatility to the interest rate volatility. We test how the restricted models $A_{(}3)$USV and $A_1\left(4\right)$USV respond to the tension between cross-sectional and time series dynamics of bond prices.

The tension between time series and yield curve dynamics does not go away completely. Evidence of this is highlighted by a negative correlation between the curvature versus both volatility and variance. Model $A_1\left(4\right)$USV performs better than $A_1\left(3\right)$USV in capturing the time time series dynamics. There is a neccesary trade-off between the time series or yield curve fitting, which depends on the application — option pricing, portfolio management, risk management and hedging and policy formulation.

Model flexibility and parameters uncertainity issues are reduced by the Bayesian MCMC estimation strategy. The result is that parameters become locally and globally identifiable, and economically meaningfull such that they are easily compared to parameters used in data from other countries.

Our study did not assess the auto-correlation of residuals in bias and RMSE due to computational burden. We did not factor the macro variables into a joint modelling but only used them as a robustness test mechanism. Future research should consider a joint modelling of macro variables as a means to compensate the role of options data in a joint estimation and evaluate their impact on parameter estimation in the context of Bayesian modelling.