Framework for Simulation Study Involving Volatility Estimation: The GAS Approach

1. Introduction

In finance and econometrics, modelling of volatility has long been used to address a variety of issues concerning the risk and uncertainties of an asset [1,2]. This study is an extension of an earlier study on volatility estimation through the Generalized Autoregressive Conditional Heteroscedasticity (GARCH) model created by Bollerslev [3] (see Samuel et al. [4]). That is, the GARCH model was used in our earlier study to begin a "simulation framework for volatility estimation" work that was introduced in Samuel et al. [4] with the Monte Carlo simulation (MCS) design and guidelines. The Generalized Autoregressive Score (GAS) model is used in this instance to extend the robust design process. As a result, this study identifies a vital framework that is pertinent to the simulation of financial time series data for volatility estimation using the GAS model. The framework can be applied for volatility estimation in financial and non-financial areas.

The GAS is a dynamic model used for time-varying parameter modelling, much like the GARCH model (see [5]). Hence, this study’s main contribution is to show how the GAS model using the GAS package can be effectively applied to enhance time varying scale parameter (volatility) estimation through MC simulation, with empirical verification. Moreover, simulation through any model largely revolves around generating an appropriate dataset that can produce a reasonable output approximation. Thus, one further contribution of this study is to investigate how to use the GAS model to generate (simulate) a dataset that can adequately reflect the behaviour of financial time series data that is relevant for volatility estimation.

MCS studies are computer-based experiments that involve the use of known probability distributions for creating data by pseudo-random sampling. The data may be simulated through parametric models or repeated resampling [6]. MCS is used for approximating real-life scenarios, where it uses random quantities (or a probabilistic approach) to provide deterministic outcomes [7]. Since its introduction by three scientists in the field of physics [8], MCS has become an increasingly applied technique in a wide range of fields including finance, supply chain, engineering, science [9], and in all portfolio and investment types [10]. One of the reasons for its wide applications is because it is based on the concept of the central limit theorem where the distribution of its estimates follows (or converges to) the Normal distribution [11].

The raw data used for this study are the daily closing Standard & Poor South African sovereign bond index, abbreviated S&P SA bond index. They are the bond market’s data in US Dollars from Datastream [12] for the period 4th January 2000 to 11th June 2021 with 5598 observations. The period under study includes the covid-19 pandemic crisis period that began in 2020. The rest of this paper is structured as follows: Section 2 reviews the theories underlying the GAS model, the true parameter recovery (TPR) measure, and the description of the simulation design. Section 3 shows how the simulation framework is illustrated using financial bond return data, with empirical verifications. Section 4 discusses the main findings and Section 5 presents concluding remarks on the novelty of the study.

2. Materials and Methods

2.1. The GAS Model

Volatility can be effectively estimated through a Score Driven (SD) model known as the GAS introduced by Harvey [13] and Creal et al. [14] (see [5,15,16]). The GAS model is also called the Dynamic Conditional Score (DCS) model, and its applications are used in many financial econometrics analytics (see [17]). It provides a connection between the stochastic volatility (SV) models [18] and GARCH models (see [19]).

The GAS model uses the score of the conditional density function to determine the time-variation in the parameters, i.e., the use of the score function is applied as the driver of the parameters’ time variation. Also, the score functions are robust to outliers because they discount extreme values by reducing the number of extreme observations, which takes place because the GAS model is conditional on time t [13]. For instance, if errors are modelled using fat-tailed innovations under a GAS model, a big absolute return swing will not lead to a large increase in volatility [14]. Hence, the GAS model is indeed robust and quite suitable for modelling skewed and/or fat-tailed time series data like financial returns [13,20]. Moreover, like the other observation-driven models, likelihood estimation using the GAS model is simple and direct [5]. Observation-driven modelling exist in the presence of time-varying parameters, where parameters are functions of lagged dependent variables, concurrent variables and lagged exogenous variables (see [14,21] for details).

Let a random vector of M-dimension at time t with a conditional distribution be represented by $r_t$∈${\Re }^M$ such that:

$$r_t|r_{1:t-1}\sim f(r_t;{\vartheta }_t),$$

(1)

