South African Government Bond Yields and the Specifications of Affine Term Structure Models

1. Introduction

The yield curve continues to be a useful tool for the financial market participants as a guide to levels of interest rates and their economic expectations. Synonymously, it is referred to as the term structure of interest rates, a representation of the yields plotted against a cross-section of unexpired maturities. Bonds and yields are used as a relative pricing benchmark for other financial instruments. Traders of various fixed income securities such as credit default swaps, loans, and related derivatives among others, rely on the yield curve to price their trades. The yield is a discount rate at which cash flow of fixed income securities is equated to the present value of an instrument a point in time. The yield is dependent on short-term interest rates, which in our SA context is a short-rate such as a repo set by the Monetary Policy Committee of the SA Reserve Bank is a proxy.

The yield is linked to a combination of either observable and unobservable (latent) factors, which may be either macroeconomic such as inflation, gross domestic product, yield curve factors, or the three statistical factors level, slope, and curvature as per the empirical evidence by [2]. They evaluate methods of constructing yield curve factors and their loadings such as principal component analysis, the dynamic Nelson-Siegel and affine term structure models (ATSMs). The yield curve factors in particular, refer to the instantaneous or short rate of interest changes and its volatility.

Our study, which motivates the choice for ATSMs, seeks to address the cross-sectional yield dynamics with special attention to conditional volatility and correlation in a multi-factor setting. This would enable us to study several active SA Treasury bonds of different maturities simultaneously. Moreover, these dynamics have a critical impact on the pricing of bond options, risk management (hedging), and measurement such as calculation of value-at-risk. Stochastic volatility models of [1,3], and [4] among others, consider the interest rate volatility to be a latent factor from which many properties of the yield may be observed.

Recent study by [5] models the SA government bond yield using the interpolation method of Nelson-Siegel. Their report cites a remarkable statement that SA bond market is considered to be a leader within the emerging markets, and also ranked the highest in terms of liquidity. In addition, our review of the BIS quarter 4 of 2023 statistics highlights an amount of $\$238.6$ billions of debt securities issued by the SA central government, in all markets at all original maturities denominated in domestic currency at nominal value; see [6]. Whereas the Nelson-Siegel method produced good-fitting results in their study, our aim is to exploit the econometric modelling capabilities within ATSMs to model the SA bond yields. ATSMs provide a tractable and flexible framework to describe this dynamic behaviour of the entire term structure of interest rates by modelling the evolution of short-term and long-term rates simultaneously. They are known to have certain mathematical properties that make them analytically tractable. In particular, they can accommodate stochastic volatility, jumps and correlation among risk factors that drive asset returns. Affine representations of state variables in the DAPMs are popular because they lead to computationally tractable pricing relations and moment equations that can be used in estimation; see [7].

The traditional practice of capturing the bond yield movement over time together with some macroeconomic variables has been conducted by using methods such as the vector autoregression. However, there are limitations when it comes to treatment of certain aspects of the bond and yield. First, a bond is an asset, and that the same bond with several different maturities may be traded at the same time. Second, there is a risk of holding long-dated bonds for a short-term as the investors may require compensation for such risk. There is an expectation that long-term yields represents a risk-adjusted average of a cross-section of short-term yields plus a risk premium in the absence of arbitrage. Third, the yields are not normally distributed, thus making it difficult to compute the expected values of future short rates. ATSMs provide a solution for this problem; see [8].

ATSMs commence with the presumption that the instantaneous short rate $r\left(t\right)$ is an affine function of an $N-$dimensional state vector Y, $r\left(t\right)={\varphi }_0+{\varphi }_Y^{\prime }Y\left(t\right)$, and that $Y\left(t\right)$ followed Gaussian and square-root diffusions, respectively. Affine models were extended by [9] into consistent and arbitrage-free multifactor model of the term structure of interest rates in which yields at selected fixed maturities follow a parametric multivariate Markov diffusion process with "stochastic volatility." [1] explore the structural differences and relative goodness-of-fit of the ATSMs. In particular, they consider a trade-off between the following modelling issues. First, the economic representation of both short and long-term dynamics of the state variable when studied simultaneously. Second, the computational burden of estimation and curve fitting. This trade-off is formalized by their classification of $N-$factor affine family into $N+1$ non-nested subfamilies of models. Their special attention to three-factor ATSMs suggests, based on theoretical considerations and empirical evidence, that some subfamilies of ATSMs are better suited than others to explaining historical interest rate behaviour. One example is the case where the conditional probability density of yields in closed-form is non-existent, rendering the maximum likelihood of no use. To circumvent this, a feasible choice should come from estimation methods such as Fourier-based, generalised methods of moments (GMM), Markov chain Monte Carlo (MCMC) and simulated methods of moments (SMM) among others. It is known that the conditional likelihood function of the latent state vector $Y\left(t\right)$ may not be known, as a result [1] follow the SMM of [10].

The remainder of the paper is structured as follows: Section 2 capturers the review of relevant literature. Section 3 describes the ATSMs. Section 4 defines a canonical representation of the ATSM and restrictions that are imposed on the parameters. 5 An overview and descibes the three-factor ATSMs. 6 Reviews briefly the characteristic function-based estimation methods for the ATSMs. Section 7 - 10 Discusses data collection, scenario determination, model implementation and analysis of the results. In section 11 we conclude.

2. Literature Review

This paper explores the behaviour of the SA interest rates in terms of historical time series and a cross-section of yields across a maturity spectrum. Inspired by the seminal work of [1], we proceed by implementing their model and its maximal counterpart. Among three options of models tested on US Treasury swaps, their $A_1\left(3\right)$ was found to perform better than other models, followed by their $A_2\left(3\right)$. Initially, they consider a comprehensive framework for the specification, analysis, and classification of ATSMs. They provide a complete characterisation of admissible and identified ATSMs from which it is required that sufficient general conditions exist; see [11], who describe the regular affine process. They also characterise the sufficient general conditions that must be met for a process to be affine; see [7,8] among others.

ATSMs are among popular models in the vast literature on interest rates term structure and bond pricing. Few examples are the early generation consisting of a single-factor Gaussian of [12], and a square-root process by [13], extended by [14] into a multi-factor. The next generation are the correlated mixture affine models of [1,9], among others. The reason for their popularity is the ability to accommodate stochastic volatility, jumps and correlations among factors driving the asset returns; and lead to computationally tractable closed-form prices, and estimation through moment equations; see [7]. Among research problems addressed using ATSMs is the description and treatment of the co-movement of short and long-term bond yields. An affine process Y is defined as one in which a conditional mean $\mu$ and variance $\sigma {\sigma }^{\prime }$ are affine functions of Y. The process is further defined and characterised by [11] as regular affine process, a class of time-homogeneous Markov processes. They consider a state space $D={\mathbb{R}}_{+}^M\times {\mathbb{R}}^N$, for integers $M\ge 0$ and $N\ge 0$, from which the logarithm of a characteristic function of a transitional probability $p\left(x_t\right)$ of such process is affine with respect to the initial state $x\in D$. [7] conveniently formalises it in terms of their exponential-affine Fourier (for continuous-time) and Laplace (for discrete-time) transforms. The affine relationship is defined by coefficients which are solved by a family of ordinary differential equations (ODEs). These ODEs are the essence of tractability of regular affine processes. [9] apply the ODEs as time-dependent drivers of the solution to a zero-coupon bond, provided the parameters are admissible. An inverted form of these zero-coupon bond gives rise to a yield as a state variable. They also exploit the idea of a yield-only analysis without including additional economic variables as latent factors.

