Multi-objective real-time integrated solar-wind-thermal power dispatch by using meta-heuristic technique

1.   Introduction

The risk-return relationship topic has been studied from as early as the 1950s, with growing enthusiasm to date (Ma et al. 2020). The foremost model used to capture the risk premium is the Generalized Autoregressive Conditional Heteroscedasticity-in-Mean (GARCH-M) by Engle et al. (1987). The GARCH approach, originated by Engle (1982), assumes that price data and innovations follow a normal distribution. However, such an assumption is based on "simplicity and convenience" because in the real world, the nature of the financial market is subject to nonlinear properties, such as asymmetric volatility and asymmetric effects (Kuang, 2020; Ayed et al. 2020). As a result, the standard GARCH model evolved to asymmetric GARCH models such as EGARCH, APARCH and GJR, to capture asymmetric volatility and the asymmetric effects - volatility feedback and the leverage effect. While the GARCH approach has the ability to capture the volatile nature of financial data, the model is unable to effectively capture asymmetry, as sufficed by Mangani (2008), Mandimika and Chinzara (2012), Ilupeju (2016), Naradh et al. (2021) and Dwarika et al. (2021).

Another extension that has gained popularity over the years, to assist with capturing asymmetry, is the combination of GARCH models with various innovation distributions (Delis et al. 2021; Kuang, 2020). According to Alqaralleh et al. (2020), innovation distributions that capture the empirical nature of returns, such as higher moment properties, skewness and excess kurtosis, make "excellent tools" in the estimation of accurate forecasts. While there is existing literature on different probability distributions, only one study by Delis et al. (2021), to the best of the authors knowledge, applied a different distribution to investigate the risk-return relationship. Drawn from Ma et al. (2020) and Tsay (2013), a reason for the lack of literature is the limited ability to manipulate the modelling of the time-varying risk-return relationship. According to Tsay (2013), when modelling the GARCH approach in R, if the model is invalid by failing to converge or is insignificant, the mean and/ or constant can be dropped to make the model valid. Since the GARCH-M, as the name suggests, contains the risk premium within the mean, it cannot be dropped (Ma et al. 2020; Tsay, 2013). Hence, this presents a gap in literature for the exploration of different innovation distributions with the GARCH-M approach to investigate the risk-return relationship.

By targeting this specific issue, this can improve the overall efficiency of the risk estimation ability of the GARCH approach. The problem statement is that due to the GARCH modelling assumption of normality, the risk-return relationship is underestimated, leading to inaccurate forecasts. This is especially relevant in the case of an emerging market, such as South Africa, which is subject to heavy tails due to high and persistent levels of market volatility (Dwarika et al. 2021; Naradh et al. 2021). At the same time, such a market presents "huge investment opportunities" for investors that seek a superior return (Ayed et al. 2020). In other words, the higher the risk taken, the higher the expected rate of return. The relationship between risk and return provides valuable information to structure investment strategies, assist budget plans and allocate resources. Given the wide use of the GARCH approach in risk management, as documented by Saddah and Sitanggang (2020) and Ayed et al. (2020), acknowledging and implementing improvements in stochastic volatility models can improve forecast accuracy and lead to more financially sound decisions.

There is no existing study, to the best of the authors knowledge, that has conducted a comparative analysis of different innovation distributions in the context of the GARCH-M approach. Therefore, the aim of this study is to investigate the optimal GARCH-M model and innovation distribution, to determine the risk-return relationship. Hence, the primary aim is to determine the optimal GARCH-M type model to investigate the risk-return relationship. The secondary aim is to find an innovation distribution to optimize the GARCH-M tails, i.e., optimally capture risk. To ensure optimal model diagnostics, this study follows the international standard procedure, established by the regulatory framework - the Basel Accords. The Basel Accord I recommends Value at Risk (VaR) and backtesting methods to ensure the GARCH model's efficiency (Ayed et al. 2020; Kuang, 2020). VaR computes the maximum loss of a portfolio over a specific period of time for a given confidence interval (Kuang, 2020). The backtesting procedure is used to test whether the actual losses are in line with the forecasted values (Ayed et al. 2020). Essentially, it is used as a tool to determine the optimal forecast power of the VaR estimates.

This study consists of five sections. First, the Introduction provides an overview of the background and what this study aims to achieve. Second, the Literature Review outlines existing South African studies on the GARCH-M approach with different innovation distributions and then extends to international work. This section is concluded by identifying gaps found in the body of literature. Third, the Methodology covers the data, GARCH type models and innovation distributions employed in this study. Fourth, the Empirical Results and Discussion consist of the implications of the results alongside the model output. Finally, the key findings of this study are summed up in the Conclusions.

2.   Literature review

According to Markowitz (1952), risk and return hold a positive and linear relationship in theory. Therefore, when an investor makes an investment in an emerging market that is characterized by high-risk, the investor expects a higher rate of potential return, in compensation for investing in a volatile environment that does not guarantee a superior return. However, early studies by Mangani (2008), Mandimika and Chinzara (2012) and Adu et al. (2015) found no risk-return relationship in the largest South African stock market at the time - the Johannesburg Stock Exchange (JSE). All three studies employed the GARCH approach, but found limitations in the model's ability to efficiently capture risk, especially by the model's innovations.

Dwarika et al. (2021) consolidated several limitations of the GARCH approach. The study analyzed the daily JSE All Share Index (ALSI) returns for the sample period from 15 October 2009 to 15 October 2019. The models employed were GARCH (1, 1), EGARCH, GJR and APARCH, with the standard normal (NORM), Student-t (Stud-t), Skewed Student-t (Skew-t) and the Generalized Error distribution (GED). The study found EGARCH to be the optimal model, but the innovations showed that asymmetry remained uncaptured, in line with Mangani (2008), Mandimika and Chinzara (2012) and Ilupeju (2016). While the overall finding was a positive risk-return relationship, contrasting results were found for the Skew-t distribution, as GARCH-M found no relationship, whereas EGARCH-M found a positive relationship. Emenogu et al. (2020) also found that the Skew-t performance varied with the different GARCH models. Thus, this inconsistency highlighted the significance of the innovation distribution choice, as it has the ability to affect parameter estimation.