where $r_{1:t-1}$≡${(r_1^{\top },...,r_{t-1}^{\top })}^{\top }$ are the past time-values of $r_t$ up till $t-1$, that indicate the set of information available prior to time t [22], and ${\vartheta }_t$∈$\mathsf{\Theta }$⊆${\Re }^I$ is a time-varying parameters’ vector with complete characterization of $f(·)$ depending only on $r_{1:t-1}$ and a group of additional time-invariant parameters $\zeta$, such that ${\vartheta }_t\equiv \vartheta (r_{1:t-1},\zeta )\forall t$. The key attribute of GAS models is that the score of the conditional distribution (i.e., the conditional probability density function of $r_t$, as stated in Equation 1) drives the evolution in the time-varying parameter vector ${\vartheta }_t$, in relation to an autoregressive element:

$${\vartheta }_{t+1}\equiv \kappa +\mathbf{A}{c}_t+\mathbf{B}{\vartheta }_t,$$

(2)

where $c_t$ is the innovation term denoting a possibly scaled score of $f(r_t;{\vartheta }_t)$ with respect to ${\vartheta }_t$, estimated in $r_t$ (see [14,15,23]). The intercept $\kappa$, vectors $\mathbf{A}$ and $\mathbf{B}$ are proper dimensional matrices of coefficients that control for the time-varying parameters’ $\left({\vartheta }_t\right)$ evolution [5,15]. These are functions of the static parameters vector $\zeta$, that is $\zeta$ = ($\kappa$, $\mathbf{A}$, ${\mathbf{B})}^{\top }$ [14,17,23]. Various classes of observation-driven models can be obtained with different choices of the scaled score $c_t$, hence these types of models are referred to as score-driven models [22]. The autoregressive Equation (2) that acts as a mechanism for updating ${\vartheta }_t$ can be generally stated (see [14]) as:

$${\vartheta }_{t+1}=\kappa +\sum _{i=1}^p{A}_i{c}_{t-i+1}+\sum _{j=1}^q{B}_j{\vartheta }_{t-j+1}.$$

(3)

The $c_t$ is a vector relative to the score of Equation (1):

$$c_t\equiv C_t\left({\vartheta }_t\right){\nabla }_t(r_t,{\vartheta }_t).$$

(4)

For any realized observation $r_t$, the time-varying ${\vartheta }_t$ is updated to the next period $t+1$ using Equation (3) with $c_t$ as defined in Equations (4) and (5) [14]. The matrix $C_t$ denotes an I × I (valued) positive definite scaling matrix identified at time t, and:

$${\nabla }_t(r_t,{\vartheta }_t)\equiv \frac{\partial logf(r_t;{\vartheta }_t)}{\partial {\vartheta }_t}$$

(5)

represents the score of Equation (1) estimated at ${\vartheta }_t$. The variance of ${\nabla }_t$ can be accounted for by setting the $C_t$ (scaling matrix) to a power $\omega >0$ of the inverse of the Information Matrix, $J_t$, of the time-varying parameters ${\vartheta }_t$ [5,14], i.e.,

$$C_t\left({\vartheta }_t\right)\equiv J_t{\left({\vartheta }_t\right)}^{-\omega },for:J_t\left({\vartheta }_t\right)\equiv {\mathtt{E}}_{t-1}\left[{\nabla }_t(r_t,{\vartheta }_t){\nabla }_t{(r_t,{\vartheta }_t)}^{\top }\right],$$

(6)

with the expectation taken on $r_t|r_{1:t-1}$. For $\omega$, values can be taken in the set $\{0,\frac{1}{2},1\}$. If $\omega$ = $\frac{1}{2}$ or $\omega$ = 1 is fixed by a user, the inverse of the square root of the conditional score ${\nabla }_t(r_t,{\vartheta }_t$) pre-multiplies the conditional score’s covariance matrix $J_t\left({\vartheta }_t\right)$. However, when $\omega$ = 0, there will be no scaling since $C_t$ = $\mathbf{I}$ (identity matrix). It should be noted that for any choice of $\omega$, $c_t$ follows a Martingale Difference with respect to the distribution of $r_t|r_{1:t-1}$ since the score’s expectation is zero, that is, E${}_{t-1}$ [$c_t$] = 0 ∀t. More so, we can easily derive the additional moment condition ${\vee }_{t-1}$ [$c_t$] = $\mathbf{I}$ when $\omega$ = $\frac{1}{2}$ [5]. Because $c_t$ is a Martingale Difference, the time-varying parameters’ process can revert to the mean near its long-term value of mean ${(\mathbf{I}-\mathbf{B})}^{-1}\kappa$ once the spectral radius of $\mathbf{B}<1$. This means that the I-valued vector $\kappa$ controls for the level, while the I × I matrix $\mathbf{B}$ controls for the persistence of the process [5,22]. In Equation (2), the extra $I\times I$ matrix of $\mathbf{A}$’s coefficients of the scaled score $c_t$ controls the effect of $c_t$ on ${\vartheta }_{t+1}$, where $c_t$ indicates the direction for updating the parameters vector from ${\vartheta }_t$ to ${\vartheta }_{t+1}$ [5,14].

