Framework for Simulation Study Involving Volatility Estimation: The GARCH Approach

1. Introduction

A simulation-based experiment is not often included in research because many upcoming researchers do not have an adequate understanding of the nitty-gritty involved. Although the details involved in simulation modelling are generally inexhaustible, this study, however, unveils a crucial framework relevant for the simulation of financial time series data using the Generalised Autoregressive Conditional Heteroscedasticity (GARCH) model for volatility estimation. Volatility is a measure of the variability in an asset’s price over time [1]. The ultimate goal of the study is to familiarise researchers with the concepts and art of simulation modelling through this model. This framework is amply general for broad applicability and can be more easily verified in typical situations in other non-financial sectors concerned with volatility estimation. The framework utilises the robust simulating resources of the GARCH model, through set parameters, to generate data that are analysed, and the estimates from the process are then used by chosen metrics to explain the behaviour of selected statistics of interest.

Monte Carlo simulation (MCS) studies are computer-based experiments that use known probability distributions to create data by pseudo-random sampling. The data may be simulated through a parametric model or via repeated resampling [2]. MCS applies the concept of imitating a real-life scenario on the computer through a certain model that can hypothetically generate the scenario. By simulating or repeating this process a considerably large number of times, it is possible to obtain outcomes that can enable precise computation of desired issues of concern, like the possible assumed error distribution/s that can suitably describe a given stock market for volatility estimation. Simulation approaches can be applied flexibly to numerous problems and used to derive appropriate solutions to questions that may be non-derivable by other methods [3]. In preparing for a simulation experiment, reasonably ample time is needed to organise a well-written and readable computer code and for simulated data generation. Implementation of a good simulation experiment and reporting the outcomes require adequate planning. There is currently no general one-size-fits-all approach to simulation modelling and to the choice of an adequate number of sample size/s, but the process can be narrowed down to individual models.

Series of R application software packages like the rugarch [4], GAS [5], SimDesign [6], tidyverse [7], to mention but a few, are currently available for simulation studies. This study exemplifies how the GARCH model through the rugarch package can be effectively used to improve volatility estimation through MCS experiment, with outcomes verified by real data empirical modelling. Although there are good books on simulation approaches in general (see [8,9,10]), but up until now, to the best of our knowledge, there has not been any monograph with direct step by step comprehensive layout on simulation framework using the GARCH model. Hence, this study rolls out an inclusive simulation design that is summarily required for a robust simulation practice in finance to estimate volatility using this model, and the knowledge can be applied in any other field. Since the rugarch package does not make provision for calculating the coverage probability 1, this study also computes the MCS estimator’s recovery levels through the "true parameter recovery (TPR)" measure as a proxy for the coverage. The results show that the MCS estimates considerably recover the true parameters.

The raw data used for this study are the daily closing S&P South African sovereign bond index, abbreviated S&P SA bond index. They are Standard & Poor data for the bond market in US Dollars from Datastream [12] for the period 4th January 2000 to 11th June 2021 with 5598 observations. The rest of the paper is organised as follows: Section 2 reviews the theories underpinning two heteroscedastic models, the TPR measure, and the description of the design of the simulation framework. Section 3 presents the practical illustration of the simulation framework, with empirical verification, on financial bond return data. Section 4 discusses the key findings and Section 5 concludes.

2. Materials and Methods

2.1. The GARCH Model

The GARCH model was developed by Bollerslev [13] as a generalisation of the Autoregressive Conditional Heteroscedasticity (ARCH) model introduced by Engle [14]. It is a classical model that is normally defined by its conditional mean and variance equations for modelling financial returns volatility [15]. The mean equation is stated as:

$$r_t={\mu }_t+{\epsilon }_t,$$

(1)

where $r_t$ is the return series, ${\epsilon }_t=z_t{\sigma }_t$ denotes the residual part of the return series that is random and unpredictable, where $z_t\sim N(0,1)$ are the standardised residuals which are independent and identically distributed (i.i.d.) random variables with mean 0 and variance 1 [16,17], ${\mu }_t$ is the mean function that is usually stated as an Autoregressive Moving Average (ARMA) process,

$${\mu }_t=\sum _{i=1}^p{\varphi }_i{r}_{t-i}+\sum _{i=1}^q{\theta }_i{\epsilon }_{t-i},$$

(2)

where ${\varphi }_i$ (i = 1, …, p) and ${\theta }_i$ (i = 1, …, q) are unknown parameters. The variance equation of the GARCH$(u,v)$ model is defined as:

$${\sigma }_t^2=\omega +{\alpha }_1{\epsilon }_{t-1}^2+\cdots +{\alpha }_v{\epsilon }_{t-v}^2+{\beta }_1{\sigma }_{t-1}^2+\cdots +{\beta }_u{\sigma }_{t-u}^2,$$

(3)

where $\omega >0$ is the intercept (white noise), coefficients ${\alpha }_j\ge 0$ (j = 1, …, v) and ${\beta }_i\ge 0$ (i = 1, …, u) respectively measure the short-term and long-term effects of ${\epsilon }_t$ on the conditional variance [18]. The non-negativity restrictions on the unknown parameters, ${\alpha }_j$ and ${\beta }_i$, are imposed for ${\sigma }_t^2>0$. The equation shows that the conditional variance ${\sigma }_t^2$ is a linear function of past squared innovation ${\epsilon }_{t-j}^2$ and past conditional variances ${\sigma }_{t-i}^2$. The GARCH model is a more parsimonious specification [19,20] since it is an equivalence of a certain ARCH(∞) model [21]. When $u=0$ in Equation (3), the GARCH model changes to the ARCH model with conditional variance stated as:

$${\sigma }_t^2=\omega +{\alpha }_1{\epsilon }_{t-1}^2+\cdots +{\alpha }_v{\epsilon }_{t-v}^2.$$

(4)