Ayed et al. (2020) investigated the optimal GARCH model and innovation distribution by analyzing daily data, over the sample period from 10 August 2006 to 14 December 2014, of the Islamic stock market. The study employed the standard GARCH and APARCH with the distributions - the NORM, Stud-t and Skew-t. The study found APARCH and the Skew-t to be optimal. From the VaR and backtesting procedure, at a 99% level of VaR, the Stud-t and Skew-t were optimal, whereas for a 95% VaR, the forecast of NORM was more robust. However, the GARCH (1, 1) is critiqued by Saddah and Sitanggang (2020) for underestimating risk. This is because a simple GARCH model is unable to effectively capture higher moment properties (Mandimika and Chinzara, 2012; Kuang, 2020). Similarly, applying GARCH type models with the conventional and normal type distributions, to a volatile emerging market, is considered ineffective (Ayed et al. 2020; Bautista and Mora, 2020; Saddah and Sitanggang, 2020). Thus, other asymmetric properties need to be taken into account.

For instance, Delis et al. (2021) found a positive risk-return relationship, in line with theoretical expectations, confirming the significance of accounting for skewness to assist with capturing asymmetry. The authors modelled the daily returns of the S & P500 index, over the sample period from 02 January 1980 to 14 October 2020, using GARCH-M with a Skewed GED. The study concluded that skewness is vital in the estimation of the risk-return relationship, especially during periods of excess volatility such as the COVID-19 pandemic. However, Khan et al. (2021) investigated the recent effects of the pandemic on the Indian stock market using a different approach, namely, the Extreme Value Theory (EVT). This approach specifically focuses on tail behavior, where extreme events such as the 2008 financial crisis and COVID-19 occur. Khan et al. (2021) used the Generalized Pareto distribution (GPD) to model the high-frequency 5-min data of the NIFTY 50 returns for the sample period 11 March 2020 to 30 September 2020. During the pandemic, it was found that returns dropped by 9%, and that 5% of the data was above the specified threshold, indicating an asymmetric return distribution, in line with empirical expectations. Khan et al. (2021) concluded the efficiency of the GPD in the field of EVT to provide useful information during the period of a crisis.

In contrast, Kuang (2020) found that the Pearson Type Ⅳ distribution (PIVD) outperformed the Stud-t, Generalized Extreme Value distribution (GEVD) and GPD using VaR. The GEVD is another innovation distribution from the EVT. The PIVD is flexible due to its ability to account for higher moment properties, particularly skewness and excess kurtosis. The author analyzed the daily data of the indices S & P500, FTSE100, CAC40, and DAX30 from 03 May 2002 to 31 December 2009. The study concluded the superiority of the PIVD in capturing skewness, which affects the accuracy of risk estimation, in line with Delis et al. (2021). Although not investigated in the study by Kuang (2020), the author mentioned another robust distribution, the Stable distribution, which accounts for heavy tails, skewness and excess kurtosis. Bautista and Mora (2020) applied VaR to the oil sector using GARCH with the distributions - NORM, Generalized Stud-t and Stable. The authors found that the NORM underperformed because the oil returns followed a heavy-tailed distribution. The Generalized Stud-t was more optimal than the NORM, but overall Stable was the most robust distribution, in line with Ilupeju (2016), Bautista and Mata (2020) and Naradh et al. (2021).

In conclusion, it can be noted that there are two major research gaps in literature. Firstly, the South African studies are limited to the conventional innovation distributions: NORM, Stud-t, Skew-t and GED, in the context of the risk-return relationship, as sufficed by Mangani (2008), Mandimika and Chinzara (2012) and Dwarika et al. (2021). While the Stud-t, Skew-t and GED can account for some degree of asymmetry, it is vital to consider more flexible distributions that can take into account heavy tails and other higher moment properties such as skewness and excess kurtosis. In turn, this will enhance the risk efficiency of the GARCH approach and improve the forecast accuracy of the risk-return relationship. Secondly, only one existing study by Delis et al. (2021), investigated the GARCH-M with a different innovation distribution - the Skewed GED. A possible reason for the lack of literature is due to limited model manipulation arising from software constraints. To overcome this limitation, it would be advisable to investigate several GARCH type models. In line with international literature, this study extends to the VaR and backtesting procedure since none of the South African studies employed this international regulatory framework, in the context of the GARCH-M approach.

3.   Methodology

Figure 1 shows a flow chart diagram of the methodology employed in this study to determine the optimal GARCH-M and innovation distribution.

The JSE ALSI market returns are modelled by the GARCH-M (1, 1) type models - standard GARCH, GJR, EGARCH and APARCH with the innovation distributions: NORM, Stud-t, Skew-t and GED. The optimal GARCH-M model is selected by information criteria. The innovations are then extracted to determine if risk remains uncaptured by the "risk check". If risk remains uncaptured, the more robust innovation distributions are employed: PIVD, GEVD, GPD and Stable. Model diagnostics are used to determine if the distributions are adequate in capturing asymmetry. Thereafter, the VaR and backtesting methods are employed to determine the optimal innovation distribution.

3.1.   Data

Following Dwarika et al. (2021) and Naradh et al. (2021), this study will obtain the daily price data of the JSE ALSI from the Integrated Real-time Equity System (IRESS) database, previously known as the McGregor Bureau for Financial Analysis (BFA). The IRESS database is a reliable source of obtaining data and is hosted on the IRESS Research Domain, where market data is available for financial analysis and academic research purposes (IRESS, 2022).