[15,16] are among several authors who have approached the application of ATSMs in discrete-time although they are known to have less popularity compared to their continuous-time counterparts. Earlier models exhibited a tendency of perfectly correlated returns of bonds of all maturities, which is an unrealistic behaviour and unsuitable for hedging; see [17]. Several authors extended these one-factor Markov representation of a short-rate by introducing a range of multi-factor models with the long-run mean $\theta \left(t\right)$, and the stochastic volatility $\upsilon \left(t\right)$ of $r\left(t\right)$ that are affine functions of $\left(r\right(t),\theta (t),\upsilon (t\left)\right)$ for which [1] explores several specifications. [18] endorse a parsimonious representation of the yield curve matching the time series and cross-sectional variation of bond yields through three-factor models. They develop a simple estimation approach by exploiting the exponential-affine structure of these models; see also [19] on the stochastic mean and stochastic volatility and three-factor model of the term structure of interest rates and its applications in derivatives pricing and risk management

A specification of an ATSM should be "admissible" and therefore lead to well-defined bond prices. The admissibility property is completely characterised by [11] in the "canonical" state space $D={\mathbb{R}}_{+}^M\times {\mathbb{R}}^N$, with a non-negative diagonal matrix. However, this property has a problem of imposing parameter restrictions on the affine process to ensure that it is well defined. One typical scenario is the restriction of parameters to ensure that the conditional variance of a state variable remains non-negative. The requirements for admissibility become more complex as the number of state variables determining conditional variances increase; see [7]. The admissibility condition ensures that the process does not exit the domain $D\in {\mathbb{R}}^N$. A family of $A_M\left(N\right)$ models with a domain $D={\mathbb{R}}_{+}^M\times {\mathbb{R}}^{N-M}$ are a common admissible family of models; where M factors evolve in a positive state space while $N-M$ evolve in an unrestricted space; see [20]. [1] verifies this easily through admissible $N-$factor ATSMs that are uniquely classified into $N+1$ non-nested subfamilies.

Admissible models should also be canonical, meaning that they are economically identified, and maximally flexible; see [7]. As a result, the $A_M\left(N\right)$ benchmark ATSM models should have a canonical representation and also satisfy the non-negative and non-explosive solution of ([21]). Their drift should satisfy a Lipschitz condition, and the diffusion should satisfy the uniqueness condition of [22]; see [8]. These conditions have an effect of restricting the correlation structure of the affine diffusions. Exploiting the Gaussian and square-root form of diffusions, there still appear to be non-satisfaction of the regularity conditions of non-explosive growth and uniqueness, giving rise to need for a Feller condition;1 see [8]. A multi-dimensional extension of a Feller condition was implemented by [9], which was found to handle the general correlated affine diffusions. The condition ensures that only positive factors enter the volatility $\sigma \left(y\right)$. This involves restrictions on the state variables that prevent the instantaneous conditional variances $S_{ii}\left(t\right)$ from becoming negative. This condition is sufficient for the existence of a unique solution to the affine SDE according to [11].

For each of the $N+1$ subfamilies, there exists a maximal model that is econometrically plausible for all other models within this subfamily. They describe further the maximal models in relation to the $N+1=4$ classification; and highlight an interaction within the family of ATSMs between the dependence of the conditional variance of each $Y_i\left(t\right)$ on $Y\left(t\right)$ and the admissible structure of the correlation matrix for Y. A key advantage of maximal models is that of overcoming the overidentifying restrictions that are imposed on yield curve dynamics; see [1]. The admissibility property is also confirmed by the no-arbitrage solution for a zero-coupon bond following [9].

[1] specification applied the continuous-time approach to the ATSMs which is popular to a majority of empirical literature. They explore the structural differences and relative goodness-of-fits of ATSMs. They refer to a trade-off between flexibility in modeling the conditional correlations and volatilities of the risk factors. They classify a family of $N-$factor affine into $N+1$ non-nested subfamilies of models. From their three-factor ATSMs, empirical analysis suggests that some subfamilies of ATSMs are better suited than others to explaining historical interest rate behaviour.

The focus of the research is to implement the specifications of [1] to test the pricing of zero-coupon bonds and forecasting the yield curve dynamics when using the SA bond yield. It also attempts to extract the latent factors from the yield itself, without any consideration for other economic factors; see [9]. ATSMs are proven to dominate both theoretical and empirical frameworks in term structure modelling; see [8]. A link between the cross-sectional and time series properties is made consistent by the ATSMs. Evolution of unobserved factors from the risk-neutral dynamics of the yield are proved to have both the drift and the diffusion coefficients as affine functions of such factors by the ATSMs; see [8]. Several methods of estimation are available and require mostly the knowledge of the joint conditional density of yields. In this study, we follow the estimation method of Fourier inversion for the characteristic function of a state variable, which is assumed to lead to a conditional density. This method leads to a closed-form solution where the maximum likelihood is an efficient estimator.

3. Model Establishment

We discuss the model in the context of admissibility of ATSMs. In the absence of arbitrage opportunities, a zero-coupon bond that matures at time T is priced as

$$\begin{array}{c}\hfill P(t,\tau )=E_t^Q\left[e^{-{\int }_t^{\tau }r\left(s\right)ds}\right]\end{array}$$

(1)

where:

• $P(t,\tau )$ is the price of a bond at time t maturing at time $\tau$.
• t is the current or initial time at which the bond is evaluated.
• $\tau$ is the maturity date, at which the bond pays its face value.
• s is a continuous time variable at which the interest rate process $r\left(s\right)$ evolves.2
• $E_t^Q[·]$ denotes the conditional expectation under the risk-neutral measure Q given the information available at time t.

To obtain an $N-$factor ATSMs it is assumed that an instantaneous short-rate $r\left(t\right)$ is an affine function of a vector of N unobservable state variables $Y\left(t\right)=Y_1\left(t\right),Y_2\left(t\right),...,Y_N\left(t\right)$, written as

$$\begin{array}{c}\hfill r\left(t\right)={\delta }_0+\sum _{i=1}^N{\delta }_i{Y}_i\left(t\right)\\ \hfill ={\delta }_0+{\delta }_y^{\prime }Y\left(t\right)\end{array}$$

(2)

where ${\delta }_0\in \mathbb{R}$ and ${\delta }_y\in {\mathbb{R}}^N$.

Another assumption is that $Y\left(t\right)$ follows an affine diffusion

$$dY\left(t\right)={\mathcal{K}}^Q({\theta }^Q-Y\left(t\right))dt+\mathrm{\Sigma }\sqrt{S\left(t\right)}d{W}^Q\left(t\right)$$

(3)

${\mathcal{K}}^Q$ and ${\theta }^Q$ represent the reversion rate and central tendency (long-term mean) parameters under a risk-neutral measure, respectively. $W^Q\left(t\right)$ is an $N-$ dimensional independent Brownian motion under the risk-neutral measure Q; $\mathcal{K}$, and $\mathrm{\Sigma }$ are $N\times N$ matrices, which may be asymmetric or non-diagonal. $S\left(t\right)$ is a diagonal matrix with the $ith$ diagonal elements written as

$${\left[S\left(t\right)\right]}_{ii}={\alpha }_i+{\beta }_i^{\prime }Y\left(t\right)$$

(4)

where ${\alpha }_i\in \mathbb{R}$ and ${\beta }_i\in {\mathbb{R}}^N$. The parameter $\alpha$ can be interpreted as an intercept, which represents the base or long-run level of the variance for the $i^{th}$ component. ${\beta }_i$ represents the sensitivity of the variance for the $i^{th}$ component to the state vector $Y_t$. Together, ${\alpha }_i$ and ${\beta }_i$ ensure that the conditional variance is always positive. The non-negativity in (4) is the core requirement for admissibility in this framework.

The drifts in (3) and conditional variances in (4) are both affine in $Y\left(t\right)$. [9] has the following time dependent solution to the price of a zero-coupon bond, provided that parameters are admissible.

$$P(t,\tau )=e^{A\left(\tau \right)-B{\left(\tau \right)}^{\prime }Y\left(t\right)}$$

(5)

and the related yield is computed as

$$y(t,\tau )=-{\displaystyle \frac{logP(t,\tau )}{\tau }}={\displaystyle \frac{A\left(\tau \right)}{\tau }}+{\displaystyle \frac{B{\left(\tau \right)}^{\prime }Y\left(t\right)}{\tau }}$$