2.1.1. Reparameterisation of the Time-Varying Parameter

The parameter vector ${\vartheta }_t$ in Equation (2) is unbounded because it has a linear specification, and its parameter space is regularly restricted (i.e., $\mathsf{\Theta }$⊂${\Re }^I$) in practice. Therefore, the GAS structure uses a nonlinear link function $\mathsf{\Lambda }(·)$, as an alternative to the positivity imposing constraints, to map the reparametrised parameters vector $\breve{{\vartheta }_t}$∈${\Re }^I$ into ${\vartheta }_t$ such that $\breve{{\vartheta }_t}$∈${\Re }^I$ also provides a linear effective specification experienced in Equation (2). Based on this, the ${\vartheta }_t$ updating equation now takes the GAS recursion of the form:

$$\begin{array}{ccc}\hfill {\vartheta }_t& \equiv & \mathsf{\Lambda }\left(\breve{{\vartheta }_t}\right)\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \breve{{\vartheta }_t}& \equiv & \kappa +\mathbf{A}{\breve{c}}_{t-1}+\mathbf{B}{\breve{\vartheta }}_{t-1},\hfill \end{array}$$

(8)

where $\kappa$≡${({\kappa }_{\mu },{\kappa }_{\sigma },{\kappa }_{\nu })}^{\top }$ represents the vkappa vector, $\mathbf{A}\equiv \mathrm{diag}(a_{\mu },a_{\sigma },a_{\nu })$, $\mathbf{B}\equiv \mathrm{diag}(b_{\mu },b_{\sigma },b_{\nu })$ and $\mathsf{\Lambda }(·)$ is a mapping function [5,16]. The estimated coefficients {${\kappa }_{\mu }$, ${\kappa }_{\sigma }$, ${\kappa }_{\nu }\}$, {$a_{\mu }$, $a_{\sigma }$, $a_{\nu }$} and {$b_{\mu }$, $b_{\sigma }$, $b_{\nu }$} are the elements of vectors $\kappa$, $\mathbf{A}$ and $\mathbf{B}$ respectively, where $\mu$, $\sigma$, and $\nu$ are the location, scale and shape parameters in that order. The initialisation of the process is at ${\vartheta }_1={(\mathbf{I}-\mathbf{B})}^{-1}\kappa$, and ${\breve{c}}_t\equiv {\breve{C}}_t\left({\breve{\vartheta }}_t\right){\breve{\nabla }}_t(r_t,{\breve{\vartheta }}_t)$, where ${\breve{\nabla }}_t(r_t,{\breve{\vartheta }}_t)$ is the score of Equation (1) with respect to ${\breve{\vartheta }}_t$ (further details on this can be found in [5,16]).

Based on the dependence of the driving mechanism in Equation (3) on the scaled score vector in Equation (4), Equations (1) to (5) describe the GAS model of p and q as GAS($p,q$), and the candidate model GAS$(1,1)$ is used for the time-varying scale parameter, or volatility, estimation in this study. This is because the time varying parameter ${\vartheta }_t$ evolves in the form of an autoregressive (AR) process of order 1 [23], and the model typically features as GAS$(1,1)$ [24]. Hence, GAS$(1,1)$ model is used as the data generating process (DGP) for the simulation modelling, and it also serves as the only candidate model for the empirical data modelling. The volatility persistence of the GAS$(1,1)$ model is defined as $B_1$ (see [5,14]). Volatility persistence is a measure of how fast the existing shock to volatility will fade away [25] and revert to the mean (normal) state.