Following Emenogu et al. (2020), this study investigates a duration of 17 years made up of a total of 4251 observations. The updated sample period analyzed is from 05 October 2004 to 05 October 2021. This period included the 2008 financial crises and COVID-19 pandemic. These are global economic events, but more importantly, extreme events that affect tail behavior which this study takes into account, in line with Kuang (2020) and Khan et al. (2021).

Following Dwarika et al. (2021), the risk-free rate proxy is the South African T-bill, which shall be obtained from the South African Reserve Bank (SARB). The ALSI price data is converted to returns by ${R}_{m} = \mathrm{l}\mathrm{n}\left(\frac{P{\text{­}}_{t}}{{P}_{t-1}}\right)$, where ${R}_{m}$ = market returns and ${P}_{t}$ = share price for the current day $t$ and ${P}_{t- 1}$ = share price for the previous day $t-1$. The annual risk-free rate is converted to a daily value by $\mathrm{D}\mathrm{a}\mathrm{i}\mathrm{l}\mathrm{y}\ {R}_{f} = {\left(1+\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{l}\mathrm{y}\ {R}_{f}\right)}^{\left(\frac{1}{365}\right)}-1$, given by Brooks (2014). Excess returns are computed as the difference between market returns and the risk-free rate (Delis et al. 2021).

Following convention in literature, the data is tested for stationarity to ensure a valid time series, by employing the Augmented Dickey-Fuller (ADF), Phillips-Perron (PP) and Kwiatkowski, Phillips, Schmidt and Shin (KPSS) test. The data is also tested for normality and heteroscedasticity to show the asymmetric nature of returns. Thus, substantiating the employment of the GARCH extensions which can account for such higher moment properties and consequentially improve forecast accuracy.

3.2.   GARCH-M Approach

Following Dwarika et al. (2021), Naradh et al. (2021), Ilupeju (2016) and Mandimika and Chinzara (2012), the foremost GARCH-M type models are estimated by the Maximum Likelihood (ML) method, in combination with the four conventional innovation distributions: NORM, Stud-t, Skew-t and GED. The GARCH-M models follow the mean model of the standard GARCH by Engle (1982):

where ${y}_{t}$ = conditional mean, $\mu$ = constant, ${\delta }_{j}$ = risk premium, ${{\sigma }^{2}}_{t-j}$ = variance and ${u}_{t}$ = innovation term.

The conditional variance of the asymmetric GARCH-M models is listed respectively. GJR by Glosten et al. (1993):

where ${\alpha }_{0}$ = constant term, ${{u}^{2}}_{t-1}$ = squared lagged innovation term, ${\alpha }_{1}$ = ARCH effect (short-term volatility persistence), ${\beta }_{j}$ = GARCH effect (long-term volatility persistence), $\gamma$ = asymmetry parameter (to capture asymmetric volatility and asymmetric effects) and ${I}_{t-1}$ = interaction variable (accounts for positive and negative shocks).

EGARCH by Nelson (1991):

where ${\alpha }_{1}[\left|{z}_{t-1}\right|-E|{z}_{t}\left|\right]$ = magnitude of the innovation and $\gamma {z}_{t-1}$ = sign effect (accounts for positive and negative shocks).

APARCH by Ding et al. (1993):

where $\delta$ = is an exponent that allows the APARCH model to take on several GARCH models. For example, in this study, the GJR-GARCH is implemented by setting δ of the APARCH model to 2, in line with Tsay (2013).

The optimal GARCH-M is selected by model diagnostics, information criteria, Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC), following standard literature. The standardized innovations of the optimal GARCH-M and innovation distribution are then extracted.

Following Dwarika et al. (2021), which built on Mangani (2008), Mandimika and Chinzara (2012) and Ilupeju (2016), the standardized innovations of the optimal GARCH-M are tested for normality, heteroscedasticity and randomness by a framework of preliminary tests. This study will formalize such testing as a "risk check" to determine the extent of uncaptured risk, more specifically, due to which property of risk is being underestimated (the asymmetric, volatile or random nature of the innovations). Table 1 shows the "risk check" framework.

Based on the aforementioned South African studies, the GARCH-M model's innovations are expected to show uncaptured risk, to substantiate the application of the more robust probability distributions: PIVD, GEVD, GPD and Stable. Following Bautista and Mora (2020), Bautista and Mata (2020), Naradh et al. (2021) and Ilupeju (2016), the innovation distributions are estimated by the ML method.

3.3.   Innovation Distributions

3.3.1.   PIVD

According to Kuang (2020), the PIVD is described as flexible due to its ability to account for higher moment properties, particularly skewness. The PIVD is given as:

where $\lambda$ = location, $\theta$ = scale, $k$ = normalization constant, $m$ = kurtosis and $\nu$ = skewness.

The location parameter $\lambda$ is the mean and the scale $\theta$ is the standard deviation. The scale determines the width, whereas the location represents the shift of the mode of the distribution. The sign indicates the direction in which the distribution shifted. For example, if positive, then the distribution is shifted to the right. The purpose of the normalization constant $k$ is to ensure the function is a probability distribution function.

In the context of this study, focus is placed on the higher moment properties of the distributions. The kurtosis parameter $m$ controls the tail thickness, where low values indicate thin tails and large values indicate heavy tails. Specifically, values more than 2.5 indicate heavy tails (Ilupeju, 2016). The direction of skewness is indicated by the sign of the parameter $\nu$, where $\nu > 0$ is positively skewed and $\nu < 0$ is negatively skewed (Kuang, 2020).

3.3.2.   GEVD