(6)

where $A\left(\tau \right)$ and $B\left(\tau \right)$ are coefficients whose solution satisfies the following ODEs (Ricatti equations)

$${\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[{\mathrm{\Sigma }}^{\prime }B\left(\tau \right)\right]}_i^2{\alpha }_i-{\delta }_0$$

(7)

$${\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[{\mathrm{\Sigma }}^{\prime }B\left(\tau \right)\right]}_i^2{\beta }_i+{\delta }_y$$

(8)

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}$. Risk-neutral dynamics of the short rate $r\left(t\right)$ in (2) through to (4) determine this specification of the ODEs.

To use the closed-form representation of (1) in the empirical study of ATSMs, it is required that the distributions of $P(t,\tau )$ and $Y\left(t\right)$ under actual physical measure P be known. To this end, a market price of risk $\mathrm{\Lambda }\left(t\right)$ is introduced as

$$\mathrm{\Lambda }\left(t\right)=\sqrt{S\left(t\right)}\lambda$$

(9)

where $\lambda$ is an $N\times 1$ vector of constants. The process $Y\left(t\right)$ under physical measure P, therefore also has an affine form 3

$$dY\left(t\right)=\mathcal{K}(\mathrm{\Theta }-Y\left(t\right))dt+\mathrm{\Sigma }\sqrt{S\left(t\right)}dW\left(t\right)$$

(10)

Note that a superscript Q has been removed. $W\left(t\right)$ is an $N-$dimensional vector of independent Brownian motion under P, $\mathcal{K}={\mathcal{K}}^Q-\mathrm{\Sigma }\Phi$, $\mathrm{\Theta }={\mathcal{K}}^{-1}({\mathcal{K}}^Q{\theta }^Q+\mathrm{\Sigma }\psi )$. $\Phi$ comprises of ${\lambda }_i{\beta }_i^{\prime }$ in its $ith$ row, and $\psi$ is an $N-$vector with ${\lambda }_i{\alpha }_i$ as its $ith$ element.

[1] acknowledge that their main purpose of (9) is to preserve the affine structure of $Y_t$ under $\mathcal{Q}$. They do not pursue the impact of the market price of risk on the forecast but only focus on the correlation and volatility dynamics of the state variables. The square-root process is followed here to ensure non-negativity to the price of risk. Their form is referred to as "completely affine", which is found to have limitations as far as pricing of risk is concerned. As a results, some authors extended the completely affine form of market price of risk to address various factors that impact on the price of risk itself; see [25,26] and [27] among others. Our approach focusses on the econometric representation of state variables, which we assume to incorporate a detailed market price of risk. In spite of this, we follow in the footsteps of [1] where the main focus was to address the specification problems considering both correlation and volatility on the cross-section of yield data. Our workings would not focus on analysing the market price of risk but rather leave it to the future research which might incorporate other forms suitable for SA yield curve. We also support the idea of extending a detailed market price of risk to the emerging markets to respond to factors such as volatility, liquidity and bond credit risk; see the recent report by [28]. Their report addresses a joint modelling of liquidity and credit risk for SA bond market, where they recommend their model or similar approach for emerging markets and corporate bonds.

4. A Canonical Representation of ATSMs

According to [1], a general specification for (10) may not always lead to a positive conditional variances over a range of Y, given an arbitrary set of parameters $\psi =\left(\mathcal{K},\mathrm{\Theta },\mathrm{\Sigma },\mathcal{B},\alpha \right)$. However, admissibility requires that parameters restrict $S_{ii}\left(t\right)$ in (10) to be strictly positive for all i; where $\mathcal{B}$ denotes the matrix of coefficients on Y in $S_{ii}\left(t\right)$.

From (4), there is a special case where there is no admissibility problem when ${\beta }_i=0$, for all i, since the instantaneous conditional variances are all constant. Outside the special case, it is necessary to impose constraints on the drift parameters $\mathcal{K}$ and $\mathrm{\Theta }$, and diffusion coefficients $\mathrm{\Sigma }$ and $\mathcal{B}=\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_N\right)$. The requirements for admissibility become more restrictive as a number of state variables determining $S_{ii}\left(t\right)$ increases.

They consider a case where there are M state variables driving the instantaneous conditional variance of the $N-$vector Y, such that $M=rank\left(\mathcal{B}\right)$. They further propose a set of $N+1$ benchmark models $A_M\left(N\right)$ as the most flexible econometrically identified affine DTSM on the state space $R_{+}^M\times R^{N-M}$; see also([11]). It is only when the admissibility conditions are met that a canonical representation may be defined.

Definition 3.1: 

For each M, $Y_t$ is partitioned as $Y^{\prime }=\left(Y^{\prime V},Y^{\prime D}\right)$ where $Y^V$ is $M\times 1$, and $Y^D$ is $(N-M)\times 1$; where V and D represent the volatility sources and the dependent factor, respectively. The canonical representation of the benchmark model $A_M\left(N\right)$ is defined as a special case of (3) with

$$\mathcal{K}=\left(\begin{array}{cc}{K}_{M\times M}^{VV}& 0_{M\times (N-M)}\\ K_{(N-M)\times m}^{DV}& K_{(N-M)\times (N-M)}^{DD}\end{array}\right)$$

(11)

for $M>0$, and $\mathcal{K}$ is either lower or upper triangular for $M=0$.

The canonical representation of $\mathcal{K}$ is the mean-reversion matrix, with diagonal terms expected to pull the mean level to non-negativity, thus influencing positive variances. Its off-diagonal terms on the other hand reflect how different state variables influence each other, indicating potential dependencies or interactions that could affect the overall system behaviour. The matrix $\mathcal{K}$ therefore captures both the stabilising effects of the mean-reversion rates and the dynamic interplay between different state variables. In the three-factor analysis, this trade-off between non-negative variance and correlations requires a special attention. It also has an impact on the choice of M, number of state variables entering volatility and the interactions among $N=3$ factors

$$\mathrm{\Theta }=\left(\begin{array}{c}{\mathrm{\Theta }}_{M\times 1}^V\\ 0_{(N-M)\times 1}\end{array}\right)$$

(12)

$$\mathrm{\Sigma }=I$$

(13)

$$\alpha =\left(\begin{array}{c}{0}_{M\times 1}\\ 1_{(N-M)\times 1}\end{array}\right)$$

(14)

$$\mathcal{B}=\left(\begin{array}{cc}{I}_{M\times M}& V_{M\times (N-M)}^{VD}\\ 0_{(N-M)\times M}& 0_{(N-M)\times (N-M)}\end{array}\right)$$

(15)

The following parameter restrictions are imposed:

$${\delta }_i\ge 0,M+1\le i\le N$$

(16)

$${\mathcal{K}}_i\mathrm{\Theta }\equiv \sum _{j=1}^m{\mathcal{K}}_{ij}{\mathrm{\Theta }}_j>0,1\le i\le M$$

(17)

$${\mathcal{K}}_{ij}\le 0,1\le j\le M,j\ne i$$

(18)

$${\mathrm{\Theta }}_i\ge 0,1\le i\le M$$

(19)

$${\mathcal{B}}_{ij}\ge 0,1\le i\le MM+1\le j\le N$$

(20)

[1] define a subfamily $A_M\left(N\right)$; of affine DTSM as nested special cases of the $Mth$ canonical model or its invariant transformation; where $M=(0,\cdots ,N)$ . Equivalent affine models are obtained under invariant transformations that preserve admissibility and identification and leave the observable quantities like short rate unchanged. Details on invariant transformation are discussed in Appendix A of [1].

• The following issues are further noted from [1]:

The assumed structure of $\mathcal{B}$ ensures that $rank\left(\mathcal{B}\right)=M$ for the $Mth$ canonical representation. To verify that M resides in $A_M\left(N\right)$, instantaneous conditional correlations among $Y^V\left(t\right)$ are zero, whereas the instantaneous correlations among $Y^D\left(t\right)$ are determined by parameters ${\mathcal{B}}_{ij}$ because $\mathrm{\Sigma }=I$. Admissibility is established provided (20) holds, and that the conditional covariance matrix of Y depends only on $Y^V$. Zero restrictions in the upper right $M\times (N-M)$ block of $\mathcal{K}$ and the constraints in (18) and (19) ensure that $Y^V$ is positive. Stationarity is also assured by ensuring that all the eigenvalues of $\mathcal{K}$ are strictly positive; see also Appendix B in [1].

In addition to an admissible canonical representation, in which the minimal known sufficient condition for admissibility were imposed, minimal normalisations for econometric identification are imposed to derive a "maximal" model in $A_M\left(N\right)$. A more unique class of maximal $A_M\left(N\right)$, referred to as the equivalence class of $A{M}_M\left(N\right)$ model is obtained by invariant transformation of the canonical representation; see Appendix A in [1].

[1] further points that the canonical representation of $A_M\left(N\right)$ models may not always be a practical way for analysing state variables in ATSMs. Often, existing literature opted for parameterising ATSMs with the riskless rate r as a state variable, resulting in "affine in r" (Ar) representation. This can be rewritten as an "affine in Y" $\left(AY\right)$, where $r\left(t\right)$ can be expressed as an unobserved state vector $Y\left(t\right)$. As a result, a thorough specification analysis for $N-$factor ATSMs necessitates evaluating $N+1$ non-nested, maximal models, and ensuring that a thorough understanding of the model’s structure and implications is obtained.

5. The Three-Factor ATSMs

Three-factor models are used to describe the historical behaviour of the term structure of interest rates. Traditionally, these factors are unobserved (latent) and can only be defined statistically using techniques such as principal component analysis to convey economic meaning. Popular yield-curve fitting approaches such as the dynamic Nelson-Siegel model apply the principal component analysis (PCA) loadings to fit a yield curve; see [2]. These approaches appear to fit and forecast well but lack the theoretical rigour to enforce some no-arbitrage restrictions. Contrary to the yield-curve fitting approaches, the empirical approaches to the factor models such as the ATSMs are worth pursuing. They consider the maximal parameterisation through which in general the economic identification of factors can be revealed. [18] is among early works that are based on enforcing the no-arbitrage restrictions by implementing the three-factors models. They constructed a simple affine model with short term interest rate, mean rate and volatility as three factors, which are easy to estimate. They further conclude that the short rate plays and important role in yield curve modeling, following their observation that it could not be dominated by any other factor across all maturities.

[1] explores various forms of the canonical ATSMs and their maximal counterparts, as influenced by the number conditional volatility and correlation of factors. Fixing these factors into $N=3$ gives rise to their three-factor models which posit mainly the representation of the short-rate itself, its mean rate and volatility as the three-factors. Analysis and comparisons are made of the Gaussian versus the square root diffusion forms of the models, even though the latter appear to be preferred as it imposes the non-negative variance restrictions.

Three-factor models were derived from the notation $A_M\left(N\right)$ where M is the number of state variables that enter volatility $S\left(t\right)$ according to [1]. Emphasis has been put on the trade-off between conditional volatility and correlation as a focus for the analysis of the term structure of interest rates. As previously discussed, [9] introduced a multi-dimensional Feller condition. It ensures that negative state factors do not enter the volatility $S_{ii}\left(t\right)$ by restricting correlations; see also [8]. We have previously also discussed in a similar context, the role of the mean-reversion rate matrix $\mathcal{K}$, its non-negative diagonal terms restrictions and interactions among state variables through its off-diagonal terms.

A number M of the factors that drive the process which enters a volatility $S\left(t\right)$ become the main argument on the choice of an $A_M\left(N\right)$ model, depending on the purpose of the study. [29] point out that more volatility factors result in less flexibility in allowing risk premium and correlation structure. As a result, they are in favour of a one conditional volatility factor models $A_1\left(N\right)$ such as the $A_1\left(3\right)$ by [1]. [30] also favour the $A_1\left(N\right)$ with $N=3$ and $N=4$ for the same purpose of allowing flexibility for risk premium and correlation. Their focus is to impose restrictions on the parameters of $A_1\left(3\right)$ such that the volatility factor $\nu \left(t\right)$ disappears from the bond pricing equation. In our approach, we analyse the admissibility of parameters and cross-equation restrictions that result from interactions among the factors $Y\left(t\right)=\left(\nu \right(t),\theta (t),r(t\left)\right))$. As previously discussed, the mean-reversion rate matrix $\mathcal{K}$ has elements with either negative of positive magnitudes playing a role of ensuring that factors are pulled from entering the variance only for non-negative values, otherwise non-negative correlations are the result. This is also applicable in the case of a three-factor model.

In this study we focus on the $A_1\left(3\right)$, $A_2\left(3\right)$ models and their maximal counterparts $A{M}_1\left(3\right)$ and $A{M}_2\left(3\right)$ to determine both the fit and estimation when applied to the SA bond yield curve.

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

These models are characterised by one factor Y as a source of conditional volatility. As a result, $M=1$ gives rise to the model form of $A_1\left(3\right)$. From the original [18] BDFS model, the $A_1\left(3\right)$ according to the notation of [1] is specified as

$$\begin{array}{r}\hfill d\nu \left(t\right)=\mu (\overline{\nu }-\nu \left(t\right))dt+\eta \sqrt{\nu \left(t\right)}d{B}_v\left(t\right)\\ d\theta \left(t\right)=\nu (\overline{\theta }-\theta \left(t\right))dt+\zeta d{B}_{\theta }\left(t\right)\\ dr\left(t\right)=\kappa (\theta \left(t\right)-r\left(t\right))dt+\sqrt{\nu \left(t\right)}d{B}_r\left(t\right)+{\sigma }_{r\nu }\eta \sqrt{\nu \left(t\right)}d{B}_{\nu }\left(t\right)\end{array}$$

(21)

where:

• ${\sigma }_{r\nu }={\displaystyle \frac{{\rho }_{r\nu }}{\eta \sqrt{1-{\rho }_{r\nu }^2}}}$

The state variables $\nu \left(t\right)$, $\theta \left(t\right)$ and $r\left(t\right)$ are the stochastic volatility for $r\left(t\right)$, central tendency or long-run mean of $r\left(t\right)$ and short rate processes, respectively. The volatility affects the short rate through its volatility factor $\eta$. The coefficient $\kappa$ represents the rate at which the short rate reverts to the central tendency. The stochastic volatility $\nu \left(t\right)$ also enters $r\left(t\right)$ and it is also instantaneously correlated with $r\left(t\right)$ as noted in the last term ${\sigma }_{r\nu }\eta \sqrt{\nu \left(t\right)}d{B}_{\nu }\left(t\right)$.

The maximal model is best suited for interpreting the parameter restrictions. As a result, [1] prefer the following model in (23) as a maximal $A{M}_1\left(3\right)$ which is affine in r. They determine their $A_1\left(3\right)$ by relaxing the parameters ${\sigma }_{\theta r}$ and ${\sigma }_{r\theta }$ in order to accommodate a non-zero correlation between the short rate and central tendency. All the other parameters inside the square boxes are set to zero to impose significant restrictions on the dynamics of interest rates and their volatility.

(22)

(23)

where:

• $\nu \left(t\right)$ serves as stochastic volatility for $r\left(t\right)$, but also enters the drift of $r\left(t\right)$, and correlated to $r\left(t\right)$ as noted in the term ${\sigma }_{r\nu }$; $\theta \left(t\right)$ is the central tendency of r and $\kappa$ is the rate at which the short rate reverts to its central tendency. Appendix E in [1] describes the transformation framework from which a test for admissibility and canonical representations in $A{M}_1\left(3\right)$ can be achieved.