The maximum likelihood estimation (MLE) of $\widehat{\vartheta }$ can be implemented through the prediction error decomposition, given the static parameter vector $\zeta$ and the past information $r_{1:t-1}$. That is, for a sample of N realisations of $r_t$, gathered in $r_{1:N}$, the MLE of the vector of parameters $\zeta$ is a solution of:

$$\widehat{\zeta }\equiv \underset{\zeta }{\arg \max }L(\zeta ;r_{1:N}),\mathrm{where}L(\zeta ;r_{1:N})\equiv \log f(r_1;{\vartheta }_1)+\sum _{t=2}^N\log f(r_t;{\vartheta }_t)$$

(9)

with ${\vartheta }_1\equiv {(\mathbf{I}-\mathbf{B})}^{-1}\kappa$, and the time-varying parameters’ vector ${\vartheta }_t\equiv \vartheta (r_{1:t-1},\zeta )$ for $t>1$, such that ${\vartheta }_t$ depends on $r_{1:t-1}$ and $\zeta$.

2.2. The True Parameter Recovery Measure

The true parameter recovery (TPR) measure in Equation (10) was developed by Samuel et al. [4] as a means of measuring the performance of the MCS estimator in recovering the true parameter. It is used as a proxy for the coverage of the MCS experiment to determine how much of the true parameter is recovered by the MC simulation estimator.

$$\mathrm{TPR}=\left(K-\left[\frac{(\vartheta -\widehat{\vartheta })}{\vartheta }\times K\right]\right)\%,$$

(10)

where $\vartheta$ is the true data-generating parameter, K = 0, 1, 2, …, 100 is the nominal recovery level, and $\widehat{\vartheta }$ is the estimator from the simulated data. For instance, a TPR estimated value of 95% or 100% indicates that the MCS estimator recovers 95% or 100% of the true parameter. A full recovery of the true parameter $\vartheta$ can be achieved by the MCS estimator $\widehat{\vartheta }$ when $\widehat{\vartheta }$ = $\vartheta$, where $\vartheta >0$, such that the TPR output equals the given nominal recovery level K (i.e., TPR = K%).

2.3. Simulation Design

The design of the simulation framework involves "background (optional), defining the aim, research questions, method of implementation, and summarised conclusion". The method of implementation is a workflow that consists of writing the code, setting the seed, setting the true parameter a priori, data generation process, and performance evaluation through meta-statistics (see [4]). These steps, as summarized by the flowchart in Figure 1, are relevant for successful simulations through the GAS model. The details on each design step are as follows:

2.3.1. Aim of the Simulation Study

After optionally stating the background of the study, the next step is to define the aim of the study and it must be concisely, clearly and unambiguously stated for reader’s understanding. The focus of MCS studies generally dwells on estimators’ capabilities in recovering the true parameters $\vartheta$, such that E$\left(\widehat{\vartheta }\right)$ = $\vartheta$ for unbiasedness, $\widehat{\vartheta }\to \vartheta$ as the sample size $N\to \infty$ for consistency, and root mean square error (RMSE) or standard error (SE) tends to zero as $N\to \infty$ for good efficiency or precision of the true parameter’s estimator. Hence, the aim of the study may revolve around those properties, like the bias or unbiasedness, consistency, efficiency or precision of the estimator. The aim can also evolve from comparisons of multiple entities, like comparing the efficiency of various error distributions, or comparing the performance of multiple models or modification of an existing method.

2.3.2. State the Research Questions

After defining the study aim, relevant questions with regards to the purpose of the simulation should be outlined. These will be pointers to the objectives of the study. The intricacies of some statistical research questions make them better resolved via simulation approaches. Simulation provides a robust procedure for responding to a wide range of theoretical and methodological questions and can offer a flexible structure for answering specific questions pertinent to researchers’ quest [26].

2.3.3. Method of Implementation

The simulation and empirical modelling of this study are implemented in R Statistical Software, version 4.0.3, with RStudio version 2023.03.1+446, using the GAS [5,16], SimDesign [27], tidyverse [28], patchwork [29], aTSA [30], rugarch [31,32], and forecast [33] packages. Computation is executed on Intel(R) Core(TM)i5-8265U CPU @ 1.60GHz 1.80 GHz. The method of implementing the simulation is as follows:

• Write the code: Carrying out a proper simulation experiment that mirrors real-life situation can be very demanding and computationally intensive, hence readable computer code with the right syntax must be ensued. MCS code must be well organised to avoid difficulties during debugging. It is always safer to start with small coding practices, get familiar with them and ensure they run properly with necessary debugging of errors before embarking on more intensive and complex ones. Code must be efficiently and flexibly written and well arranged for easy readability.
• Set the seed: Simulation code will generate different sequence of random numbers each time it is run unless a seed is set [34]. A set seed initialises the random number generator [32] and ensures reproducibility, where the same result is obtained for different runs of the simulation process [35]. The seed needs to be set only once, for each simulation, at the start of the simulation session [6,32], and it is better to use the same seed values throughout the process [6]. The seed is essentially used for reproducibility. Samuel et al. [4] used the GARCH model to show that as the sample size N increases, the arrangement (or pattern) of the seed values generally does not influence the efficiency and consistency of an estimator. This however may depend on the quality of the model used.
• After setting the seed, the true parameter value/s of the true sampling distribution (or true model) are then set a priori [22,36].
• Next, simulated observations are generated using the true sampling distribution or the true model given some sets of (or different sets of) fixed parameters. Generation of simulated datasets through the GAS model is carried out using the R package "GAS". Random data generation involving this package can be implemented using either of two approaches [5,16]. The first approach is to carry out the data generating simulation directly on a fitted object "fit" (or uGASFit) using the UniGASSim (univariate GAS simulation) function for the simulated random data. The second approach uses the function UniGASSim with "fit = null". This latter approach involves full specification of a GAS model that includes selection of the conditional error distribution of the time series process, and specifying the static parameter $\zeta$ = ($\kappa$, $\mathbf{A}$, $\mathbf{B}$) that controls the time variation in ${\vartheta }_t$, as described in Section 2.1.
  The simulation or data generating process can be run once or replicated multiple times. However, data generation through the GAS process is generally designed to run only once. In their study through autoregressive process, Samuel et al. [4] showed that the outcomes of MCS experiments using the GARCH model are the same for datasets simulated with the same seed value, regardless of whether the simulation or data generating process is run once or replicated multiple times. Hence, the outcome of a single run is the same as the average of multiple runs or replicates. Their findings further revealed that in any multiple generated (simulated) datasets in the family GARCH model, each simulated series is distinct in randomness and shape, and it is different from every other series in the datasets.
• The generated (simulated) data are analysed and the estimates from them are assessed using classic methods through meta-statistics to derive relevant information about the estimators. Meta-statistics (see [27]) are performance measures or metrics for evaluating the modelling outcomes by judging the closeness between an estimate and the true parameter. A few of the frequently used meta-statistical summaries, as described below, include the bias, root mean square error (RMSE), standard error (SE) and relative efficiency (RE). For more meta-statistics, see [6,27,37].

Bias

The bias, on average, measures the tendency of the simulated estimators $\widehat{\vartheta }$ to be smaller or larger than their true parameter value $\vartheta$. It is defined as the average difference between the true (population) parameter and its estimate [38]. The optimal value of bias is 0 [37,39]. An unbiased estimator, on average, yields the correct value of the true parameter. Bias with a positive (negative) value indicates that the true parameter value is over-estimated (under-estimated). However, in absolute values, the closer the estimator is to 0, the better it is. Bias is stated mathematically as $\mathtt{E}(\tilde{\vartheta }-\vartheta )$, but can be presented in Monte Carlo simulations (see [40]) as $\mathrm{bias}=\frac{1}{L}{\sum }_{i=1}^L(\widehat{{\vartheta }_i}-\vartheta )$. The two formulae are connected as:

$$\mathtt{E}(\tilde{\vartheta }-\vartheta )=\mathtt{E}\left(\tilde{\vartheta }\right)-\vartheta =\frac{1}{L}\sum _{i=1}^L{\widehat{\vartheta }}_i-\vartheta =\frac{1}{L}\sum _{i=1}^L({\widehat{\vartheta }}_i-\vartheta ),$$