The aim of the EVT is to examine the behavior of outliers (extreme values) of a random variable. The advantage is that focus is placed on tail behavior instead of the entire distribution (Echaust and Just, 2020). Such an approach is relevant, given the recent COVID-19 pandemic, as sufficed by Omari et al. (2020), Khan et al. (2021) and Delis et al. (2021).

Following Echaust and Just (2020), Kuang (2020) and Khan et al. (2021), the EVT is characterized by two modelling approaches, the Block Maxima Model (BMM) and Peaks Over Threshold (POT). The POT approach uses a predetermined threshold to determine extreme values, whereas the BMM chooses a maximum value for a given block which is defined as a specified period of time. The GEVD falls under the BMM, whereby the innovations are divided into equal lengths, representing consecutive blocks. The maximum innovation values are then selected from each block. The choice of block size is vital as a small value can indicate estimation bias, whereas a large value can result in estimation variance. A block value of five or ten is sufficient to ensure accurate results (Ilupeju, 2016). Echaust and Just (2020) give the GEVD as:

where $\lambda$ = location, $\theta$ = scale and $\xi$ = shape.

The shape $\xi$, also known as the tail index, describes the nature of the distribution. If $\xi < 0$, this indicates a Weibull distribution, a short-tailed distribution defined as having a finite right tail. If $\xi > 0$, the distribution is a Frechet type (inverse Weibull distribution) that is heavy-tailed, where examples include Stud-t and Stable. If $\xi = 0$, this indicates a thin-tailed distribution such as NORM and log-NORM.

Following Khan et al. (2021), model diagnostics in the form of normality tests are employed to determine whether the GEVD is an adequate fit for the innovations. This study employs the Anderson Darling (AD) test, probability plot (PP), quantile-quantile (QQ) plot, return level plot and density plot. The AD normality test focuses on tail behavior. If the AD test statistic is greater than the critical value, the null hypothesis that the innovations are normally distributed can be rejected. With respect to the plots, a major deviation between the empirical and theoretical distribution, would indicate nonnormality, thus an inadequate fit of the GEVD.

3.3.3.   GPD

According to Echaust and Just (2020), the POT approach is considered more advantageous than the BMM, because it focuses on all possible events greater than the threshold and not just the largest events. The aim of POT is to assess the data and model all the data points that are above the specified threshold $u$. Making an accurate selection of the threshold is vital for analysis.

The conventional approach is to use graphical analysis - mean residual life plot and parameter stability plot - to determine the threshold value. The mean residual life plot shows the mean excesses against the threshold of the standardized innovations for both tails. The selected threshold would occur where there is an upward linear trend in the slope of the mean excess values (Khan et al. 2021). A parameter stability plot is employed to show the optimal shape and scale parameters against the threshold values that range from 1.0 to 2.0 (Naradh et al. 2021). Such methods consist of identifying stable regions in the graphs which are considered highly subjective, by Echaust and Just (2020). In order to confirm the threshold, a Pareto quantile plot is employed to state the respective number of positive and negative innovations above the threshold. The innovations are fitted to the GPD and then the shape and scale parameters are estimated by the ML method. Echaust and Just (2020) give the GPD as:

where $\theta > 0$, $x\ge 0$ for $\xi \ge 0$ and $0\le x\le -\frac{x}{\xi }$ for $\xi < 0$. The GPD has just two parameters, $\theta$ = scale and $\xi$ = shape.

A thin-tailed distribution is indicated by $\xi < 0$, whereas a heavy-tailed distribution such as Stable, is indicated by $\xi > 0$. Similar to the GEVD procedure, the model diagnostics employed are a PP, QQ, return level and density plot.

3.3.4.   Stable

The school of $\alpha$-Stable or Stable distributions is a widely documented class of probability distributions that can effectively model excess kurtosis, skewness, volatility clustering, asymmetry and heavy tails (Bautista and Mata, 2020; Kuang, 2020). Bautista and Mora (2020) give the Stable distribution as:

where $-\infty \le \lambda \le \infty \in R$ = location, $\theta > 0$ = scale, $-1\le v\le 1$ = skewness, and $0 < \alpha \le 2$ = tail exponent.

The sign of $v$ indicates the direction of skewness and is subject to the given constraint above. The tail exponent, also known as the index of stability, represents the rate at which the tails diminish or come to an end. When $\alpha > 2$, the mean is equivalent to the location $\lambda$, but no mean exists when $\alpha < 1$. If $\alpha < 2$, this results in an infinite variance and the tails are thus heavy.

Following Bautista and Mata (2020) and Naradh et al. (2021), the model diagnostics include the Kolmogorov-Smirnov (KS) and AD normality tests. KS computes the difference between the theoretical and empirical distribution, whereas AD computes the mean of square differences. A variance stabilized PP plot and density plot are employed to show a graphical analysis of the distribution fit.

The four fitted innovation distributions are analyzed by model diagnostics, thereafter the optimal distribution is determined by VaR and backtesting methods.

3.4.   VaR and backtesting

3.4.1.   VaR

VaR is a well-known measure of market risk and is used in the fields of portfolio management and risk management (Bautista and Mata, 2020; Emenogu et al. 2020; Omari et al. 2020). Following Kuang (2020), the estimation of the VaR is made up of the GARCH-M in combination with the innovation distributions. VaR is given as $P\left({RR}_{t}\le VaR{\text{­}}_{t}\left(\alpha \right)\right) = \alpha$, which is defined as the probability of the risk-return relationship $RR$ at time $t$ less than or equal to the VaR estimate at the percentage of time $\alpha$. Delis et al. (2021) stated that there are both upper and lower tails of a distribution that correspond to $\alpha$. In the context of this study, the upper tails are known as the short position, associated with gains, where the confidence intervals are defined as 97.5 and 99%. Similarly, the lower tails are known as the long position, associated with losses and are defined as 1 and 2.5%. The highest VaR estimates indicate high-risk potential losses and biased values, where the distribution is inadequate; the lowest VaR estimates indicate the opposite (Ayed et al. 2020; Bautista and Mata, 2020).