5.2. $A_2\left(3\right)$

These models are characterised by two factors of Y as a source of conditional volatility. As a result, $M=2$ gives rise to the model form of $A_2\left(3\right)$. The [19] model is the member of this sub-class of models, and it is represented as

$$\begin{array}{c}\hfill d\nu \left(t\right)=\mu (\overline{\nu }-\nu \left(t\right))dt+\eta \sqrt{\nu \left(t\right)}d{W}_1\left(t\right)\\ \hfill d\theta \left(t\right)=\nu (\overline{\theta }=\theta \left(t\right))dt+\zeta \sqrt{\theta \left(t\right)}d{W}_2\left(t\right)\\ \hfill dr\left(t\right)=\kappa (\theta \left(t\right)-r\left(t\right))dt+\sqrt{\nu \left(t\right)}d{W}_3\left(t\right)\end{array}$$

(24)

$W_1$,$W_2$ and $W_3$ are independent Brownian motions. The $\theta$ follows a square-root diffusion unlike in the case of the BDFS model. Other parameters $\nu$, $\eta$ and $\kappa$ are the same as in the above models. These leads us to the convenient maximal model for $A_2\left(3\right)$ which is represented as

(25)

(26)

[1] relaxes the restrictions on ${\kappa }_{\theta \nu }$, ${\kappa }_{r\nu }$ and ${\sigma }_{r\nu }$, while other parameters within a square box are restricted to zero.

6. Estimation for Affine Models

Several estimation strategies such as the maximum likelihood, generalised method of moments, simulated method of moments, Markov chain Monte Carlo, and the characteristic function-based method are discussed by many authors; see [7,8,31]. We mention three among a possible list of issues to consider when selecting an estimation strategy for the affine models. First, an infinite set of moment conditions can cause a stochastic singularity problem, which leads to constraints for the GMM. This result from a cross-sectional yield data with many maturities. Second, a choice between inclusion or non-inclusion of a measurement error in a representation that links the observed yield with a state variable. Third, maximum likelihood efficiency is dependent on the conditional density of the state variable which is not always known; see [32].

In contrast to the maximum likelihood estimation which requires that the density functions must be computed, CCF-based methods are straight-forward. They depend on the knowledge of functional form of the CCF for variables that are observed from affine diffusions. CCFs are the foundation for computationally tractable and asymptotically efficient estimators of the parameters of affine diffusions and asset pricing models representing the affine state variables; see [31].

It is generally known that the conditional density function f of $Y_{t+1}$ has a solution up to an inverse Fourier transform of ${\psi }_t(\varphi ;\gamma )$

$$f\left(Y_{t+1}\right|Y_t)={\displaystyle \frac{1}{{\left(2\pi \right)}^N}}{\int }_{R^N}{e}^{-i{\varphi }^{\prime }{Y}_{t+1}}{\psi }_t(\varphi ;\gamma )d\varphi$$

(27)

The charaterisic function ${\psi }_t(\varphi ;\gamma )$ for $Y_{t+1}$ given $Y_t$

$$\psi \left(Y_{t+1}\right)=E\left[e^{i{\varphi }^{\prime }{Y}_{t+1}}\right|Y_t]$$

(28)

From the Proposition 1 of [33], it can be shown that under suitable regulations (5) is the conditional characteristic function of $Y_{t+1}$, with $A\left(\tau \right)$ and $B\left(\tau \right)$ derived from the solution of (7) and (8) for $\tau =1$. Therefore, the conditional characteristic function becomes

$$\psi \left(\tau \right)=e^{A\left(1\right)+B{\left(1\right)}^{\prime }{Y}_{t+1}}$$

(29)

The log-likelihood form for (27) becomes

$$l_T\left(\gamma \right)={\displaystyle \frac{1}{T}}\sum _{t=1}^{T}log\left\{{\displaystyle \frac{1}{{\left(2\pi \right)}^N}}{\int }_{R^N}{e}^{-i{\varphi }^{\prime }{Y}_{t+1}}{\psi }_t(\varphi ;\gamma )d\varphi \right\}$$

(30)

By conjecturing the parameters $\gamma$ and computing the Fourier inversion, maximum likelihood can be obtained by maximising (30), to obtain a maximum likelihood estimator by characteristic function (ML-CCF); see [7].

[7] considers densities of individual columns $Y_i$, for $i=1,...,N$. A selector vector $l_j$ is assigned an entry 1 and zero elsewhere. The density f of $y_{j,t+1}=l_j·y_{t+1}$ given the entire $y_t$ is the inverse Fourier transform of ${\psi }_t(\varphi l_j;\gamma )$

$$f\left(Y_{t+1}\right|Y_t)={\displaystyle \frac{1}{{\left(2\pi \right)}^N}}{\int }_{R^N}{e}^{-i\varphi l_j^{\prime }{Y}_{t+1}}{\psi }_t(\varphi l_j;\gamma )d\varphi$$

(31)

Estimation of (31) is based on one-dimensional N integrations instead of $N-$dimensional integrations.

Alternatively, the general method of moments (GMM) using a characteristic function is achieved by the residual

$$ϵ(t+1,\varphi ;\gamma )=e^{-i{\varphi }^{\prime }Y(t+1)}-\psi (t,\varphi ;\gamma )$$

(32)

For an arbitary instrument $z(t,\varphi )$; the estimator becomes

$$l_T={\displaystyle \frac{1}{T}}\sum _t{\int }_{R^N}z(t,\varphi )e^{-i{\varphi }^{\prime }{Y}_{t+1}}-\psi (t,\varphi ;\gamma )$$

(33)

The GMM approach is a beter alternative to a multi-dimensional Fourier inversion. However, as a grid of $\varphi$ becomes finer, correlations among moments become increasingly large, leading to a singular distance matrix; see [7].

For an affine DTSM, a link between a set of N-dimenstional yields $y_t$ of several maturities $\tau$, it follows from 6 that $y_t=A\left({\gamma }_0\right)+B\left({\gamma }_0\right)Y_t$; where $Y_t$ follows an affine diffusion, A is an $N\times 1$ vector, B is an $N\times N$ matrix. Vector ${\gamma }_0$ is a set of parameters linking $Y_t$ to the affine representation under the risk-neutral measure $\mathcal{Q}$.

A solution to the latent variable $Y_t$ can be solved, provided $B\left({\gamma }_0\right)$ is invertible, as

$$Y_t=B{\left({\gamma }_0\right)}^{-1}(y_t-A\left({\gamma }_0\right))$$

(34)

By the standard change-of-variable analysis, the conditional density function of $y_t$ under P becomes

$$f_y\left(y_{t+1}\right|y_t:\gamma )=f_Y\left(B{\left(\gamma \right)}^{-1}[y_{t+1}-A\left(\gamma \right)]\right|y_t^{\prime }\gamma )abs\left|B{\left(\gamma \right)}^{-1}\right|$$

(35)

If the conditional density function $f_Y$ of a state vector is known, it is easy to continue with the estimation of parameters $\gamma$. For special cases of continuous-time Gaussian and independent square-root processes, $f_Y$ is known as long as the choice of the market price of risk is chosen to suite an affine process under P; see [7] and references therein.

For the continuous-time affine models in the family $A_M\left(N\right)$$(M\ne 0,N)$, the unknown $f_y$ can be computed from the CCF. It is easy to express the CCF of $y_t$ in terms of the CCF of $Y_t$, since $B\left(\gamma \right)$ is nonsingular and both $Y_t$ and $y_t$ can generate the same information set; see [7].

For an affine diffusion it can be shown that CCF of $y_t$ is

$${\psi }_{y_t}(u,\gamma )=e^{i{u}^{\prime }A\left(\gamma \right)}{\psi }_{Y_t}\left(B{\left(\gamma \right)}^{\prime }u\right)$$

(36)

The Fourier inversion becomes