(11)

where $\tilde{\vartheta }$ = ${\widehat{\vartheta }}_i(i=1,...,L)$ is a finite $i^{th}$ number of the sample estimate $\widehat{\vartheta }$ for the datasets, L is the number of replications, and $\vartheta$ is the true parameter.

Standard Error

Sampling variability in the estimation can be evaluated via the standard error (SE) as stated (see [40,41]) in Equation (12). Also called Monte Carlo standard deviation, it is a measure of the efficiency or precision of the true parameter’s estimator, where it is used to estimate the long-run standard deviation of ${\widehat{\vartheta }}_i$ for finite repetitions. It does not require knowing the true parameter $\vartheta$ but depends on its estimator ${\widehat{\vartheta }}_i$ only. The smaller the sampling variability, the more the efficiency or precision of the $\vartheta$’s estimator (see [6]). Sampling variability decreases with increased sample size [37].

$$\mathrm{SE}=\sqrt{\frac{1}{L}\sum _{i=1}^L{({\widehat{\vartheta }}_i-\overline{\vartheta })}^2},\mathrm{where}\overline{\vartheta }=\sum _{i=1}^L\widehat{\vartheta }/L.$$

(12)

RMSE

The root mean square error (RMSE) is an accuracy measure for evaluating the difference between a model’s true value and its prediction. RMSE measure indicates the sampling error of an estimator when compared to the true parameter value [37] and it is stated as:

$$\mathrm{RMSE}=\sqrt{\frac{1}{L}\sum _{i=1}^L{\left({\widehat{\vartheta }}_i-\vartheta \right)}^2}.$$

(13)

Its computation involves the true parameter $\vartheta$. An estimator with smaller RMSE is more efficient in recovering the true parameter value [37,41], and minimum RMSE produces maximum precision [42]. Consistency of the estimator occurs when RMSE decreases such that $\widehat{\vartheta }\to \vartheta$ as the sample size $N\to \infty$ [6,31]. RMSE is related to bias and sampling variability as:

$$\mathrm{RMSE}=\sqrt{{\mathrm{bias}}^2+{\mathrm{SE}}^2}.$$

(14)

That is, the RMSE is an inclusive measure that combines bias and SE, such that low SE can be penalised for bias. The mean squared error (MSE) is obtained by squaring the RMSE. MCS is highly dependent on the law of large numbers, and it is expected that the distribution of a suitably large sample should converge to that of the underlying population as the sample size increases [43]. It is also expected that increasing the sample size should relatively decrease the sampling error, but this is not always the case. That is, the number of sample size cannot always be sufficiently increased to limit the sampling error to a bearable level [43]. Using sample sizes between 100 and 7000 inclusive, Kwiatkowski [11] graphically shows the effect of varying sample sizes on simulation outcomes. The author reveals that as the sample size increases, the spread of an estimated value decreases proportionately. However, at some point, additional sample increase does not further reduce the spread [11].

Relative Efficiency

Sampling efficiency between one estimator and the other can be compared using the relative efficiency (RE) statistic by creating ratios between the estimators’ RMSE values as shown in Equation (15). In the formula, the first RMSE or RMSE${}_1$ is taken as the reference estimator, such that other RMSEs are assessed in relation to it. This performance measure can be used when the estimators of interest are unbiased. The RE ≡ 1 when i = 1, and values above (less than) 1 indicate less (greater) efficiency. Less efficiency means more variability than the statistic that is referenced, while greater efficiency indicates less variability and higher precision for recovering the true parameter $\vartheta$ [37,40]. For instance, $RE<1$ suggests that estimator i is the preferred one, i.e., if i = 2 for instance, estimator 1 is less efficient in relation to estimator i = 2 [44].

$$\mathrm{RE}={\left(\frac{{\mathrm{RMSE}}_i}{{\mathrm{RMSE}}_1}\right)}^2.$$

(15)

2.4. Discussion and Summary

After the method of implementation, the last stage in the framework steps is the conclusion, which should reflect a summary discussion of all logical findings from the experiments, with answers to the research questions. The conclusion brings out the novelty of the research and may also include the limitations experienced and opportunities for future work. In addition, relevant information on simulation results can be conveyed through graphics, tabular presentation, or both.