3.4.2.   Violation ratio

Following several studies, such as Kuang (2020), Ayed et al. (2020) and Echaust and Just (2020), the backtesting procedure is used to analyze the forecast power of the VaR estimates. In this study, the backtesting procedure consists of the violation ratio (VR), the Kupiec unconditional coverage test, Christoffersen conditional coverage test and VaR duration-based test. The calculation of the violation ratio (VR) is given by Emenogu et al. (2020) as:

where if VR < 1, risk is underestimated and if VR > 1, risk is overestimated. If 0.8 ≤ VR ≤ 1.2, the forecasted risk is optimal, whereas 0.5 < VR < 1.5 is inadequate.

The VR serves as a preliminary test as it is a very simple method to check model adequacy; therefore, more powerful tests are employed. According to Kuang (2020) and Ayed et al. (2020), an optimal VaR model satisfies the unconditional and conditional coverage properties. The unconditional property states that the observed number of violations should be in line with expected values at a given confidence interval. If the observed and expected values are inconsistent, this phenomenon would be referred to as the failure rate which is defined as the number of violations. The conditional property states that such violations are independent of each other. Thus, the more robust tests of the backtesting evaluation criteria consist of the Kupiec unconditional coverage test and Christoffersen conditional coverage test.

3.4.3.   Unconditional test

To investigate the unconditional property, the proportion of failure test by Kupiec (1995) is employed. The null hypothesis is that the observed number of violations is consistent with the expected values at a given confidence interval. The Kupiec test statistic is given as:

where p = the assumed probability of occurrence, m = number of hits of the model (number of times the loss is greater than the VaR estimate) and n = number of tests.

The test statistic ${K}_{UC}$ follows a chi-squared distribution with one degree of freedom $\left({X}^{2}\left(1\right)\right)$. If ${K}_{UC}$ is less than the critical value, do not reject the null hypothesis and conclude model adequacy. If ${K}_{UC}$ is greater than the critical value, reject the null and conclude model inadequacy.

3.4.4.   Conditional test

To investigate the conditional property, the conditional joint test by Christoffersen (1998) is employed, where the test statistic is given as:

The test statistic consists of the Kupiec unconditional test and Christoffersen (1998) independence test. For the latter, the null hypothesis is that the failure rate is independent, indicating model adequacy. The test statistic is given as:

where ${\widehat{\pi }}_{ij}$ = number of times that state $i$ follows state $j$, state 0 = independence and state 1 = dependence.

The decision procedure for the test statistic ${C}_{IND}$ follows the Kupiec unconditional test. It should be noted that for the joint decision procedure, the test statistic ${L}_{CC}$ follows a chi-squared distribution with two degrees of freedom $\left({X}^{2}\left(2\right)\right)$. A limitation of the independence test is that it only tests for one-day violations instead of over multiple days (Echaust and Just, 2020). Thus, a VaR duration-based test is employed to investigate the time dynamics between the past and future violations.

3.4.5.   Duration-based test

For the VaR duration-based test by Christoffersen and Pelletier (2004), the expectation is that given any day, future violations should not occur because of the past. The null hypothesis is that the probability of a violation on any day is subject to the no memory property, whereby it is independent of the effects of the previous day(s) and has a mean distribution of $1/p$.

To explain the latter, consider the Weibull parameter = $b$, where the Weibull distribution is of a school of continuous probability distributions that can take on several shapes. When $b = 1$, it follows an exponential distribution which is the only distribution that is memory free. Thus, if the time between violations has no memory of each other, the average is $1/p$ as defined by the exponential distribution. The Weibull distribution is a function of the failure rate to an index of time, given as:

where $d$ = number of days, $b$ = shape parameter and $a$ = scale parameter.

The test statistic is defined as:

where $L$ = likelihood, $\widehat{a}$ = estimated scale parameter and $p$ = probability of no memory violation.

The decision procedure for the test statistic $D$ follows the Kupiec unconditional test.

Following standard literature, this study uses the estimated p-values to make the decisions for the hypothesis tests. The highest p-value indicates the optimal model (Ayed et al. 2020; Omari et al. 2020).

4.   Empirical results and discussion

4.1.   Data exploration

Following Dwarika et al. (2021) and Naradh et al. (2021), the daily closing price data of the JSE ALSI are obtained from IRESS for the sample period from 04 October 2004 to 05 October 2021. The ALSI price data is converted to market returns by taking the natural log form of the difference between closing prices of the current and previous day. Following Dwarika et al. (2021), the annual T-bill obtained from the SARB is converted to a daily value to compute excess returns. The ALSI (excess) return data is made up of a total of 4251 observations. The ALSI returns is confirmed to be stationary, by the ADF, PP and KPSS tests, to avoid a spurious regression and ensure a valid time series. The statistical properties of the ALSI returns are investigated by standard preliminary tests. Table 2 shows the preliminary test results of the ALSI returns.

From Table 2, since the p-values are less than 5%, the null hypothesis that the ALSI returns follow a normal distribution is rejected at a 5% level of significance. The null hypothesis that the ARCH effect is absent within the ALSI return series is also rejected at a 5% level of significance. In conclusion, the series of the ALSI returns are asymmetric, nonnormal and volatile in nature; thus, substantiating the employment of the GARCH-M type models.

4.2.   GARCH-M approach

Following Dwarika et al. (2021), Naradh et al. (2021), Ilupeju (2016) and Mandimika and Chinzara (2012), the foremost GARCH-M type models and innovation distributions are investigated. The GARCH-M models are shown for the four innovation distributions: NORM, Stud-t, Skew-t and the GED. Table 3 shows the relevant significant ML parameter estimates of the GARCH-M (1, 1) models with the different innovation distributions.