GARCH(1,1) is the simplest model specification with u = 1 and v = 1 in Equation 3, and it is conceivably the best candidate GARCH model for several applications [21,22]. The volatility persistence of the GARCH$(1,1)$ is defined as ${\alpha }_1$ + ${\beta }_1$ (see [23,24]). Volatility persistence is used to evaluate the speed of decay of shocks to volatility [15]. Volatility exhibits long persistence into the future if $\alpha$ + $\beta$→ 1, hence the closer the sum of the coefficients is to one (zero), the greater (lesser) the persistence. However, if the sum is equal to one, then shocks to volatility persist forever and the unconditional variance is not determined by the model. This process is called "Integrated-GARCH" [23,25]. If the sum is greater than one, the conditional variance process is explosive, suggesting that shocks to the conditional variance are highly persistent. Covariance stationarity of the GARCH model is ensued when ${\sum }_{j=1}^v{\alpha }_j+{\sum }_{i=1}^u{\beta }_i<1$, while the unconditional variance of ${\epsilon }_t$ is ${\sigma }^2=\mathtt{E}\left({\epsilon }_t^2\right)=\omega /\{1-({\sum }_{j=1}^v{\alpha }_j+{\sum }_{i=1}^u{\beta }_i)\}$ [21].

For the maximum likelihood estimation (MLE), the log-likelihood function for maximising the likelihood of the unknown parameters given the observations is stated as:

$$L\left(\vartheta \right|\epsilon )=\sum _{t=1}^N\frac{1}{\sqrt{2\pi {\sigma }_t^2}}\exp \left(-\frac{{\epsilon }_t^2}{2{\sigma }_t^2}\right),$$

(5)

where $\vartheta ={(\mu ,\omega ,{\alpha }_1,\dots ,{\alpha }_v,{\beta }_1,\dots ,{\beta }_u)}^{\top }$ is a vector of parameters, and $\epsilon =({\epsilon }_1,\dots ,{\epsilon }_N)$ is a realisation of length N. The quasi-maximum likelihood estimation (QMLE) based on the Normal distribution and MLE have the same set of instructions for estimating $\widehat{\vartheta }$, the only difference, however, is in the estimation of a robust standard deviation of $\widehat{\vartheta }$ (see [26,27,28]).

The maximised log-likelihood function with Student’s t distribution (Duda and Schmidt, 2009) is stated as:

$$lnL\left(\vartheta \right|\epsilon )=\sum _{t=1}^N\left(ln\left[\frac{\Gamma \left(\frac{\nu +1}{2}\right)}{\Gamma \left(\frac{\nu }{2}\right)\sqrt{(\nu -2)\pi }}\right]-\frac{1}{2}ln\left({\sigma }^2\right)-\left(\frac{\nu +1}{2}\right)ln\left[1+\frac{{\epsilon }_t^2}{{\sigma }^2(\nu -2)}\right]\right),$$

(6)

where $\Gamma (·)$ and $\nu$ are the gamma function and degree of freedom, respectively.

2.2. The fGARCH Model

The family GARCH (fGARCH) model, developed by Hentschel [29], is an inclusive model that nests some important symmetric and asymmetric GARCH models as sub-models. The nesting includes the simple GARCH (sGARCH) model [13], the Absolute Value GARCH (AVGARCH) model [30,31], the GJR GARCH (GJRGARCH) model [32], the Threshold GARCH (TGARCH) model [33], the Nonlinear ARCH (NGARCH) model [34], the Nonlinear Asymmetric GARCH (NAGARCH) model [35], the Exponential GARCH (EGARCH) model [19], and the Asymmetric Power ARCH (apARCH) model [36]. The sub-model apARCH is also a family model (but less general than the fGARCH model) that nests the sGARCH, AVGARCH, GJRGARCH, TGARCH, NGARCH models, and the Log ARCH model [37,38]. The fGARCH$(u,v)$ model is stated as:

$${\sigma }_t^{\gamma }=\omega +\sum _{j=1}^v{\alpha }_j{\sigma }_{t-j}^{\gamma }\left(\right|z_{t-j}-{\lambda }_{2j}{|-{\lambda }_{1j}\{z_{t-j}-{\lambda }_{2j}\})}^{\delta }+\sum _{j=1}^u{\beta }_j{\sigma }_{t-j}^{\gamma }.$$

(7)

This robust fGARCH model allows different powers for ${\sigma }_t$ and $z_t$ to drive how the residuals are decomposed in the conditional variance equation. Equation (7) is the conditional standard deviation’s Box-Cox transformation, where the transformation of the absolute value function is carried out by the parameter $\delta$, and $\gamma$ determines the shape. The ${\lambda }_{2j}$ and ${\lambda }_{1j}$ control the shifts for asymmetric small shocks and rotations for large shocks, respectively. The fit of the full fGARCH model can be implemented with $\gamma$ = $\delta$ (see [24]). Volatility clustering in the returns can be quantified through the model’s volatility persistence stated as:

$$\widehat{P}=\sum _{j=1}^u{\beta }_j+\sum _{j=1}^v{\alpha }_j{\varrho }_j,$$

(8)

where ${\varrho }_j$, expressed in Equation (9), is the expected value of $z_t$ in the absolute value asymmetry term’s Box-Cox transformation. Volatility clustering implies that large changes in returns tend to be followed by large changes and small changes tend to be followed by small changes. The persistence is obtained in this study through the "persistence()" function in R rugarch package. See [4,29] for details on fGARCH and the nested models.

$${\displaystyle {\varrho }_j=\mathtt{E}\left(\right|z_{t-j}-{\lambda }_{2j}|-{\lambda }_{1j}(z_{t-j}-{\lambda }_{2j}){)}^{\delta }={\int }_{-\infty }^{\infty }\left(\right|z-{\lambda }_{2j}{|-{\lambda }_{1j}(z-{\lambda }_{2j}))}^{\delta }f(z,0,1,\dots )dz}$$

(9)

2.3. The True Parameter Recovery Measure