$$f_y\left(y_{t+1}\right|y_t;\gamma )={\displaystyle \frac{1}{{\left(2\pi \right)}^N}}{\int }_{R^N}{e}^{i{u}^{\prime }B\left(\gamma \right)Y_{t+1}}{\psi }_{Y_t}\left(B{\left(\gamma \right)}^{\prime }u\right)$$

(37)

Fourier transforms are more suitable for low dimension problems, as they become numerically burdensome as N increases. The burden could be avoided by selecting a method of moments, though at the expense of econometric efficiency; see [7].

From (34) it is clear that the measurement errors in the yields were excluded. As a result, $B\left(\gamma \right)$ can be inverted to compute the state vector $Y_t$, with Kalman filtering becoming very usefull; see [8] and references therein. Kalman filtering provide the best solution for extracting nonlinear state vectors from the affine diffusions. We do not discuss the Kalman filter and its different forms in this paper even though we apply it to filter out the state vector $Y_t$. There are many sources for a detailed discussion on the Kalman filter; see [23].

7. Data Collection

A sample of weekly yields for active SA government treasury bonds over the periods October 2013 to September 2024 with maturities 3 months, 5, 10, 12, 20, 25, and 30 years, were retrieved from the Thomson Reuters database. Our in-sample and out-sample data were based on the periods October 2013 - Sep 2023 and October 2023 - September 2023, respectively. The out-sample will be best-suited for forecasting and validation. A summary of descriptive statistics for yields across maturities is presented in Table 1. Mean values range between 6.6% to 10.4% exhibiting an upward slope which is also convex in shape. Recent study by [5] reported a similar behaviour for the average yields. The highest weekly standard deviation of 1.3% is observed for the 3-months maturity which is typical in the short end of the yield curve, suggesting that yields may react quickly to changes in monetary policy or market sentiment. It is followed by a drop to 0.7% and 0.6% for the 5-year and 10-year maturities, respectively. The 20-year maturity exhibits a rise in weekly standard deviation to 0.9%, remaining constant towards the long-end. This suggests variation and higher volatility in the short end, followed by relative stability for the long end of the yield curve.

Table 2 presents the correlations across maturities of yields. The short end of the yield curve is characterised by weak correlations. The 3-month and 5-year terms exhibit negative correlations with their long end counterparts. This may suggest a difference and diversity in dynamics between the short and long end of the yield curve. It is also possible among reasons for the negative correlation that high liquidity in SA bond market has a portion of foreign investors who are favoured by falling currency exchange rates. SA bond is also found to have long maturities when compared to its emerging market counterparts; see [28]. Other possible pertaining to positive correlations towards the long end should be associated with lower volatility and large portions of pension fund portfolios investing the same bonds with similar maturities but less trading activity. From the 10-year maturity, higher correlations ranging from 0.612 to 0.999 are observed, suggesting stability as the yield curve approaches its long end.

Figure 1 presents the first three principal components (PC) of the yields over maturities. They are explained by the variance of about 99.70%, which is close but slightly abovet the 98% empirical finding according to [34]; see also [2,8]. The first PC represents a key rate shift or level change in rates. It is the result of volatility causing rates of all maturities to fluctuate by almost the same amount. The observation is that short end is associated with high volatility and increasing rates. At mid-term around the 10-year maturity they reach a peak, followed by stability as they approach the long end of the yield curve. The second PC represents a slope which exhibits a downward slope but with its highest level in the short end which might be associated with rate increases and volatility, followed by a drop in rates in mid-term region and stable but falling rates towards the long end. Volatility forces a fall in rates in the short end followed by a drop towards mid-term, 10-year maturity, thereafter, stabilises towards the long end. Both the slope and curvature exhibit the downward and convex slope in the same direction, suggesting a stylised fact of high volatility in the short end and low towards the long end.

The behaviour of PCs is also associated with the correlations as discussed earlier, where short end is associated with weak correlations whereas in the long end of the yield curve, strong correlations were observed. These patterns are suitable for trading in swaps and correlation-based hedging strategies. Our focus being the ATSMs, we believe that these PCs are somehow closely related to the latent factors derived by the solution of coefficients $A\left(\tau \right)$ and $B\left(\tau \right)$ in (6). The three-factor models with three labels short-rate, volatility and central tendency should exhibit nearly a similar pattern to PCs; see [8].

Figure 2 presents a time series of three principal components level, slope and curvature as derived from the observed yields. They are compared to the Kalman filtered time series of state vectors, namely factor 1, factor 2 and factor 3. Both time series are plotted against the observations period of October 2013 - Sep 2023. Both plots have similar patterns with factors on the right hand exhibiting smoother shapes than the principal components on the left. The analysis is extended to conducting a regression between the Kalman-filtered factors (dependent variable) and principal components (PC), with results presented in Table 3. We note a high cumulative variance (94.3%) , suggesting that the first few PCs explain the yields variance. Intercept coefficients are close to zero for all the PCs, suggesting less impact in explaining behaviour of movements relative to the PCs. A sharp jump from 84.1% to 95.8% is a confirmation that the first few PCs captured most of the information in the yields data, with the later components contributing progressively less. The standard errors also vary across rows, with some being quite small, which indicates that the estimates of the corresponding coefficients are relatively precise. It can be concluded that the level, slope and curvature display similar characteristics to the factors from an affine model, therefore confirming the empirical findings by [34].

Figure 3 presents a selection of observed yields from the SA Treasury bonds plotted against maturities of up to 30 years. A spread between 5 and 20 years is also plotted and expected to represent a slope. Crossovers are observed among individual plots from time to time, indicative of either positive or negative (inverted) yield curves. Initial unobserved inputs to the three-factor simulation of $\nu \left(t\right),\theta \left(t\right),r\left(t\right)$ are set to the initial values of the first three principal components. The length of a full matrix of yields comprising of 418 weekly observations and seven maturities is based on the maturity column with a minimum length. No gaps in data were discovered, otherwise the ommisions would be filled by applying an average of any preceeding two values. We apply the first vector of the yield matrix as an initial input values together with the initially guessed parameters to simulate the state variables from (10). These are further used as inputs to the solution of ODE (7) and (8) from which coefficients $A\left(\tau \right)$ and $B\left(\tau \right)$ are obtained. Thereafter, a model-based set of zero-coupon bonds and zero-yields are obtained from (6). The selection is also guided by the observations from the PCA, suggesting that our selection is a proxy for level, slope and curvature, taking into consideration the PCA features for short end, mid-term and long end.

8. Scenario Determination

Our objective is to fit the ATSMs using the SA bond yield curve to extract the latent factors by way of zero-yields. From these zero-yields, we pursue a three-factor approach to model specifications in the following scenarios.

• Evaluate the performance of model $A_1\left(3\right)$ and its maximal counterpart from in-sample. Ensure that parameters are admissible and that models meet the canonical form in-sample
• Evaluate the performance of model $A_2\left(3\right)$ and its maximal counterpart from in-sample. Ensure that parameters are admissible and that models meet the canonical form in-sample
• Estimate both models and their counterparts out-sample and evaluate their performance.

Further, we draw more insight from optimisation results and statistical analysis. In all the scenarios, we assume the market price of risk is completely affine and of the form $\mathrm{\Lambda }\left(t\right)=\sqrt{S\left(t\right)}\lambda \left(t\right)$.

9. Model Implementation

9.1. Three-Factor Models

The first three principal components are selected as proxies to the required three factors. Their initial values are inputs to both the $A_1\left(3\right)$ and $A_2\left(3\right)$ and their maximal counterparts. To calibrate the three-factor parameters we select the 3-month, 10-year and 20-year observed yields. Solution to SDE is initialised by the initial state vector and parameters to obtain $\nu \left(t\right)$, $\theta$ and $r\left(t\right)$. The market price of risk associated with each of these factors are also computed by derived equations, following a transformation process; see Appendix E in [1].