From Table 3, the GARCH-M models are valid since the ARCH and GARCH effects are significant at all the levels of significance (1, 5 and 10%). The sum of the ARCH and GARCH effects ${(\widehat{\alpha }}_{1}+{\widehat{\beta }}_{1})$ is positive and high, indicating strong and persistent levels of volatility in the South African market. The only exception is for APARCH-M Skew-t, with a sum greater than one, indicative of an overestimation of risk, in line with Dwarika et al. (2021). All the models have an asymmetry parameter $\widehat{\gamma } > 0$, which indicates the presence of asymmetric volatility and asymmetric effects. The presence of high levels of volatility, asymmetric volatility and asymmetric effects are in line with the empirical expectations of an emerging market and prior South African studies.

The risk premium $\widehat{\delta }$ is statistically significant at a 5 and 10% level of significance for the GARCH-M approach. The majority of the GARCH-M models have a positive $\widehat{\delta }$, indicating a positive risk-return relationship. This finding is in line with Dwarika et al. (2021) but contrasts to earlier studies, by Mangani (2008), Mandimika and Chinzara (2012) and Adu et al. (2015), who found no relationship. The significance of the risk premium of the asymmetric GARCH-M models improved in comparison to the standard GARCH-M. This is due to the more robust fit of the asymmetric GARCH-M models to the asymmetric nature of the financial data, as suggested by Saddah and Sitanggang (2020).

In terms of theoretical expectations, it follows that a positive risk-return relationship will exist. That is, since the South African market is characterized as high-risk, an investor can expect a superior rate of return. This was confirmed by the positive risk premium found by the majority of the GARCH-M test results. The only exceptions are GJR-M Stud-t and APARCH-M Skew-t, where the risk premium is negative. The change in sign is due to modelling adjustments made, such as the dropping of a constant, to ensure the model is valid. The optimal GARCH-M is investigated by model diagnostics. Table 4 shows the information criteria of the GARCH-M models with the different innovation distributions.

From Table 4, the consistent optimal fitting innovation distribution to the respective GARCH-M models is Skew-t, in line with Dwarika et al. (2021), Ayed et al. (2020) and Omari et al. (2020). This finding is in contrast to Emenogu et al. (2020), who found that the Skew-t performance varied with different GARCH type models. The NORM is the least optimal fit for EGARCH and APARCH, as expected by Mandimika and Chinzara (2012), Kuang (2020) and Saddah and Sitanggang (2020). The asymmetric financial data would make a more robust fit with distributions designed to account for heavy tails and other higher moment properties.

EGARCH is found to be the optimal model in this study, in line with Dwarika et al. (2021), Naradh et al. (2021), Saddah and Sitanggang (2020) and Omari et al. (2020). This is in contrast to APARCH which was found to be optimal by Ilupeju (2016) and Ayed et al. (2020). Overall, Skew-t is the optimal fitting innovation distribution in combination with EGARCH-M, in line with Dwarika et al. (2021) and Omari et al. (2020). The risk check is employed to investigate the properties of uncaptured risk.

Following Dwarika et al. (2021), which built on Mangani (2008), Mandimika and Chinzara (2012) and Ilupeju (2016), the standardized innovations of the EGARCH-M Skew-t model are investigated by the risk check to determine the extent of uncaptured risk. Table 5 shows the risk check test results for the innovations of the EGARCH-M Skew-t model.

From Table 5, the normality tests show that the innovations are asymmetric in nature, indicating that asymmetry has been inadequately captured. Since the p-values of the heteroscedasticity tests are significant, this means that volatility is present within the innovations and remains uncaptured. This finding is in contrast to Dwarika et al. (2021) and Ilupeju (2016). It could be due to this study's selected sample period, which included highly volatile periods, such as the 2008 financial crisis and recent COVID-19 pandemic.

The randomness tests reveal that the innovations exhibit random behavior, hence, an indication of uncaptured risk. In conclusion, the asymmetric, volatile and random nature of the innovations is inadequately captured by EGARCH-M Skew-t. The overall finding of uncaptured risk is in line with Mangani (2008), Mandimika and Chinzara (2012), Ilupeju (2016) and Dwarika et al. (2021). The failed risk check motivates the employment of more robust nonnormal innovation distributions

4.3.   Different innovation distributions and model diagnostics

4.3.1.   PIVD

The different innovation distributions are fitted to the standardized innovations extracted from EGARCH-M Skew-t. Table 6 shows the ML parameter estimates of the PIVD.

From Table 6, $\widehat{m}$ is 8.9423 which indicates a positive and high kurtosis value. More specifically, a heavy-tailed distribution, as expected, in line with a PIVD and the market characteristics of an emerging market. The sign of the skewness parameter 6.3842 indicates a positively skewed distribution. The AD test indicates that the innovations follow a normal distribution. Thus, it can be concluded that the PIVD outperforms the Skew-t distribution and is a more robust fit for the innovations, in line with Ilupeju (2016) and Kuang (2020).

4.3.2.   GEVD

Following Ilupeju (2016), the block size used in this study is five to ensure accurate results. Table 7 shows the ML parameter estimates of the GEVD with the Standard Errors (SE) in brackets. Let ${z}_{t}$ = positive innovations and ${z}_{t}^{\mathrm{*}} = -{z}_{t}$ = negative innovations.

From Table 7, the tail index represented by the shape parameter $\widehat{\xi }$ is less than zero for both the positive and negative innovations. This indicates a Weibull distribution, a short-tailed distribution defined as having a finite right tail. The AD test confirms that the GEVD is an adequate fit, since the innovations follow a normal distribution, in line with the PIVD above. Figures 2 and 3 show the model diagnostics of the GEVD for the standardized innovations.