Since the focus of MCS studies involves the ability of the estimator to recover the true parameter (see [39]), this study applies the "true parameter recovery (TPR)" measure in Equation (10) to compute the level (degree) of recovery of the true parameter through the MCS estimator. The TPR measure is a means of evaluating the performance of the MCS estimates in recovering the true parameter. That is, it is used to determine how much of the true parameter value is recovered by the MCS estimator.

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

(10)

where K = 0, 1, 2, …, 100 is the nominal recovery level, $\vartheta$ is the true data-generating parameter and $\widehat{\vartheta }$ is the estimator from the simulated (synthetic) data. For instance, a TPR estimated value of 95% or 100% denotes that the MCS estimator recovers the complete 95% or 100% of the true parameter.

2.4. Simulation Design

The design of the simulation framework includes "background (optional), defining the aim, research questions, method of implementation, and summarised conclusion". The method of implementation is a workflow that involves writing the code, setting the seed, setting the true parameter/s a priori, data generation process, and performance evaluation through meta-statistics. As summarised by the flowchart in Figure 1, these crucial steps are relevant for successful simulations through the GARCH model. The details of each design step are as follows:

2.4.1. Aim of the Simulation Study

After optionally stating the background that explains crucial underlying facts about the study, the next step is to define the aim of the study, and it must be clearly, concisely and unambiguously stated for the 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 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 on improvement to an existing method.

2.4.2. State the Research Questions

After defining the study’s aim, relevant questions concerning 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 [3].

2.4.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 2022.12.0+353, using the rugarch [4,24], SimDesign [6], tidyverse [7], zoo [40], aTSA [41] and forecast [42] 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 situations 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 a different sequence of random numbers each time it is run unless a seed is set [43]. A set seed initialises the random number generator [4] and ensures reproducibility, where the same result is obtained for different runs of the simulation process [44]. The seed needs to be set only once, for each simulation, at the start of the simulation session [2,4], and it is better to use the same seed values throughout the process [2].
  Now, through the GARCH model, this study carries out an MCS experiment to ascertain whether the seed values’ pattern or arrangement affects the estimators’ efficiency and consistency properties. Two sets of seeds are used for the experiment, where each set contains three different patterns of seed values. The first set is $S_1$ = {12345, 54321, 15243}, while the second set $S_2$ = {34567, 76543, 36547}. In each set, the study tries to use seed values arranged in ascending order, then reverses the order, and finally mixes up the ordered arrangement. The simulation starts by using GARCH(1,1)-Student’s t, with a degree of freedom $\nu$ = 3, as the true model under four assumed error distributions of a Normal, Student’s t, Generalised Error Distribution (GED) and Generalised Hyperbolic (GHYP) distribution. Details on these selected error distributions can be seen in [24,45]. The true parameter values used are ($\mu$, $\omega$, $\alpha$, $\beta$) = (0.0678, 0.0867, 0.0931, 0.9059), and they are obtained by fitting GARCH(1,1)-Student’s t to the SA bond return data.
  Using each of the seed patterns in turn, simulated dataset of sample size N = 12000, repeated 1000 times are generated through the parameter values. However, because of the effect of initial values in the data generating process, which may lead to size distortion [46], the first N = {11000, 10000, 9000, 8000} sets of observations are each discarded at each stage of the generated 12000 observations to circumvent such distortion. That is, only the last N = {1000, 2000, 3000, 4000} are used under each of the four assumed error distributions, as shown in Table A1, Appendix A. These trimming steps are carried out following the simulation structure of Feng and Shi [28] 2. An observation-driven process like the GARCH can be size distorted with regards to its kurtosis, where strong size distortion may be a result of high kurtosis [47]. The extracts of the RMSE and SE outcomes for the GARCH volatility persistence estimator $\widehat{\alpha }+\widehat{\beta }$ are shown in Table A1. For $S_1$ in Panel A of the table, as N tends to its peak, the performance of the RMSE from the lowest to the highest under the four error distribution assumptions is Student’s t, GHYP, GED and Normal in that order, while that of SE from the lowest to the highest is GHYP, Student’s t, GED and Normal in that order, for the three arrangements of seed values.
  For $S_2$ in Panel B of the table, as N reaches its peak at 4000, the performance of the RMSE from the lowest to the highest is Student’s t, GHYP, GED and Normal in that order, while that of SE from the lowest to the highest is GHYP, GED, Student’s t and Normal in that order, for the three $S_2$ patterns of seed values. Hence, efficiency and precision in terms of RMSE and SE are the same as the sample size N becomes larger under the three seeds, regardless of the arrangement of the seed values under $S_2$, as also observed under $S_1$. In addition, the flows of $\sqrt{N}$ consistency of the estimator under the seed values in $S_1$ are roughly the same; this is also applicable to those of the seed values in $S_2$. The plotted outcomes can be visualised as displayed by the trend lines within the 95% confidence intervals in Figure 2 for the three seed values of sets $S_1$ in Panel A and $S_2$ in Panel B, where the efficiency and consistency outcomes are roughly the same with increase in N.
  To summarise, this study observes that, as $N\to \infty$, the pattern or arrangement of the seed values does not affect the estimator’s overall consistency and efficiency properties, but this may likely depend on the quality of the model used. The seed is primarily used to ensure reproducibility. Panels C and D of the figure further reveal that the RMSE/SE → 0 as $N\to \infty$ for the four error distributions in $S_1$ and $S_2$.
  Table A1 further shows that the MCS estimator $\widehat{\alpha }+\widehat{\beta }$ considerably recovers the true parameter $\alpha +\beta$ at the 95% nominal recovery level, where some of the estimates even recover the complete true value (0.9990) with TPR outcomes of 95%. These recovery outcomes can be seen in the visual plots of Figure 3 (or as shown in Panels A and B of Figure A1, Appendix B), where Panels A(i) and B(i) reveal that the MCS estimates perform quite well in recovering the true parameter as shown by the closeness of the TPR outcomes to the 95% (i.e., 0.95) nominal recovery level for $S_1$ and $S_2$, respectively. The bunched up TPR outcomes in Panels A(i) and B(i) are clearly spread out as shown in Panels A(ii) and B(ii) for $S_1$ and $S_2$, respectively. From these recovery outputs, two distinct features can be observed. First, the TPR results do not depend on the sample size as shown in Panels A and B of Figure 4 for $S_1$ and $S_2$, which is a feature of coverage probability (see [11]); second, the closer (farther) the MCS estimate is to zero, the smaller (larger) the TPR outcome, as revealed in Panels C and D of the figure.