The three-factors models are further calibrated using the ML-CCF which is an efficient estimator for the log-likelihood function. The calibration process identifies optimal parameters and the maximum of the log-likelihood function. In addition, the scores of likelihood functions are used to identify the market price of risk parameters; see [31,35]. The same assumptions of completely affine price of risk are made regarding market price of risk. Individual parameters ${\lambda }_{\upsilon }$, ${\lambda }_{\theta }$ and ${\lambda }_r$ are assumed to be easy to identify in the case where time-varying volatility exists. The scores of likelihood function are also applicable to identify the market price of risk as the parameters ${\lambda }_{\upsilon }$, ${\lambda }_{\theta }$ and ${\lambda }_r$ are expected to be non-constant nor collinear. From the invariant transformations, [1] derives the equations to compute these parameters in Appendix E.

10. Analysis of Results

10.1. Fitting the Yield Curve from the Instantaneous Rate

The main output from the three-factor models is the instantaneous rate $r_t$ from (22) and (25) for $A_1\left(3\right)$ and $A_2\left(3\right)$, respectively. The process also applies to their maximal counterparts. From the short rate $r_t$ we bootstrap zero-coupon bond prices using a discount function (1) which are thereafter inverted using (6) to obtain the model-implied yields. Performance analysis is based on the root mean square error (RMSE) and the residual ${\widehat{ϵ}}_t$, written as

$${\widehat{ϵ}}_{t,i}=y_{t,i}-Y_{t,i}$$

(38)

where $y_{t,i}$ and $Y_{t,i}$ are simulated yields and observed yields at time i for maturity t, respectively.

$$RMSE=\sqrt{{\displaystyle \frac{1}{NK}}{\widehat{ϵ}}^2}$$

(39)

where K= number of observations and N the number of maturities. Pairwise analysis between models is done using the above metrics.

10.2. In-Sample Analysis

10.2.1. $A_1\left(3\right)$ and $A{M}_1\left(3\right)$ Models

In-sample data were applied to implement both models through SDE’s, from which variables $\nu \left(t\right)$, $\theta \left(t\right)$ and $r\left(t\right)$ for stochastic volatility, central tendency and short rate, respectively were produced. Market price of risk parameters ${\lambda }_{\nu }$, ${\lambda }_{\theta }$ and ${\lambda }_r$, are also obtained from the results of these SDE according to the formulation from Appendix E in [1]. A linear combination of these variables result in a short rate that is fitted according to the affine function (22). The process is followed by implementation of optimisation techniques from which an efficient maximum likelihood and parameters estimates are obtained.

Table 4 presents the estimation results for both $A_1\left(3\right)$ model and its maximal counterpart $A{M}_1\left(3\right)$. Estimation results are based on the Fourier inversion of the characteristic function ML-CCF of the factor $Y\left(t\right)$ from (30). Variables $\nu \left(t\right)$, $\theta \left(t\right)$ and $r\left(t\right)$ were simulated from (23) where initial guess parameters, and the initial vector of factors $Y(i,t)$ for $i=1,3$, at time t, were inputs. According to [1] the BDFS model is a special case of (23) with the parameters in square brackets restricted to zero. However, in their specification, parameters ${\sigma }_{\theta r}$ and ${\sigma }_{rv}$ are relaxed and can be assumed take any other non-zero value.

The Nelder-Mead algorithm was used as it does not require a computation of gradients or second derivatives. This is despite its significant limitation that it cannot directly produce standard errors for the parameter estimates due to its inability to compute the Hessian matrix. To evaluate model goodness-of-fit and performance, we relied on the Chi-square (${\chi }^2$) statistics. Both AIC and BIC indicate how model complexity was penalised, taking into account both fit and the number of parameters.

AIC and BIC are both lower for $A{M}_1\left(3\right)$ compared to $A_1\left(3\right)$ indicating $A{M}_1\left(3\right)$ to be slightly better at balancing model complexity while it also has slightly more parameters. However, differences between AIC and BIC values, such as -1557.54 vs.-1615.92 for AIC, and -1460.68 vs. -1519.07 for BIC are quite small, suggesting that the models might be very similar in performance. Nonetheless, $A{M}_1\left(3\right)$ seems to be marginally better. The model with a lower ${\chi }^2$ of 16.38, which is the $A{M}_1\left(3\right)$ appears to have a better fit than the $A_1\left(3\right)$ with a ${\chi }^2$ of 24.98.

Parameters ${\kappa }_{r\nu }$, ${\sigma }_{\theta \nu }$, ${\alpha }_r$ and ${\beta }_{\theta }$ are restricted to zero for the $A_1\left(3\right)$ model, whereas they are relaxed for the maximal model $A{M}_1\left(3\right)$. In addition, some negative correlations were introduced for ${\sigma }_{r\theta }$ and ${\sigma }_{\theta r}$. A small movement in ${\kappa }_{r\nu }$ after calibration from zero to 0.035 for $A{M}_1\left(3\right)$ is exhibited. This suggest that the model has incorporated a weak positive relationship between short rate and volatility, which could reflect realistic market dynamics, though the effect is not large enough to dominate the model’s outcomes.

For model $A_1\left(3\right)$ there is no change in ${\sigma }_{\theta r}$ from -0.094 to -0.094 after calibration which is interpreted as an almost unchanged relationship between the volatility of the central tendency $\theta$ and short rate r. For $A{M}_1\left(3\right)$ a change to 0.089 reflects a slightly weak inverse relationship between the volatility of $\theta$ and r. $A{M}_1\left(3\right)$ is therefore slightly less sensitive to the impact of central tendency volatility on short rate volatility.

In the case of ${\sigma }_{r\theta }$, $A_1\left(3\right)$ there has been no change after calibration from-3.420 suggesting no increase in the negative relationship between the volatility of the short rate r and and the central tendency $\theta$. For $A{M}_1\left(3\right)$ a change from -3.420 to -3.770 indicates a significant increase in the negative correlation between the volatility of r and $\theta$. In summary, $A_1\left(3\right)$ exhibits a weaker relationship with volatilities while in the case of $A{M}_1\left(3\right)$ a slightly stronger negative correlation between short rate and central tendency is observed.

10.2.2. $A_2\left(3\right)$ and $A{M}_2\left(3\right)$ Models

Table 5 presents the results of parameter estimation for both models $A_2\left(3\right)$ and $A{M}_2\left(3\right)$. The second column list the initial guesses of parameters4. The third and fourth columns report the estimates of parameters that result from the optimisation.

Both AIC and BIC are almost the same, suggesting that either one of them is worth a selection from a model complexity perspective. However, individual differences of -1750.42 vs -1721.53 for AIC, and -1645.50 vs -1616.61 for BIC, suggest model $A{M}_2\left(3\right)$ to be slightly better at balancing model complexity while it also has slightly lower parameters. Marginally, the differences are too small implying both models to be balancing the parameter and model complexity in the same manner. ${\chi }^2$ statistics of 4.83 and 5.15 for models $A{M}_2\left(3\right)$ and $A_2\left(3\right)$, suggest a slightly better fit for the data$A_2\left(3\right)$. The same applies to the p-values of 1 for each model, which may suggest a good fit, except for the issue of possible overfitting.

Parameter ${\kappa }_{\theta \upsilon }$ under model $A_2\left(3\right)$ moves from -33.900 to -33.962 suggests a stronger mean reversion which may lead to lower stochastic volatility and a more stable behaviour of the short rate. For model $A{M}_2\left(3\right)$ it moves to -12.400, indicative of weaker mean reversion, higher stochastic volatility and more variability in the short rate. For $A_2\left(3\right)$, ${\kappa }_{r\nu }$ exhibit a movement from -35.300 to -35.071 (a difference of 0.229). This suggest an increase in the negative relationship between short rate and volatility, and an increase in volatility. A large negative value of -273.996 in ${\kappa }_{r\nu }$ for $A{M}_2\left(3\right)$ is observed, this implies that even a small change in the short rate would have a large effect on volatility.