From Figures 2 and 3, the PP plot shows a perfect fit of the innovations to the GEVD. However, from the QQ plots, the data points do show some departure from the theoretical straight line. From the return level plots, most of the points that deviate from the straight line still fall on or within the confidence bands; therefore, the fit is adequate. The density plot confirms the adequate fit as most of the points lie on the GEVD. According to Echaust and Just (2020), the GPD is considered more advantageous than the GEVD, because it focuses on all possible events greater than the threshold and not just the largest events.

4.3.3.   GPD

The threshold selection is made, by graphical methods, before the GPD analysis. Figures 4 and 5 show the mean innovation life plots for the standardized innovations.

From Figures 4 and 5, the selected threshold is estimated to be approximately 2. The parameter stability plots are further employed to explore the selected threshold values. Figures 6 and 7 show the parameter stability plots for the standardized innovations.

From Figures 6 and 7, the stable region is from the threshold values 1.1 to 1.4, respectively. The Pareto quantile plots are employed to confirm the threshold values. Figures 8 and 9 show the Pareto quantile plot for the standardized innovations.

From Figures 8 and 9, the threshold for ${z}_{t}$ is $u$ = exp (0.1782694) = 1.1951 and ${z}_{t}^{*}$, $u$ = exp (0.3915099) = 1.4792. The selected threshold values are used to estimate the ML parameter estimates of the GPD. Table 8 shows the ML parameter estimates of the GPD.

From Table 8, the number of observations above the threshold values is 438 and 326, respectively. These values are higher than those found in the studies by Naradh et al. (2021), Omari et al. (2020) and Ilupeju (2016), due to a larger sample period analyzed in this study. The scale parameters are significant and positive, in line with Naradh et al. (2021), Omari et al. (2020) and Ilupeju (2016). The shape parameters are negative, $\widehat{\xi } < 0$, indicating a thin-tailed Weibull distribution. This finding contrasts with the empirical expectations of an emerging market and the result of Naradh et al. (2021). However, such an exception was also found in the studies by Omari et al. (2020) and Ilupeju (2016). Figures 10 and 11 show the model diagnostics of the GPD for the standardized innovations.

From Figures 10 and 11, the PP plot shows that the empirical and theoretical distributions are in alignment. However, similar to the case of the GEVD above, the data points from the QQ plots show deviation from the theoretical straight line, especially for the negative innovations. From the return level plots, the data points either lie on or fall within the confidence bands, indicating an adequate fit. The density plot confirms the adequate fit as most points lie along the GPD. In conclusion, the model diagnostics of the EVT distributions, the GEVD and GPD, are in line with Naradh et al. (2021), Khan et al. (2021) and Ilupeju (2016).

4.3.4.   Stable

Table 9 shows the ML parameter estimates of the Stable distribution and the normality test statistics.

From Table 9, skewness represented by $\widehat{v}$ is negative, indicating that the innovations are negatively skewed. Since $\widehat{\alpha }$ < 2, this results in an infinite variance and heavy tails, in line with empirical expectations. Both normality tests show that the innovations follow a normal, or rather, a Stable distribution. Hence, Stable is a robust fitting probability distribution to the ALSI returns. Figures 12 and 13 show the graphical model diagnostics of the Stable distribution.

From Figure 12, there appears to be an almost perfect fit between the empirical and theoretical fit of the innovations to the Stable distribution. Figure 13 confirms an adequate fitting distribution by the empirical and theoretical lines being in perfect alignment. The optimal fit of Stable is in line with Naradh et al. (2021), Bautista and Mora (2020), Bautista and Mata (2020) and Ilupeju (2016).

4.4.   VaR and backtesting

Following several studies, such as Kuang (2020), Ayed et al. (2020) and Echaust and Just (2020), the forecast power of the VaR estimates of EGARCH-M with the different innovation distributions are analyzed by the backtesting procedure, which consists of the VR, Kupiec unconditional test, Christoffersen conditional test and VaR duration test. Table 10 shows the relevant VaR estimates for different levels (1, 2.5, 97.5 and 99%).

From Table 10, the highest VaR estimates are produced by the GPD and GEVD at the short and long positions, respectively. In contrast, the robust VaR estimates are produced by Stable at 1%, PIVD at 2.5% and the GPD at the long position. Overall, GEVD has the highest VaR estimate and Stable the lowest. The EVT distributions are considered inadequate in fitting the innovations for the short position. Khan et al. (2021) examined only the long position and found the highest VaR estimate for the GPD at 95% in the Indian stock market. With respect to South Africa, this finding is in line with Naradh et al. (2021), who found that the GPD had the highest VaR estimate, in the short position and long position, for the JSE Mining Index and ALSI, respectively.

This finding supports Ilupeju (2016), who found the highest VaR estimate for the GPD at 1%, but contrasts to GEVD which had the lowest VaR estimates at 2.5%. In the study by Bautista and Mata (2020), the authors found that Stable had the lowest VaR estimates at the examined long position for the Mexican Stock Exchange. Naradh et al. (2021) also found Stable to be optimal at the short and long positions. This contrasts with Ilupeju (2016), who found Stable with the highest VaR with the exception of the 1% level. The VaR estimates are examined by the backtesting procedure, as recommended by the Basel Accord, to determine the optimal innovation distribution. Table 11 shows the VR estimates for the different levels.

From Table 11, the majority of the VRs indicate an optimal forecast since they are between 0.8 and 1.2. The GEVD and GPD in the short position are the exceptions, as the VRs are well above 1, indicating that risk has been overestimated. This is in line with the robustness of the EVT distributions as indicated by the VaR estimates above. Tables 12 and 13 show the p-values of the Kupiec unconditional and Christoffersen conditional joint tests, respectively. Table 14 shows the p-values of the VaR duration-based test.