• After setting the seed, the true parameter values of the true sampling distribution (or true model) are then set a priori [48,49].
• 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 GARCH model is carried out using the R package "rugarch". Random data generation involving this package can be implemented using either of two approaches. The first approach is to carry out the data generating simulation directly on a fitted object "fit" using the ugarchsim function [4,24] for the simulated random data. The second approach uses the ugarchpath function, which enables simulation of desired number of volatility paths through different parameter combinations [4,24].
  The simulation or data generating process can be run once or replicated multiple times. This study carries out another MCS investigation through the GARCH model to determine the effect (on the outcomes) of running a given GARCH simulation once or replicating it multiple times. That is, for a given sample size and seed value, the outcome of running the simulation once is compared to that of running it with different replications like 2500, 1000, and 300. This MCS experiment uses GARCH(1,1)-Student’s t, with $\nu$ = 3, as the true model under four assumed error distributions of a Normal, Student’s t, GED and GHYP. However, it should be understood that any non-normal error distributions (apart from the Student’s t that is used here) can also be used with GARCH(1,1) model for the true model. The GARCH(1,1)-Student’s t fitted to the SA bond return data yields the true parameter values ($\mu$, $\omega$, $\alpha$, $\beta$) = (0.0678, 0.0867, 0.0931, 0.9059).
  Using these parameter values, datasets of sample size N = 12000 are generated in each of the four distinct simulations (i.e., simulations with 1, 2500, 1000, and 300 replicates). After necessary trimmings in each simulation, to evade initial values effect, the last N = {1000, 2000, 3000} sets of observations are used at each stage of the generated 12000 observations under the four assumed innovation distributions. That is, datasets of the last three sample sizes, each simulated once, then replicated L = {2500, 1000, and 300} times are consecutively generated. From the modelling outputs, it is observed that the log-likelihood (llk), RMSE, SE and bias outcomes of $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\alpha }+\widehat{\beta }$ estimators for each simulation under the four assumed errors are the same for the three sample-size datasets with the same seed value, regardless of whether the simulation is run once or replicated multiple times. For brevity, this study only displays the outcomes of the experiment under the assumed GED error for each run in Table 1. However, increasing the number of replications may reduce sampling uncertainty in meta-statistics [6].
• The generated (simulated) data are analysed, and the estimates from them are evaluated using classic methods through meta-statistics to derive relevant information about the estimators. Meta-statistics (see [6]) are performance measures or metrics for assessing the modelling outputs by judging the closeness between an estimate and the true parameter. A few of the frequently used meta-statistical summaries, as described below, include bias, root mean square error (RMSE) and standard error (SE). For more meta-statistics, see [2,6,50].

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 [28]. The optimal value of bias is 0 [50,51]. 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 MCS (see [39]) as $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,\dots ,L)$ is a finite ith 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 [39,52]) in Equation (12). Also called Monte Carlo standard deviation, it is a measure of the efficiency or precision of the true parameter’s estimator, which 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 [2]). Sampling variability decreases with increased sample size [50].

$$SE=\sqrt{\frac{1}{L}\sum _{i=1}^L{({\widehat{\vartheta }}_i-\overline{\vartheta })}^2},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 [50] and it is stated as:

$$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 lesser RMSE is more efficient in recovering the true parameter value [50,52], and minimum RMSE produces maximum precision [53]. Consistency of the estimator occurs when RMSE decreases such that $\widehat{\vartheta }\to \vartheta$ as the sample size $N\to \infty$ [2,24]. RMSE is related to bias and sampling variability as:

$$RMSE=\sqrt{{bias}^2+{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 reliant on the law of large numbers, and it is expected that the distribution of an appropriately large sample should converge to that of the underlying population as the sample size increases [54]. It is also expected that the Monte Carlo sampling error should decrease as the sample size increases, but this is not always the case. That is, the sample size cannot always be sufficiently increased to limit the sampling error to a tolerable level [54].

2.4.4. Discussion and Summary

After implementing the method, the last stage in the framework steps is the conclusion, which needs to 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.

3. Results: Simulation and Empirical

3.1. Practical Illustrations of the Simulation Design: Application to Bond Return Data

By way of illustration, this section practically describes how the stated steps can be applied using Monte Carlo simulations with a real data empirical verification.

3.1.1. The Background

It is believed that observation-driven models can appropriately estimate volatility when fitted with a suitable error distribution [55]. 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 [56,57] for details). Data generation using the rugarch package can be done through a variety of models that include the simple GARCH, the exponential GARCH (EGARCH), the GJR-GARCH, the Component GARCH (CGARCH) [58], the Multiplicative Component GARCH (MCGARCH) [59], among others, and two omnibus models apARCH and fGARCH (as described in Section 2.2).

The apARCH model is less robust than the fGARCH model [24], hence the latter is used for the data generation in this study. Specifically, fGARCH(1,1) model is used as the true data generating process (DGP) for the MC simulation because the first lag of conditional variability can considerably capture volatility clustering existing in the time series data. In other words, the dependence of volatility on recent past activities is more than on distant past activities [60]. Hence, this illustrative study showcases the effectiveness of the observation-driven model fGARCH for estimating volatility, where the outcomes of the model fitted with each of ten assumed innovation distributions of the Normal, skew-Normal, Student’s t, skew-Student’s t, GED, skew-GED, GHYP, Normal Inverse Gaussian (NIG), Generalised Hyperbolic Skew-Student’s t (GHST) distribution and Johnson’s reparametrised SU (JSU) distribution are compared. Details on the error distributions can be found in [24,45,61,62,63,64,65,66,67].

The DGP fGARCH(1,1) model, as stated in Equation (15) is used to generate simulated return observations using the non-Normal Student’s t error with $\nu$ = 4.1 as the true error distribution.

$${\sigma }_t^{\gamma }=\omega +{\alpha }_1{\sigma }_{t-1}^{\gamma }\left(\right|z_{t-1}-{\lambda }_{21}{|-{\lambda }_{11}\{z_{t-1}-{\lambda }_{21}\})}^{\delta }+{\beta }_1{\sigma }_{t-1}^{\gamma }$$

(15)

Here, a Student’s t with shape parameter or degree of freedom $\nu$ = 4.1 is used to ensure that $\mathtt{E}\left[z_t^4\right]$ < ∞, which enables $\sqrt{N}$ consistency of the QML estimation following the assumption of Francq and Thieu [68] (see [69]). Moreover, the Student’s t distribution is used as the true error distribution in this study because it can suitably deal with leptokurtic or fat tailed features [70,71] experienced in financial data [29], and it is also assumed that stock prices appear to have a distribution much like the Student’s t [72]. However, based on relevance and research needs, users may choose to use any leptokurtic distributions, like the GED or others, for their data generation. Simulation through the rugarch package can be carried out using the ugarchsim and ugarchpath functions, but not all the stated data generating models currently support the use of ugarchpath method (see [24]). Hence, this illustration is implemented using the ugarchsim function through the "fit object" approach.

Further background study reveals the findings of Morris et al. [2], where the authors showed that RMSE is more applicable as a performance measure where the objective of the simulation is prediction rather than estimation. The authors discussed how more sensitive RMSE is to the choice of the number of observations used during method comparisons than when only SE or bias is used. Hence, for fairness in performance assessments, the SE is used as the key metric or measure of efficiency (precision) in this illustrative study.

It is also noticed from the outcomes of the family GARCH modelling that two sets of standard error (SE) estimates are returned. That is, the default MLE standard errors (SEs) and the robust QMLE SEs [24,73,74]. This study used the robust QMLE SEs of the fGARCH modelling for the simulation illustrations because they are claimed to be consistent (but not efficient) and asymptotically normally distributed if the volatility and mean equations are well specified [27,75].

3.1.2. Aim of the Simulation Study

This study aims at obtaining the most appropriate assumed error distribution for conditional volatility estimation when the underlying (true) error distribution is unknown.

3.1.3. Research Questions

This simulation study should result in responses to the following questions:

• Which among the assumed error distributions is the most appropriate from the fGARCH process simulation for volatility estimation?
• Financial data are fat-tailed [76], i.e., non-Normal. Hence, will the combined volatility estimator $\widehat{\alpha }+\widehat{\beta }$ of the most suitable error assumption still be consistent under departure from Normal assumption?
• What type (i.e., strong, weak or inconsistence) of $\sqrt{N}$ consistency, in terms of RMSE and SE, does the fGARCH estimator $\widehat{\alpha }+\widehat{\beta }$ exhibit?
• How is the performance of the MCS estimator $\widehat{\alpha }+\widehat{\beta }$ in recovering the true parameter?

3.1.4. Method of Implementation

To initiate the implementation method, the written code is first used to fit the true model fGARCH(1,1)-Student’s t to the SA bond return data through the ugarchfit function of the fGARCH fit object (see [4,24] for details on the fit object). Next, through the ugarchsim function, using seed 12345 in the code, the outputs from the fit are set (or used) a priori as the true parameter values ($\alpha$, $\beta$) = (0.0748, 0.9243) for the simulation process as shown in Table 2. These parameter values with other estimates from the fit object are utilised directly or unaltered to generate N = 15000 sample size observations, replicated 1000 times. However, after trimming down the simulated dataset, following the simulation structure of Feng and Shi [28], to prevent the effect of initial values, by removing the first N = {7000, 6000, 5000} sets of observations at each stage of the simulated 15000 observations, the last N = {8000, 9000, 10000} are processed under each of the ten assumed innovation distributions as shown in Table 2.

Figure 5 displays the visual outlooks of the simulated returns and volatilities for the first three series in the 1000 replicated series for sample size N = 8000. These sampled visuals show that each of the 1000 replicated series of the simulated (synthetic) data has a unique randomness and shape that make them different from every other series. Hence, the estimate from the family GARCH simulation is the average of all the estimates from the different replicated series.

After generating simulated observations, the fGARCH(1,1) model is fitted to each simulated dataset under the ten error assumptions. The full code that shows the command lines for the stages of the method of implementation is available from the authors on request. The parsimonious ARMA(1,1) model is also used in the code as the most suitable among the tested candidate ARMA models to remove serial correlation in the simulated observations. However, consistency can still be achieved in simulation modelling regardless of correlated sample draws. That is, the sampled variates do not need to be independent to achieve consistency [77].

Next, the selected meta-statistics are now used to evaluate the estimators. The most suitable assumed error distribution for volatility estimation will be obtained from the estimator/s with the best precision and efficiency from the meta-statistical comparisons done under all the selected error assumptions. Three meta-statistical summaries that include the bias, RMSE and SE, are used in this illustration. The computation of the metrics is direct but may sometimes be nerve-racking, and manual programming may even cause unanticipated coding errors and other abrupt setbacks. To circumvent this, SimDesign statistical package [6,50] with in-built meta-statistical functions for computational accuracy is used in this illustration, beginning from bias estimation. The log-likelihood (llk) of the estimates, with the RMSE, bias and SE for $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\alpha }+\widehat{\beta }$ estimators are displayed in Table 2, but SE is the key performance measure for efficiency and precision.

