Modelling Crude Oil Returns Using the NRIG Distribution

Over the past decade, crude oil prices have risen dramatically, making the oil market very volatile and risky; hence, implementing an efficient risk management tool against market risk is crucial. Value-at-risk (VaR) has become the most common tool in this context to quantify market risk. The main aim of the study is to propose models for accurate risk assessment in the crude oil and gasoline market using subclasses of generalised hyperbolic distributions (GHDs). We introduce the new subclass of GHDs, namely normal reciprocal inverse Gaussian distribution (NRIG), in evaluating the VaR for the crude oil and gasoline market. Furthermore, VaR estimation and backtesting procedures using the Kupiec likelihood ratio test are conducted to test the extreme tails of these models. The main findings from the Kupiec likelihood test statistics suggest that the best GHD model should be chosen at various VaR levels. Thus, the final results of this research allow risk managers, financial analysts, and energy market academics to be flexible in choosing a robust risk quantification model for crude oil and gasoline returns at their specific VaR levels of interest. Particularly for NRIG, the results suggest that a better VaR estimation is provided at the long positions.


Introduction
Crude oil is among the most important energy commodities in the world and forms the basis for many products, including transport fuels such as gasoline, diesel and jet fuel. According to the EIA 2017 report, 71% of oil is used for transport, 24% is based on industrial energy, and 3%, 2% and 1% are used for residential, commercial and electrical energy, respectively. The report also states that gasoline is the most widely used petroleum product in the United States. Finished gasoline consumption averaged about 9.33 million b/d (392 million gallons per day) in 2017, which is equivalent to about 47% of the U.S.s total oil consumption. Under normal circumstances, demand and supply are among the main drivers of crude oil prices. In addition, political events, extreme weather, and speculation in the financial markets are, among others, main factors of the crude oil market, increasing the degree of price volatility in the oil markets. Sales and revenues of major sectors worldwide and capital budgeting plans, as well as financial instability in oil-exporting and oil-consuming countries, are significantly affected by price volatility. The volatility of the energy commodity price does not only have a significant macroeconomic impact, but it also has an impact on stock markets Sadorsky (1999Sadorsky ( , 2003. Therefore, quantifying risk for oil-related commodities is very important to economic agents and policymakers, although it still remains one of the most challenging issues because of many fluctuations over time. Many risk measurement tools have been proposed in existing literature to give financial institutions, risk managers and market investors appropriate and simple technical approaches to measure the risk of financial and energy markets. One of the most common risk 206 MODELLING CRUDE OIL RETURNS USING THE NRIG DISTRIBUTION 2. Literature review

Modelling and analysing risk in the energy market
The benefits of VaR have long been discovered by many academic researchers in energy economics, and a wide range of research has been published, particularly in the last 15 years, during which the interest in energy commodities has increased dramatically due to the high financial importance of these products, along with extreme uncertainty in the energy market. Researchers have been comparing the various methodologies and strategies to decide on the best possible way to estimate VaR for a particular energy commodity or portfolio for a particular data sample. Cabedo and Moya (2003) introduced a version of the traditional HS, the historical simulation ARMA forecast (HSAF). They compared the VaR estimates produced by HS and HSAF with those produced by the basic GARCH model, assuming a normal distribution of returns, using an 8-year sample (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) of daily Brent oil prices. The researchers found that the HSAF is the most suitable methodology among these three VaR estimations, since it gives more robust VaR measurement than HS that better fits the extremely high price volatility, whereas the standard normal GARCH overestimates VaR.
Furthermore, the results of Cabedo and Moya (2003) were verified by Sadeghi and Shavvalpour (2006), and it was shown that HSAF produced the most precise VaR estimates in a 6-year sample (1997)(1998)(1999)(2000)(2001)(2002)(2003) of the weekly OPEC oil price data series compared to the basic ARCH and GARCH models. The authors observed that both ARCH and GARCH had still produced overestimated VaR, although they argue that whatever approach is used, VaR is an effective and robust instrument for measuring oil price risks.
Another improvement of the standard HS approach, the exponentially decreased frequency with the ARMA forecast approach (EDFAAF), was presented by Sadeghi and Shavvalpour (2006), based on the HSAF theory. Their approach was evaluated by comparing it with the HSAF approach, using a 12-year data sample of weekly Brent oil spot price observations, and it was concluded that the VaR forecasting was always better than the HSAF approach. Costello et al. (2008) argue that the results of Cabedo and Moya (2003) and Sadeghi and Shavvalpour (2006) are largely based on the assumption that the distribution of returns in the ARCH and GARCH models examined is computed by a normal distribution, an assumption that is empirically denied by many researchers for both the financial and, in particular, the even more volatile energy market commodities. In order to tackle this issue, the researchers integrated the idea of Barone-Adesi, Giannopoulos and Vosper (1999) and changed the standard technique for estimating the VaR using the GARCH model.
Especially Barone-Adesi et al. (1999) attempted to boost the GARCH models inadequate VaR forecasts, assuming that the normal distribution of returns is indicated by the use of historical simulations to determine the future distribution of returns. Their findings showed that the VaR estimates significantly outperform those of the standard GARCH model with their proposed GARCH-HS hybrid strategy. Costello et al. (2008) followed the key principle of the above mentioned semi-parametric technique and concluded that the GARCH semi-parametric method performs better than the HSAF of Cabedo and Moya (2003), because it captures the changing volatility of the daily Brent oil prices in a 23-year sample (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003).
The absence of capacity of the GARCH models to produce adequate VaR estimates when the innovations from the GARCH models are assumed to be normally distributed has given rise to the concern of many researchers, both in changing the standard normal GARCH model in such a manner that it would outperform other competing VaR estimation methodologies and in determining which alternative of the basic model would outperform the other GARCH models. Fan et al. (2008) suggested that GARCH models based on the generalised error distribution (GED) should be able to better capture characteristics such as fat tails and skewness, which are widely acknowledged to be found in highly volatile energy commodity returns. The researchers concluded that both GED-GARCH and GED-TGARCH models are more suitable for modelling the volatility characteristics of the daily spot price for WTI oil and Brent oil in a 20-year sample (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) as both of them outperformed the HSAF model in estimating the daily VaR at a 95% level of confidence.
Other research mainly focused on determining which particular innovation distribution will modify the standard ARCH and GARCH models and build a better VaR forecasting model. The skewed t-APARCH, skewed t-ARCH and risk metrics models were evaluated by Giot and Laurent (2003) on the basis of their forecasted VaR values. Using a 15-years sample (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) of daily spot prices for rare precious metals as well as WTI oil and Brent oil, the researchers concluded that the skewed t-APARCH model is the most suitable model for the fat-tailed and skewed distribution of WTI oil and Brent oil, and produces the best 1% VaR estimate.
Likewise, the basic GARCH model with the normal distribution, the Student GARCH model and the GARCH-HT model with the heavy-tailed distribution suggested by Politis (2004), were contrasted by Hung, Lee and Liu (2008). Using a 10-year sample (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) of daily spot prices for WTI oil, Brent oil, heating oil, propane and gasoline, the researchers noted that the GARCH-HT model is a robust model for VaR measurement, providing very precise estimates at high levels of confidence and thus making it appropriate for a more pragmatic risk management analysis.
The extreme value theory (EVT) method is another most common methodology for predicting VaR. It is either used independently or, in most cases, combined with the GARCH-type models, with the aim of removing the limitations of the standard methodology and achieving a more precise VaR forecast. In particular, as far as energy commodity risk management is concerned, a sufficient amount of research has been conducted with researchers integrating EVT; the great majority of them have tried to compare EVT models with other models that depend on other competitive techniques. For example, Krehbiel and Adkins (2005) examined the returns for WTI oil, Brent oil, heating oil, gasoline and natural gas for both spot and futures prices, and found evidence of heavy tails and negative skewness. Following the same framework as Byström (2005), they further concluded that the GARCH-EVT model (based on the POT approach) is superior to the GARCH family models, employing symmetrical distribution when modelling the risk of extreme price volatility and measuring the overall risk as in the VaR situation while also producing more precise results than the GARCH-EVT model using the BM approach.
Marimoutou, Raggad and Trableski (2009) modelled a 20-year sample (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) of daily Brent oil spot returns and made effective future VaR forecasts using the AR-GARCH model. Researchers altered methods such as HS and conditional EVT, and integrated the AR-GARCH model to produce models that were then evaluated based on their capacity to estimate VaR. The findings indicate that the AR-GARCH-EVT model, based on the POT approach and the AR-GARCH FHS model, outperform the AR-GARCH standard model and the traditional HS model. However, the AR-GARCH model with Student's t-distribution has been shown to be equally capable of producing VaR estimates, capturing and adapting to variations in volatility dynamics. Table 1 summarises some of the previous studies since 2010.

Generalised hyperbolic distribution and its subclasses in VaR estimation
The GHD of Barndorff-Nielsen (1997) has been widely used for modelling data with deviations from normality; not only does it capture skewness in the data, but it also accounts for extreme events that are not gradual in nature. It can therefore be considered particularly appropriate for modelling the energy returns under focus. Since its development, the GHDs have been successfully implemented in different disciplines such as physics and biology. Thereafter, several literature studies have been fitting the GHD subclasses to the financial data. Eberlein and Keller 1995 are among the first researchers who applied these distributions to finance; they fitted the hyperbolic subclass to the German financial data. Hence, the GHD and its subclasses generally fit the data better than the normal distribution.
Recently, studies have focused on the identification of an appropriate distribution and an adequate VaR model in the GHD class. For instance, Hu and Kercheval 2007 examined the generalised hyperbolic distributions (GHDs) for modelling equity returns compared to normal distributions. They described these GHDs and some of their properties, and tested them using a 6-years sample of daily S&P500 index prices. Their finding on the value at risk from calibrated distributions shows that the normal distribution leads to VaR estimates that significantly underestimate the realised empirical values, whereas the GHDs do not. Of the several GHD families considered (NIG, skewed t, VG and HYP), the most successful one was the skewed t-distribution. Huang et al. (2014) proposed three subclasses of the GHDs as appropriate models for JSE Mining Index returns. These models outperform the traditional assumption of normality and accommodate a number of stylised features,  such as skewness and volatility clustering, embedded within the financial data. The various criteria utilised suggested the generalised hyperbolic skewed Students t-distribution (GHST) as the most robust model for South African Mining Index returns. Chinhamu et al. (2015) evaluated the performance of the GHDs (the full GHD and its sub-classes, i.e., HYP, NIG, VG, and GHST) and the stable distribution in characterising the gold price returns. Their findings confirm that these distributions outperform the classical normality assumption of financial returns, and the full GHD and NIG were the most adequate models for the gold returns. Kemda et al. (2015) provided assessment of the adequacy of GHD for modelling the USD/ZAR exchange rate. The goodness-of fit-tests and the VaR estimates selected the VG as the optimal model over the NIG, GHST and HYP distributions. Mukhodobwane et al. (2020) evaluated the relative performance of ARMA(1, 1)-GARCH(1, 1) with normal, skewed-normal, Students t, skewed-Students t, generalized error distribution (GED), skewed-GED and the generalized hyperbolic (GHYP) distribution as innovations for BRICS stock market. They concluded that the Indian NIFTY markets volatility is best characterised by the generalized hyperbolic (GHYP) distribution.
Recently, the new subclass of GHD, the NRIG, has become the most powerful model (and one of the most popular ones) for financial data. Guo (2017) investigated the adequacy of the GARCH models in risk management on the Hong Kong stock market returns to account for conditional volatility; the author considered the NRIG and compared its empirical performance with the Students t-distribution and the NIG distribution. The results indicated that the NRIG distribution performs slightly better than the other two types of distributions. In Guo (2019), the author further introduced the NRIG distribution to the GJR-GARCH model and also compared its empirical performance with the Students t-distribution and NIG, using a variety of asset return series. Again, he concluded that there is no overwhelmingly dominant distribution in fitting the data under the GARCH framework, although the NRIG distribution performs slightly better than the other two types of distribution. For markets index series, the NRIG distribution with GJR terms improved the models performance, but it is ambiguous for individual stock price series; the author also noted that the GJR-GARCH-NRIG model has practical advantages in quantitative risk management.
Following Guos 2017 framework, there has recently been a growing interest in the NRIG distribution in finance. Several studies also considered the fit of the NRIG distribution in comparison with the Students t-or/and the skewed Students t-distributions. Maree, Carr and Howard 2017 investigated the conditional heavy tails of daily silver spot returns under the GARCH framework; they concluded that the NRIG distribution has the best empirical performance to capture the daily silver spot returns dynamics. Hong, Lee and Ding 2017 also compared the Students t-distribution and the NRIG within the GARCH framework for the daily stock market returns of South Korea (KOSPI). Their findings showed that the NRIG distribution has a better in-sample performance than the Students t-distribution. Oden, Hurt and Gentry 2017 also suggested that the NRIG distribution has superior performance in fitting the German stock market returns compared to the Students t-distribution.
However, the Students t-distribution seemed to overpower the NRIG distribution in some studies. Ma (2017) and Ma et al. (2017) showed that the NRIG cannot outperform the Students t-distribution in modelling the metals spot return. Liu et al. (2017) also suggested that the older Students t-distribution still performs better than the NRIG distribution for modelling the daily Bitcoin exchange rate returns. Additionally, some researchers suggest that the NIG and full GHD provide a better fit to some financial data than an NRIG. For instance, Day et al. (2017) suggest that the NIG has the best empirical performance in fitting the Chinese stock index returns. Kayaba et al. (2017) also took advantage of the results in Guos 2017 study and fitted the market return series through GHDs, namely the Students t-distribution, the skewed t-distribution, the NIG distribution and the full GHD distribution. Their results show that the full GHD has the best fit and generates the most suitable risk measures in modelling the Tokyo stock exchange returns.
To the best of our knowledge, none of the literature has discussed the choice of NRIG distribution and other GHDs to fit the returns of the energy market. This study is therefore considered to become the first study to take into account other subclasses of GHDs, such as HYP, GHST and VG in comparison with the NRIG distribution in risk management for crude oil prices.

Econometric Models
In this section, we present some important theory about the GHD models. The vaule-at-risk and the backtesting procedures are also briefly discussed.
The generalised hyperbolic distributions are the distributions of five parameters that form part of a larger family with the properties of normal variance-mean mixture distributions (Hu and Kercheval, 2007).

MODELLING CRUDE OIL RETURNS USING THE NRIG DISTRIBUTION
The univariate generalised hyperbolic distribution (GHD) can be parametrised in several ways, in this paper we follow Prause (1999) parameterisation. For a random variable X, the probability density of GHD is given by gh(x|λ, α, β, σ, µ) = a(λ, α, β, σ) ( where K λ denotes the modified Bessel function of the third kind with index λ. The domain of variation of the parameters is µ ∈ R and If X has the distribution given in (1), we write X ∼ gh (λ, α, β, σ, µ).
Where µ and σ are the location and scale parameters, respectively. The parameter β determines the skewness, α affects the peakedness and tail behaviour of the distribution, it is related to the kurtosis, whereas λ affects the shape of the distribution and its tail behaviour.

Hyperbolic distribution
This subclass of the generalised hyperbolic distribution occurs when the parameter λ = 1.
Definition A random variable X is said to follow a hyperbolic (HYP) distribution if its probability density function is given by where β is an asymmetrical parameter, σ is a scale parameter, γ = √ α 2 − β 2 and K 1 denotes the modified Bessel function of the third kind with index 1.

Variance gamma distribution
The variance gamma (VG) distribution is a continuous probability distribution that is defined as the normal variance-mean mixture where the mixing density is the gamma distribution. This subclass of the GHD is obtained when the parameter σ → 0.
The parameter domain is also given by λ > 0 and α > |β|. The mean and variance of this distribution are, respectively (see Kemda et al., 2015) .

Generalised hyperbolic skewed Student's t-distribution
The GHST is a limiting case of GHD obtained when λ = − v 2 and α → |β|. It is characterised by having one polynomial and one exponential tail. Following Aas and Haff 2006, a random variable X is said to follow a generalised hyperbolic skewed Students t-distribution (GHST) if its probability density function is given by where Kemda et al. (2015). When β → 0, we have the ordinary Students t-distribution. The mean and variance of a GHST-distributed random variable X are, respectively, The variance is only finite when v > 4, as opposed to the symmetric Students t-distribution, which only requires v > 2. The derivation of the skewness (S) and kurtosis (K) is relatively straightforward due to the normal mixture structure of the distribution. These are given by ] .
The skewness and kurtosis do not exist when v ≤ 6 and v ≤ 8, respectively.

Normal reciprocal inverse Gaussian distribution This definition follows from
) .
The NRIG distribution is normalised by setting which implies E(X) = 0 and V ar(X) = 1. When β = 0, the NRIG distribution is symmetric with excess kurtosis

Risk measures
The amount of asset risk capital reserved by financial institutions under the Basel Agreement is immediately correlated with the portfolio risk level; value-at-risk (VaR) and expected shortfall (ES) are two of the most popular benchmark measures for assessing such risk. ES evaluates the anticipated loss (or gain) value that exceeds the appropriate amount of VaR. The precision of the assessment of VaR and ES therefore depends on how well a chosen model shows extreme data observations McNeil et al. (2005). However, in this study, we only employ the VaR method to assess the risk in the energy market returns.

Value-at-risk (VaR)
For a random variable X (be a profit or loss) with distribution function F over a specified time period, the VaR, for a given probability p, can be defined as the p th quantile of F , given as where F −1 is the quantile function.

Backtesting procedure
The different backtesting procedures are used to examine the adequacy and efficiency of estimates of VaR and ES. In particular, VaR backtesting is carried out using the Kupiec likelihood ratio unconditional coverage test of Kupiec (1995) and the conditional coverage test of Christoffersen (1998).
The backtesting technique used in this study is the Kupiec test, which examines the frequency of losses along the tail. It is based on an understanding that an appropriate model should have its proportion of VaR estimates violations around the corresponding tail probability level Baharul-Ulum et al. (2012).
Under the null hypothesis that the expected proportion of violations is equal to α, the Kupiec statistic is given as and is asymptotically distributed according to a chi-square distribution with one degree of freedom. We therefore reject the null hypothesis if the expected proportion of exceedances x α is lower than LRuc. Where x α is the number of times the observed returns fall below or above the VaR estimate at level α, that is, r t < V aR α or r t > V aR α , and compare the corresponding failure rates to α.

Empirical Results
This section first describes the characteristics of the energy returns and then the robust VaR models are selected for each returns series.

Data and empirical properties
In this paper, we examine the relative performance of the GARCH models combined with GHDs models with regard to evaluation and forecasting VaR in energy market. The data used are the daily closing spot prices of where r t is the daily return for period t, and P t and P t−1 are the daily closing prices for day t and t − 1, respectively. The time series plots in Figure 1 indicate that the daily closing prices have many trends and they do not seem to be stationary in either mean or variance, but the log returns series seem stationary, with a mean value that hovers around zero. However, the variance does not appear to be constant over time, suggesting the pattern of volatility clustering as well as heteroscedasticity, which is highly expected as we are dealing with financial returns. The descriptive statistics, unit root and ARCH tests are reported in Table 2. In Panel B † , we present the results of the Augmented Dickey-Fuller 1979 (ADF), the Phillips-Perron 1988 (PP) unit root tests and the Kwiatkowski, Phillips, Schmidt and Shin 1992 (KPSS) stationarity test. The ADF and PP tests, undoubtedly reject the null hypothesis of unit root (non-stationary) of returns. Furthermore, the KPSS test also reveals that we cannot reject the null hypothesis of stationary at the 5% level of significance. We thus conclude that crude oil and gasoline return time series are stationary. From Panel A, it is obvious that all distributions exhibit fat-tailed and leptokurtic. As indicated by the coefficients of skewness and kurtosis; the Brent, WTI and GCCGR return series are right-skewed, while the NYHCGR presents left-skewed. This evidence is also supported the Q-Q plots in Figure 2. According to the Jarque and Bera (1980) and Shapiro and Wilk (1965) test statistics, we definitely reject the null hypothesis of Gaussian distribution for all returns and conclude that the returns are leptokurtic for the estimation period. Therefore, we propose the use of GHDs to capture the fat-tailedness (non-normality) of the distributions of these returns. Using the LjungCBox Q statistic of order 20 based on the returns, we can also reject the hypothesis of white noise and assert that time series are autocorrelated.  Q(20) is the Ljung-Box Portmanteau statistic. The Q statistic is calculated with 20 lags and is distributed as chi-squared with 20 degrees of freedom. * and ** denote statistically significance at the 5% and 1% level, respectively. The p-values of the KPSS statistics are given in parenthesis.

Parameter estimation and model adequacy test
The maximum likelihood estimation (MLE) method is used to estimate the parameters of the GHDs. The estimates are then used to analyse the goodness-of-fit for these GHDs, and are then compared to the normality assumption. The results of the maximum likelihood estimates of the GHDs and the Anderson-Darling (AD) goodness-of-fit test are presented in Table 3.
The AD goodness-of-fit test provides high p-values at most all the returns, except for HYP and VG for WTI returns (Panel B), suggesting that the GHDs provide a very good fit. The lowest AD statistics suggest the better fit of GHD member with emphasis along the tail, hence the NRIG in Panel A, and the GHST in Panel B, C and D indicate that they fit better compared to other members with an emphasis on the tails.

Value-at-Risk estimates and backtesting
The further analysis of the extreme tails is investigated by using the value-at-risk tool and the Kupiec likelihood ratio test in Tables 4 and 5 to assess the fit of data at the left and right tails.
From Table 4, we can observe that the VaR estimates of GHDs are closer to those of the empirical distribution than to those of the normal distribution at most VaR levels. This is not surprising, as we have seen that the GHDs provide a better fit compared to the normal distribution. In order to check the adequacy of these VaR estimates, the Kupiec likelihood ratio test is used. The p-values of the Kupiec likelihood ratio test at various significance levels are shown in Table 5. Table 5 confirms the inaccuracy of the normal distribution for the VaR estimation; it is rejected by the Kupiec test several VaR levels, especially that upper and lowermost tails. This was anticipated, as it was earlier seen that the normal distribution provided an inadequate depiction of the data set.
From Table 5 (Panel A), it can be seen that almost all GHD models reject the null hypothesis of model adequacy at lower VaR levels in the long position, whereas the models do not reject the null hypothesis of model adequacy at higher levels, except for the GHST distribution. This indicates that most of the VaR models are adequate at high levels in the long position. In short positions, almost all the GHD VaR models fail to reject the null hypothesis of model adequacy at 5%, and at 10% significance level for the VG model. This means that the VaR models are adequate at all levels in a short position. The best model is chosen at various levels by employing the Kupiec likelihood test. The model with the highest p-value at the given level is chosen as the best (robust) model. In the long position, VG distribution has the highest p-values at the 1% level, and the joint highest p-value with HYP distribution at the 5% level; the GHST also has the highest p-value at the 2.5% level. In the short position, the GHST distribution has the highest p-value at the 95% level and the joint highest p-value with the VG and NIG distribution at the 99% level. Finally, the HYP and NRIG distributions have the joint highest p-value at the 97.5% level.
From Table 5 (Panel B), it can be noticed that for the short position, almost all the models rejected the null hypothesis of model adequacy at some VaR levels except the GHST, that is, the NRIG rejected the null hypothesis of model adequacy at the 95% VaR level, whereas the NIG rejected the null hypothesis of model adequacy at the 97.5% VaR level. It can also be observed that the NRIG model has the highest p-values for the long position, except at the 2.5% VaR level, where it is outperformed by the GHST model. The GHST model has the highest p-values for the short positions except at the 95% VaR level, where it is outperformed by the NIG model.
At the 5% level of significance, Table 5 (Panel C) strongly suggests that the normal distribution provides a very poor fit to the NYHCGR returns, characterised by the small p-values at the 1% and 95% VaR levels. This is expected, as the returns are leptokurtic. Interestingly, the table also shows that the p-values for all the fitted GHDs models are greater than 0.05; thus, the null hypothesis of model adequacy is not rejected at the 5% significance level. Using the p-values of the Kupiec likelihood test statistics, the best model is chosen at various levels. The GHST distribution has the highest p-value at the 2.5%, 95% and 97.5% VaR probability levels. The NRIG distribution has the highest p-value at the 1% level, the NIG has the highest p-value at the 5% level, and the HYP and VG have the joint highest p-value at the 99% level.
The results from Table 5 (Panel D) strongly suggest that the normal distribution provides a very poor fit; this is characterised by the small p-values (less than 0.05). The table also shows that the NIG distribution does not fit the tails adequately at the 97.5% level with a low p-value of 0.0396. The p-values for all other fitted GHDs models are greater than 0.05; thus, the null hypothesis of model adequacy is not rejected at all levels under investigation. Using the Kupiec likelihood test statistics p-value, the best model is chosen at various levels. The GHST distribution has the highest p-value at the 2.5%, 95%, and 97.5% VaR probability levels. The VG distribution has the highest p-value at the 1% and 99% VaR levels.

Conclusion
As the volatility in the energy market increases, this market has become risky. It is therefore extremely important to implement an effective risk management system against market risk. In this context, VaR has become the most popular tool for quantifying market risk. This study examined the suitability of using the generalised hyperbolic distributions (GHDs) for modelling the VaR of crude oil and gasoline returns. In particular, our primary objective was to identify an optimal GHD subclass of models for portraying the risks associated with these energy returns, keeping more focus on the NRIG distributions performance. This was done by incorporating 5 GHDs members' methodology.
Both crude oil and gasoline returns were found to be leptokurtic and to exhibit time-varying volatility; hence, the fitting of the GHDs as the heavy-tailed distributions was of importance. The heavy-tail feature of the returns was well captured by the GHDs.
The goodness-of-fit test based on the Anderson-Darling statistics confirmed that the GHDs generally fit the data well, except HYP and VG which didn't fit well at on WTI returns. The performance of the GHDs fit was quite similar, however, overall, the GHST slightly outperformed the other GHDs for WTI, NYHCGR and GCCGR returns, and a marginal performance with the VG for the Brent returns.
The adequacy of the VaR estimate for these GHDs was assessed by in-sample backtesting using the Kupiec likelihood ratio test. Using the Kupiec likelihood test statistics, the best model was chosen at various VaR levels of confidence. The main findings of the study on the risk estimation of crude oil and gasoline returns are summarised in Table 6 below.  In conclusion, the overall fits suggest that the NRIG cannot outperform the VG and GHST in modelling the Brent return, neither the GHST in modelling the WTI and the NYHCGR returns nor VG and NIG and WTI in modelling the GCCGR returns. However, the results from the Kupiec test suggest that the NRIG provides a better VaR estimation at the following VaR levels: • From the Brent returns, the NRIG provides a better marginal fit with the HYP distribution at the 97.5% level. • From the WTI returns, the NRIG provides a better VaR estimate at the 1% and 5% levels.
• From the NYHCGR returns, the NRIG provides a better VaR estimate at the 1% level.
Further work includes the incorporation of a hybrid of the GARCH-type models and the GHD methodology to account for the volatility clustering and the leverage effect characteristics of the returns, and we also recommend an out-of-sample backtest comparison between these GHD models and and also with Generalised Extreme Value Distribution (GEVD) and Generalised Pareto Distribution (GPD) used by Chikobvu and Jakata (2020) for South African financial index. We recommend the review of these models in other emerging markets.