From Tables 12 and 13, it can be seen that the values of the conditional joint test and unconditional test are present for the same levels. Both the PIVD and Stable passed the conditional and unconditional coverage tests at the given levels of 1, 2.5, 97.5 and 99%. The optimal performance of Stable is in line with Naradh et al. (2021), Bautista and Mata (2020) and Ilupeju (2016). However, since the p-values for the PIVD are the highest relative to Stable, PIVD is the optimal innovation distribution. Model inadequacy is indicated for the remaining levels, as well as both the EVT distributions for the overall short and long position, due to an over or underestimation of the realized values. The null hypothesis is rejected for the long position of the GEVD and 99% level of GPD, consolidating model inadequacy of the EVT distributions.

From Table 14, the null hypothesis that the time between violations has no memory is rejected for all the innovation distributions of the long position. The same decision is made for both the GEVD and GPD in the short position, meaning that the EVT innovation distributions are incorrectly specified. The null is not rejected for the levels 1 and 2.5%, indicating model adequacy, where the majority of the p-values for the PIVD are higher than Stable.

In conclusion of the backtesting procedure, the PIVD is the optimal innovation distribution due to higher p-values relative to Stable and the EVT distributions which are found to be inadequate. The outperformance of the PIVD and Stable compared to the EVT distributions are in line with Naradh et al. (2021), Kuang (2020) and Ilupeju (2016). This is in contrast to Khan et al. (2021), who concluded the GPD is efficient and Omari et al. (2020), who concluded that the GARCH approach in the context of the EVT distributions is optimal. The sample periods used in the study by Khan et al. (2021) and Omari et al. (2020), were less than a year and fourteen and a half years, respectively. Therefore, the EVT distributions are considered more optimal in studies with smaller sample periods. In addition, it is vital to consider several innovation distributions and the backtesting procedure to determine which distribution optimizes the tails of any GARCH model.

5.   Conclusions

The only known study, to the best of the authors knowledge, to investigate a different innovation distribution (Skewed GED) in the context of GARCH-M, was by Delis et al. (2021). Therefore, this was the first study, on a local and international level, to investigate different nonnormal innovation distributions in combination with GARCH-M type models. The aim was to investigate the optimal GARCH-M model and innovation distribution to determine the risk-return relationship. The primary aim was to determine the optimal GARCH-M type model. This study found EGARCH-M as the optimal model, in line with Naradh et al. (2021), Saddah and Sitanggang (2020) and Omari et al. (2020). More specifically, EGARCH-M Skew-t was found to be optimal. However, the innovations failed to capture asymmetry, in line with previous South African studies by Mangani (2008), Mandimika and Chinzara (2012), Ilupeju (2016) and Dwarika et al. (2021). Therefore, this study employed more robust innovation distributions, which led to the secondary aim, which was to find an innovation distribution to optimize the GARCH-M model tails, i.e., optimally capture risk.

The nonnormal innovation distributions employed were the PIVD, GEVD, GPD and Stable. These distributions were considered more robust and flexible because they have the ability to take into account heavy tails and other higher moment properties such as skewness and excess kurtosis. This is relevant in the case of emerging markets, such as South Africa, which are characterized by heavy tails due to high levels of volatility. Model diagnostics revealed that the employment of the different innovation distributions was robust as there were no major deviations between the empirical and theoretical distributions. As the Basel Accord I recommended, more rigorous testing, VaR and backtesting methods were employed. From the VaR and backtesting analysis, it was found that the EVT distributions, the GEVD and GPD, were found to be the least robust relative to Stable and the PIVD. While the results of Stable and the PIVD were close, the PIVD was the optimal innovation distribution as it had the higher p-values.

Therefore, investors, risk managers and policymakers would opt to use the EGARCH-M in combination with the PIVD when modelling the risk-return relationship in the South African market. The robustness of EGARCH-M is supported by previous studies by Dwarika et al. (2021), Naradh et al. (2021), Saddah and Sitanggang (2020) and Omari et al. (2020). Stable was found to be an optimal innovation distribution by Bautista and Mora (2020), Bautista and Mata (2020) Naradh et al. (2021) and Ilupeju (2016). However, in the context of the risk-return relationship, this study found PIVD to be optimal, followed by Stable. This is in line with the study conducted by Kuang (2020), who found that the PIVD outperformed several other distributions, namely, the Stud-t, GEVD and GPD.

The EVT distributions, the GEVD and GPD, were found to be suboptimal in this study, in line with Naradh et al. (2021), Kuang (2020) and Ilupeju (2016). Based on prior studies, Khan et al. (2021) and Omari et al. (2020), the EVT distributions were found to be optimal given that smaller sample periods were used. Thus, it would be optimal to consider several distributions to ensure an unbiased estimation of the most robust innovation distribution, irrespective of the sample period used in a study.

To conclude, this study confirmed that the standard GARCH-M and conventional normal type innovation distributions are ineffective in fitting the asymmetric, volatile and random nature of financial data. It should be noted that this study focused on asymmetry and employed normality tests for the model diagnostics of the innovation distributions. Therefore, the first recommendation is to include heteroscedasticity and randomness tests to investigate other properties, such as the volatile and random nature of the innovations. In other words, to apply the complete risk check to the fitted innovation distributions. The second recommendation is to investigate the risk-return relationship using GARCH-M type models in direct combination with the more robust nonnormal innovation distributions. Finally, to explore other nonnormal innovation distributions in the context of GARCH-M, that were not included in the scope of this study, such as the Inverse Gaussian, Generalized Hyperbolic, Johnson's SU and Skewed GED.

Acknowledgments

The author would like to thank Prof. John Nolan for the grant received for the STABLE package in R (http://www.robustanalysis.com/).

Conflict of interest

The author declares no conflict of interest.

Data

The secondary price dataset used in this study can be obtained from the IRESS database.