Now comparing RMSE for $\widehat{\alpha }$, the results from the table show that both the skew-GED and skew-Student’s t outperform the other assumed innovations in efficiency with the least values as N tends to the peak at 10000. For $\widehat{\beta }$, the JSU, followed by the GHST and NIG, outperforms the rest of the innovation assumptions in efficiency with the least RMSE value as N tends to the peak. For $\widehat{\alpha }+\widehat{\beta }$, the GHST followed by the skew-Student’s t outperform the remaining eight innovation assumptions as N tends to the peak, but the skew-Student’s t is the best as the sample size reaches the middle at N = 9000 for both $\widehat{\alpha }+\widehat{\beta }$ and $\widehat{\beta }$.

Comparing bias for $\widehat{\alpha }$, as N approaches the peak, the absolute values of biases for the GED and skew-Student’s t outperform the rest. For $\widehat{\beta }$, the JSU and the true innovation Student’s t both take the lead as N reaches the peak. For $\widehat{\alpha }+\widehat{\beta }$, the JSU followed by the GHYP outperform the other innovations in absolute values of biases as N reaches the peak.

For precision and efficiency comparison in terms of the key performance metric SE, the skew-Student’s t relatively outperforms the others in efficiency and precision as N tends to the peak for $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\alpha }+\widehat{\beta }$, in particular, for $\widehat{\beta }$ and $\widehat{\alpha }+\widehat{\beta }$. Finally, for the llk comparison, the GHYP outperforms the rest, with the largest estimates at the three N sample sizes. To summarise, when the true innovation is Student’s t, the skew-Student’s t assumed innovation distribution relatively outperforms the other nine innovation assumptions in efficiency and precision, while the GHYP performs better than the rest through the log-likelihood. It is observed here that the SEs of $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\alpha }+\widehat{\beta }$ estimators for the assumed Normal innovation distribution are the largest when compared with those of the other nine assumed innovation distributions. This justifies the claim that the QMLE of the family GARCH model (with Normal innovation) is inefficient. Furthermore, it is observed from the tabulated outputs that the RMSE and SE of the estimators are considerably $\sqrt{N}$ consistent in recovering the true parameters under the assumed innovations. The visual illustrations of the consistency for the outputs of $\widehat{\alpha }+\widehat{\beta }$ estimator are graphically displayed in Figure 6.

The figure shows that the closer the absolute values of biases are to zero, the closer the SE is to the RMSE. Whenever bias drifts away from zero, the gap between SE and RMSE widens, but if otherwise, then their trend lines closely follow the same trajectory. The visual plots also show that RMSE and SE decrease as N increases, but bias is independent of N.

It is also observed from the table and as shown in Panels A and B of Figure 7 that the MCS estimates for the estimator $\widehat{\alpha }+\widehat{\beta }$ considerably recover the true (volatility) parameter value of 0.9991, with TPR outcomes closely clustered around the 95% (i.e., 0.95) nominal recovery level under the ten error assumptions. This indicates a good performance of the MCS experiments with suitably valid outputs. However, the non-Normal errors perform slightly better in the recovery than the Normal and skew-Normal errors, as clearly revealed in Panel B. It can also be seen from the tabular results that the TPR outcomes are independent of N, and the closer the MCS estimate is to zero, the smaller the TPR estimate.

3.2. Empirical Verification

Next, the outcomes of the MCS experiments are now empirically verified using the real return data from the SA bond market index. Among the ten assumed error distributions, the most appropriate for the fGARCH process to estimate the volatility of the bond market’s returns is examined. For the market index, the price data are transformed to the log-daily returns by taking the difference of logarithms of the price, expressed in percentage as:

$$r_t=ln\left(\frac{P_t}{P_{t-1}}\right)\times 100.$$

(16)

The $P_t$ and $P_{t-1}$ are the closing bond price index at time t and the previous day’s closing price at time $t-1$, respectively; r${}_t$ is the current return, and ln represents the natural logarithm.

3.2.1. Exploratory Data Analysis

To start with, the price index and returns are first inspected through exploratory data analysis (EDA) as displayed in Figure 8. The EDA visually sheds light on the content of the dataset to reveal relevant information and potential outliers. Figure 8 unearths some downswings or steep falls in the volatility of price (in plot a) and returns (in plot c) around the years 2002, 2008, 2016 and 2020. The most recent as shown by the plots for 2020 was due to the global covid-19 pandemic.

For further inspection, the figure is now separated into two panels: left and right. The left panels contain plots a, b, e, f for daily bond prices, while the right panels consist of plots c, d, g, h for the returns. The left panels reveal non-stationarity in the price index as observed in the price series plot, the density plot, quantile-quantile (QQ) plot and the box plot. On the other hand, the right panels show stationarity in the returns through the return series plot, the density plot, the QQ plot and the box plot. These summarily elucidate the non-stationarity in the SA daily bond prices and stationarity in the returns.

3.2.2. Tests for Serial Correlation and Heteroscedasticity

Next, linear dependence (or serial correlation) and heteroscedasticity are filtered out by fitting ARMA-fGARCH models with each of the ten innovation distributions to the stationary return series. ARMA(1,1) model, as stated in Equation (17), is found to be the most adequate, among all the examined candidate ARMA(p,q) models, to remove serial correlation from the SA bond market’s return residuals. Table 3 presents the outcomes of the Weighted Ljung-Box (WLB) test (see [78] for details) for fitting the ARMA(1,1) model. The p-values of the test at lag 5 all exceed 0.05 under each error distribution. Based on this, we fail to reject the null hypothesis of "no serial correlation" in the SA bond market’s returns. This means there is no evidence of autocorrelation in the return residuals.

$$r_t={\varphi }_0+{\varphi }_1{r}_{t-1}+{\phi }_1{\epsilon }_{t-1}+{\epsilon }_t$$

(17)

Following the filtering of linear dependence in the return series, Engle’s ARCH test (see [14]) is carried out using the Lagrange Multiplier (LM) and Portmanteau-Q (PQ) tests to check for the presence of heteroscedasticity or ARCH effects in the residuals. These tests are implemented based on the null hypothesis of homoscedasticity in the residuals of an Autoregressive Integrated Moving Average (ARIMA) model. Both tests’ outcomes show highly significant p-values of 0 as shown in Figure 9. Hence, the null hypothesis of "no ARCH effect" in the residuals is rejected, which denotes the existence of volatility clustering. Based on this, a heteroscedastic model can be fitted to remove the ARCH effects in the series. To do this, the candidate robust fGARCH$(u,v)$ models, with each of the ten error distributions, are fitted to the SA bond returns, where the fit of the parsimonious fGARCH(1,1) model as shown in Equation (18) is found to be the most suitable.