The covariance parameters ${\sigma }_{r\upsilon }$ in $A_2\left(3\right)$ exhibits a move from -182.000 to -182.301 which is a modest increase in the negative relationship between the short rate r and volatility $\upsilon$. This implies a slight increase in negative correlation between the r and $\upsilon$. For $A{M}_2\left(3\right)$ a change from -182.000 to -133.003 indicates a significant increase in the negative correlation between the r and $\upsilon$. It is apparent that $A_2\left(3\right)$ exhibits a slight increase in relationship with volatilities while in the case of $A{M}_2\left(3\right)$ a slight increase in negative correlation between r and $\upsilon$ is observed.

Table 6 reports model performance in terms of both the RMSE and mean residuals. A pairwise analysis between model $A_1\left(3\right)$ and $A_1\left(3\right)$ is accomplished by comparing the RMSE, t-statistic and p-value for each model in the upper panel. The same analysis is carried out the mean residual ${\widehat{ϵ}}_{t,i}={ϵ}_{t,1}-{ϵ}_{t,2}$ where ${ϵ}_{t,1}$ and ${ϵ}_{t,2}$ are residuals for models $A_1\left(3\right)$ and $A{M}_1\left(3\right)$, respectively. The same variables can be applied to models $A_2\left(3\right)$ and $A{M}_2\left(3\right)$, whose values are sitting on the right hand side of the table.

Model $A_1\left(3\right)$ exhibits lesser RMSE than $A{M}_1\left(3\right)$ within the short end of the yield curve, up to the region just over 10-year maturity. For MAE model $A_1\left(3\right)$ is below $A{M}_1\left(3\right)$, although with a marginally small difference. This pattern is confirmed by Figure 4 (a) and (c). Statistically, despite one significance is noted along the 10-year maturity with a p-value of 0.012 which is lesser than 0.05, other results suggest no strong evidence that the models differ consistently in their performance in terms of RMSE. In the case of mean error $\widehat{{ϵ}_t}$, four instances exhibit negative direction implying lesser errors for model $A_1\left(3\right)$. Overall, $\widehat{{ϵ}_t}$ also do not prove a significant difference between the models despite four instances found to have p-values less than 0.05.

Model $A_2\left(3\right)$ exhibits lesser errors than $A{M}_2\left(3\right)$, both in terms of RMSE and MAE. The plots in Figure 4 (b) and (d) confirm the observation. There are three instances of weak statistical evidence of significant differences in RMSE between the models, but the majority of p-values are greater that 0.05 suggesting a strong evidence of a real difference. The models might perform similarly for those instances, while showing more distinct differences in some others. However, $\widehat{{ϵ}_t}$ has a p-values far less than 0.05 which confirms that models are significantly different. The plots are in favour of the differenes in models, particularly favouring model $A_2\left(3\right)$ versus its maximal counterpart, both in terms of RMSE and MAE.

Finally, we compare models $A_1\left(3\right)$ and $A_2\left(3\right)$ to determine the best performing in terms of the RMSE using a pairwise analysis. The results are reported in Table 7 where model $A_1\left(3\right)$ displays lower RMSE when compared to $A_2\left(3\right)$. Statistically, three instances show evidence for significance but the majority show a strong support for no differences in the models. The decision on model performance is in favour of $A_1\left(3\right)$.

10.3. Instanteneous Short Rate

Figure 5 plots the model-implied instantaneous rate against maturities. Both $A_1\left(3\right)$ and $A{M}_1\left(3\right)$ in blue and magenta, respectively. The plots are based on estimated parameters fitted to the models. $A_2\left(3\right)$ and its maximal counterpart are plotted against maturities in Figure 6.

The model-implied instantaneous short rate functions from (22) and (25) are made up of a linear combination of factors $(Y_1\left(t\right),Y_2\left(t\right),Y_3\left(t\right))$ mapped to $\upsilon \left(t\right)$, $\theta \left(t\right)$ and $r\left(t\right)$ and representing volatility, central tendency and short rate, respectively. They exhibit an average 11.9% and 8.1% short rates for $A_1\left(3\right)$ and $A{M}_1\left(3\right)$, respectively. The $A_1\left(3\right)$ and $A{M}_1\left(3\right)$ models exhibits higher volatility than the $A_2\left(3\right)$ and $A{M}_2\left(3\right)$ models. As the instantaneous rate evolves, mean-reversion persists throughout the process. Both subfamilies of models suggest that a positive but stable economic growth is forecasted, exept for a higher volatility as reflected in $A_1\left(3\right)$ and $A{M}_1\left(3\right)$ plot.

10.4. Out-Sample Analysis

Table 8 presents the results of RMSE and mean error $\widehat{{ϵ}_t}$ to determine both the gootness-of-fit and model performance. RMSE figures show smaller errors for $A_1\left(3\right)$ when compared to $A{M}_1\left(3\right)$. Statistically, p-values show a strong evidence of no differences in the models. $A_2\left(3\right)$ and $A{M}_2\left(3\right)$ on the other hand display a mixture of RMSE where early sections one model has lower RMSE and higher RMSE and in the other section, vice versa. The statistical results display a strong significance differences between the models $A_2\left(3\right)$ and $A{M}_2\left(3\right)$ with all p-values far below 0.05.

For mean error analysis $A_1\left(3\right)$ tends to be slightly more consistent with a mix of positive and negative $\widehat{{ϵ}_t}$, but $A{M}_1\left(3\right)$ shows slightly better performance in pairs where the differences are statistically significant. As result, $A{M}_1\left(3\right)$ may be considered to have slightly better performance in those instances. Model $A{2}_{(}3)$ and $A{M}_2\left(3\right)$ show a different patterns with $A{M}_2\left(3\right)$ performing slightly better thatn its non-maximal counterparts, but with a strong evidence for significance.

In Table 8 the final analysis shows that model $A_1\left(3\right)$ out-performs $A_2\left(3\right)$ when measured in terms of the RMSE.

11. Conclusions

The historical behaviour of the term structure of interest rates for the SA treasury bond was analysed. The purpose was to establish whether the ATSMs specifications of the three-factor models were suitable for the data. The study also considered the conditional volatilities and correlations of the state variables as they are essential components. It has therefore been essential to thoroughly examine the interactions among the state variables—stochastic volatility, central tendency, and the short rate—within the framework of a three-factor approach to the ATSM.

Results of both AIC, BIC for model selection and parameter complexity are in favour of maximal models. The ${\chi }^2$ statistic also show bias towards the maximal models. The differences in AIC, BIC and ${\chi }^2$ statistic show small differences marginally, suggesting that both models and their maximal counteparts may fit data equally well. Evidence from a pairwise analysis between models confirm this position with RMSE, mean errors and statistical results accordingly. In both cases, despite a strong possibility of model similarities, RMSE tends to exhibit slightly better performance for $A_1\left(3\right)$ and $A_2\left(3\right)$ than their maximal counterparts. We extended the tests to a pairwise analysis on both models, which lead to a conclusion that $A_1\left(3\right)$ out-performs $A_2\left(3\right)$ both in and out-sample.

Our study does not focus more on the market price of risk despite its importance. It is only relevant to preserve the formulation of the affine structure of the SDE and the transition from risk-neutral to physical measure. SA as an emerging market is susceptible to several risk factors, such as liquidity, exchange rate, credit risk and political. It is crucial that a specific form of market price of risk be incorporated into these specifications. We could not explore any of these forms of market price of risk that are mentioned within the ATSMs literature for SA. We leave them to future research.

The specification makes a selection of a number of factors M in the $A_M\left(N\right)$ model to serve as proxies for volatility. There are limitations in this method giving rise to a need to consider unspanned stochastic volatility models (USV). USV is capable of identifying latent volatility that could not be easily detected by three-factor models. They are also found to work better for capturing jumps, option pricing and hedging as they are able to highlight the hidden volatility. Future research on suitable specifications for SA should be extended to incorporating USV.