$${\sigma }_t^{\gamma }=\omega +{\alpha }_1{\sigma }_{t-1}^{\gamma }\left(\right|z_{t-1}-{\lambda }_{21}{|-{\lambda }_{11}\{z_{t-1}-{\lambda }_{21}\})}^{\delta }+{\beta }_1{\sigma }_{t-1}^{\gamma }.$$

(18)

After fitting the fGARCH model to the returns, the weighted ARCH LM test is used to ascertain if ARCH effects have been filtered out. The p-value of the "ARCH LM statistic (7)" at lag 7 in Table 3 is greater than 5% under each of the ten innovation distributions. Hence, this indicates that heteroscedasticity is filtered out since we fail to reject the null hypothesis of "no ARCH effect" in the residuals. These outcomes show that the variance equation is well specified.

3.2.3. Selection of the Most Suitable Error Distribution

Next, selection of the most suitable assumed error distribution to describe the market’s returns, when fitted with the fGARCH model for volatility estimation is obtained from Table 3. It is observed from the table that all, but two, of the fGARCH volatility parameter estimates ($\widehat{\omega }$, $\widehat{\alpha }$, $\widehat{\beta }$, ${\widehat{\lambda }}_{11}$, ${\widehat{\lambda }}_{21}$ and $\widehat{\gamma }$) under the ten innovation assumptions are statistically significant at 1% level. This means that these parameters are actively needed in the model. The two exceptions are the insignificant $\widehat{\omega }$ for the Normal and skew-Normal, and the estimates ${\widehat{\lambda }}_{21}$ that are mostly not significant or barely significant. The strongly significant ${\widehat{\lambda }}_{11}$ indicates the dominance of asymmetric large shocks in the return series.

Comparison of the error distributions are carried out using the log-likelihood and four information criteria that include the Akaike information criterion (AIC), Bayesian information criterion (BIC), Hannan-Quinn information criterion (HQIC) and Shibata information criterion (SIC) (see [24] for details). The largest log-likelihood value with the smallest values of the information criteria under a given assumed innovation indicates that it is the most appropriate innovation distribution to describe the market for volatility estimation.

It is observed from Table 3 that the values of all four information criteria are smallest under the skew-Student’s t innovation distribution, but the GHYP innovation has the highest log-likelihood value. Hence, the skew-Student’s t is the most suitable innovation assumption strictly based on the information criteria, while the GHYP is the most appropriate if the decision is made using the log-likelihood. The GHYP and skew-Student’s t also yield better goodness of fit (GoF) outcomes when compared with the remaining eight errors, as shown by their large p-values in the table, which shows that they are the best fit among the ten error assumptions for the distribution of the SA bond’s return residuals. Hence, the volatility of the SA bond market’s returns can be most suitably estimated through ARMA(1,1)-fGARCH(1,1) model fitted with the GHYP or skew-Student’s t assumed error distribution. These empirical results are consistent with the Monte Carlo simulation modelling outcomes. This study also checked the empirical outcomes of fitting the less omnibus apGARCH(1,1) model to the bond return data and we arrived at the same results (of skew-Student’s t through information criteria and GHYP via log-likelihood) obtained by fGARCH(1,1) model (see Table 4). The table only shows the outcomes of the log-likelihood and information criteria for brevity.

For the run-time, it is observed that the skew-Student’s t is about four (approximately 4.2) times faster than the GHYP for both simulation and empirical modelling. That is, it takes the GHYP about four times the computational time it takes the skew-Student’s t to run the same process. Since the empirical and simulation run-times are approximately the same, we only present the empirical run-times for the ten innovations in Table 3 to conserve space. For both simulation and empirical runs, the GHYP has the highest runtime among the ten innovations, followed by the NIG, while the Normal has the least.

From the outputs of the ARMA(1,1)-fGARCH(1,1) model in Table 3, the mean and variance (from the conditional standard deviation’s Box-Cox transformation in Section 2.2) equations of the model fitted with each of the GHYP and skew-Student’s t are stated as:

$$\begin{array}{ccc}\hfill \mathrm{With}\mathrm{GHYP}:r_t& =& {\mu }_t+{\epsilon }_t\hfill \\ & =& 0.0156+{\epsilon }_t\hfill \\ \hfill {\sigma }_t^{\gamma }& =& \omega +{\alpha }_1{\sigma }_{t-1}^{\gamma }\left(\right|z_{t-1}-{\lambda }_{21}{|-{\lambda }_{11}\{z_{t-1}-{\lambda }_{21}\})}^{\delta }+{\beta }_1{\sigma }_{t-1}^{\gamma }\hfill \\ \hfill {\sigma }_t^{1.2086}& =& 0.0267+0.0661{\sigma }_{t-1}^{1.2086}\left(\right|z_{t-1}-0.0942{|-0.3370\{z_{t-1}-0.0942\})}^{1.2086}+0.9241{\sigma }_{t-1}^{1.2086}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \mathrm{With}\mathrm{skew-Studen}{\mathrm{t}}^{’}\mathrm{s}t:r_t& =& {\mu }_t+{\epsilon }_t\hfill \\ & =& 0.0177+{\epsilon }_t\hfill \\ \hfill {\sigma }_t^{\gamma }& =& \omega +{\alpha }_1{\sigma }_{t-1}^{\gamma }\left(\right|z_{t-1}-{\lambda }_{21}{|-{\lambda }_{11}\{z_{t-1}-{\lambda }_{21}\})}^{\delta }+{\beta }_1{\sigma }_{t-1}^{\gamma }\hfill \\ \hfill {\sigma }_t^{1.2058}& =& 0.0270+0.0661{\sigma }_{t-1}^{1.2058}\left(\right|z_{t-1}-0.0943{|-0.3445\{z_{t-1}-0.0943\})}^{1.2058}+0.9236{\sigma }_{t-1}^{1.2058}\hfill \end{array}$$

4. Discussion and Summarised Conclusion

In conclusion, it is observed that using the uGARCHsim approach for the MCS illustration, the GHYP and skew-Student’s t evolve as the most suitable assumed error distributions to reckon with and use with the fGARCH model for volatility estimation of the SA bond returns when the underlying error distribution is unknown. These outcomes are verified empirically. The persistence of volatility under these most suitable error distributions are 0.9795 for the GHYP and 0.9792 for skew-Student’s t. Hence, this indicates that the volatility of the SA bond market’s returns is considerably (highly) persistent.

The conclusion under this section continues by providing answers to the four research questions. In this study, consistency is termed "strong" when the estimator’s RMSE/SE value decreases as N increases without distortion. Otherwise, it is weak. Now, answering the questions: first, the GHYP and the skew-Student’s t distributions are the most appropriate among the stated assumed error distributions from the fGARCH process simulation for the volatility estimation. Second, the volatility estimator $\widehat{\alpha }+\widehat{\beta }$ for each of the most suitable assumed error distributions GHYP and skew-Student’s t is strongly $\sqrt{N}$ consistent for both RMSE and SE under departures from the Normal assumption as revealed in Panels G and D of Table 2. Third, there are strong $\sqrt{N}$ consistency for the RMSE and SE of the fGARCH estimator $\widehat{\alpha }+\widehat{\beta }$ under all, but one, of the ten assumed error distributions as shown in Table 2. The lone exception however is the weak consistency in the SE of the Normal assumption. Fourth, as a proxy for the coverage of the MCS experiment, the MCS estimator $\widehat{\alpha }+\widehat{\beta }$ performed well in recovering the true parameter $\alpha +\beta$ through the TPR measure for the 95% nominal recovery level as revealed in Table 2 and Figure 7. The results show that the TPR outcomes are suitably close to the 95% nominal recovery level under the ten error assumptions.

5. Conclusions

This study showcases a robust step-by-step framework for a comprehensive simulation by presenting the functionalities of the rugarch package in R for simulating and estimating time-varying parameters through the family GARCH observation-driven model. The framework hands out an organised approach to Monte Carlo simulation (MCS) study that involves "background (optional), defining the aim, research questions, method of implementation, and summarised conclusion". The method of implementation is a workflow that includes writing the code, setting the seed, setting the true parameters a priori, data generation process, and performance assessment through meta-statistics.

This novel, easy-to-understand framework is illustrated using financial return data; hence, users can easily use it for effective MCS studies. With the uGARCHsim simulation approach involved in the modelling, the implementation method is clearly explained with relevant details. Key observations are identified, and novel findings brought to light. The framework also outlays clear guidelines for data generation using the package, since data generation is without a doubt an integral part of MCS studies. The key observations and novel findings in this study include first, it is shown in the experiment that as the sample size N becomes larger, the consistency and efficiency properties of an estimator in a Monte Carlo process are generally not affected by the pattern or arrangement of the seed values, but this may depend on the quality of the model used. Hence, regardless of the arrangement of the seed values, the efficiency and consistency of an estimator generally remain the same as N tends to infinity.

Second, it is investigated and revealed in this study that the outcomes of the GARCH MCS experiments are the same regardless of whether the simulation or data generating process is run once or replicated multiple times. Third, this study derived a "true parameter recovery (TPR)" measure as a proxy for the coverage of MCS experiment. This new (original) novel measure is flexible to apply and can henceforth be used by upcoming researchers to determine the level of recovery of the true parameter value by the MCS estimates. It is also observed that the volatility estimator of the used fGARCH model displays considerably strong $\sqrt{N}$ consistency. Lastly, the outcomes of the illustrative study show that the GHYP and skew-Student’s t errors are the most suitable among the ten assumed innovations to describe the SA bond returns for volatility estimation. The fit of these two error assumptions with the fGARCH model revealed considerably high volatility persistence in the returns. On a wider scale, since volatility is a practical measure of risk, the fit of the GHYP and skew-Student’s t errors with a specification of the fGARCH (or, apARCH) model for a robust volatility estimation may benefit financial institutions and markets by enhancing the accuracy of their risk estimations. This could potentially lead to a significant reduction in asset losses. It is anticipated that researchers will leverage this study’s novel findings and robust design for improved simulation studies. It is also anticipated that the outcomes of this study will be broadly applicable in finance and other sectors.

5.1. Limitations in the Study

Three limitations or challenges are noticed in this study. The first is on how to obtain a sufficient sample size that can generate accurate outcomes. To tackle this using the illustrative example, the process involves testing a selected number of sample sizes, with each used in turn, until a pattern of efficiency and/or consistency starts to evolve under the stated error distributions. The error distribution with the most efficient in terms of the given performance measure under a particular sample size is carefully noticed. If a set of sample sizes yields the same efficiency outcomes, the outcome with the best consistency among the set can be used for a final decision on sample size determination. This is a guide to obtaining the required sample size/s.

The second is running time. It is observed that a large simulated dataset may sometimes be needed to obtain accurate computations and this may be done at the cost of a large computational (or running) time, depending on the model used. This can be very demanding, especially when dealing with different stages of large sample sizes. Third, since the rugarch package does not make provision for calculating the coverage probability, this study derived a proxy for the coverage using the TPR measure, and it is observed that the MCS estimates considerably recover the true parameters.

5.2. Future Research Interest

The authors intend to further use the ugarchpath function of the rugarch package through any of the models that support its use for the framework illustration. The authors also intend extending the simulation framework ideas to other volatility estimating models, like the Generalised Autoregressive Score (GAS) model, and to modelling and estimating multivariate processes. The future extension also includes a framework for volatility forecasting and portfolio management.