Best Fitting Fat Tail Distribution for the Volatilities of Energy Futures: Gev, Gat and Stable Distributions in GARCH and APARCH Models

Precise modeling and forecasting of the volatility of energy futures is vital to structuring trading strategies in spot markets for risk managers. Capturing conditional distribution, fat tails and price spikes properly is crucial to the correct measurement of risk. This paper is an attempt to model volatility of energy futures under different distributions. In empirical analysis, we estimate the volatility of Natural Gas Futures, Brent Oil Futures and Heating Oil Futures through GARCH and APARCH models under gev, gat and alpha-stable distributions. We also applied various VaR analyses, Gaussian, Historical and Modified (Cornish-Fisher) VaR, for each variable. Results suggest that the APARCH model largely outperforms the GARCH model, and gat distribution performs better in modeling fat tails in returns. Our results also indicate that the correct volatility level, in gat distribution, is higher than those suggested under normal distribution with rates of 56%, 45% and 67% for Natural Gas Futures, Brent Oil Futures and Heating Oil Futures, respectively. Implemented VaR analyses also support this conclusion. Additionally, VaR test results demonstrate that energy futures display riskier behavior than S&P 500 returns. Our findings suggest that for optimum risk management and trading strategies, risk managers should consider alternative distributions in their models. According to our results, prices in energy markets are wilder than the perception of normal distribution. In this regard, regulators and policy makers should enhance transparency and competitiveness in the energy markets to protect consumers.


Introduction
The last few decades have witnessed many financial market upsets, like the 1987 Stock Market Crash, 1997Asian Currency Crisis, 1998LTCM debacle, 2000 Brothers and the greatest of all since the Great Depression-the 2007 Global Financial Crisis. These events have challenged the mainstream view that extreme events in financial markets are outliers with negligible probability (Jackwerth and Rubinstein 1996). Policy makers, practitioners and academicians have thus become interested in studying financial markets, considerably deviating from mainstream financial theory. Thus, over the last many years, modeling and volatility analysis of financial time series has become an important area of research. In addition, financial decisions such as those made by fund managers, option traders and the like portray interest in forecasts of volatility for portfolio optimization, asset pricing and risk minimization. Volatility is the conditional variance of the asset returns and is not precisely and directly measurable, but is a key parameter in many financial applications like derivatives valuation, asset management and risk management (Engle et al. 2007;Ahmed and Zakaria 2011).
In financial time series, where the range of variables is large, the assumption of homoscedasticity might be unreasonable and might lead to a false sense of precision (Engle 1982;Engle et al. 2007).
Numerous models have been proposed on well-documented evidence that asset return follows a geometric Brownian motion (and thus, the concern that the homoscedastic assumption) and is not reliable in financial time series data (Jackwerth and Rubinstein 1996;Sheikh 1991). This has motivated researchers in the field of finance to use changing volatility processes to model asset return. These models can be classified into observation-driven models; for example, Auto-Regressive Conditional Heteroscedastic (ARCH) and Generalized Auto-Regressive Conditional Heteroscedastic (GARCH) proposed by Engle (1982); Bollerslev (1986), respectively, and parameter-driven stochastic volatility (SV) models proposed by Taylor (1986). Engle (1982) was the first to propose a stationary, non-linear model for conditional variance, which evolves according to an autoregressive-type process termed as ARCH model that has been extended to GARCH by Bollerslev (1986); this has been extensively used in financial literature.
These models have been refined from time to time to accommodate the two empirical irregularities of the financial time series. First, the long-known existence of 'fat tails' or 'outlier-prone' probability distribution in the financial time series data. This irregularity is commonly known as leptokurtosis in financial literature. One reason for financial time series data to exhibit this irregularity may be that the conditional variance is not constant, and the outliers appear when the variance is large. GARCH processes exhibit heavy tails even if ε t is Gaussian. Therefore, when we use GARCH models, we can model both the conditional heteroscedasticity and the heavy-tailed distributions of financial time series. GARCH processes reduce leptokurtosis by normalizing the time-varying variances, but can't completely eliminate it, nonetheless. Many financial time series have tails that are heavier than those implied by a GARCH process with Gaussian ε t . The second empirical irregularity is regarding occurrence of volatility clustering where negative returns, unlike equally large positive returns, are followed, by large increases in volatility. This irregularity was first explored by Black (1976) and subsequently by Christie (1982); Schwert (1989), and is now commonly known as the leverage effect. In order to describe this asymmetry, models such as Exponential GARCH-EGARCH (Nelson 1991), Quadratic GARCH-QGARCH (Sentana 1991;Engle 1990), Threshold GARCH-TGARCH (Zakoian 1994), Glosten-Jagannathan-Runkle GARCH-GJR-GARCH (Glosten et al. 1993) and Asymmetric Power ARCH-APARCH (Ding et al. 1993) have been introduced. Asymmetric Power ARCH (APARCH) model is one of the most promising ARCH-type models in a way that it tries to model volatility clustering by nesting at least seven ARCH-type models; this has been found to be effective in capturing several stylized facts and long-memory properties of financial time series, even when the distribution has excess kurtosis (fat tails).
In this study, by considering the aforementioned stylized facts of financial time series, we model the volatility of Natural Gas Futures (NGF), Brent Oil Futures (BOF), and Heating Oil Futures (HGF) through GARCH and APARCH models. Studying volatility of energy commodities for efficient forecasting is important in quantifying the risk due to increased volatility in the energy commodities markets. This is primarily because the energy prices are dynamic and may be explained by factors other than the demand and supply (Aloui and Mabrouk 2010;Nomikos and Andriosopoulos 2012). In volatility modeling, we use different types of heavy tailed distributions unlike most of the literature. In both models, we implement the alpha-stable distribution, generalized extreme values (gev) distribution and generalized asymmetric t (gat) distribution.

Literature Review
Since the pioneering studies of Engle (1982); Bollerslev (1986), a great amount of literature has been developed. Many models have been introduced by the researchers and each of them has some strengths and weaknesses. A large number of studies also have been implemented in energy commodities through different models.
Estimating volatility accurately not only requires modeling mean and variance, but also needs adequate capturing of conditional distribution, fat tails and price spikes. Gibson and Schwartz (1990), Cortazar and Schwartz (1994); Schwartz (1997) show that oil prices exhibit clustering and persistence and thus the assumptions of random walk may be inappropriate to use when studying volatility in energy markets. Cartea and Figueroa (2005) observed that the existence of heavy tails in energy markets is so obvious and thus the probability of extreme events is much higher than estimated by Gaussian distribution.
The existence of asymmetry has been found to play an important role in increasing the volatility of energy commodities. Fan et al. (2008) examined West Texas Intermediate (WTI) and Brent Crude oil prices with various GARCH specifications, and found an asymmetric leverage effect in WTI returns while and not in Brent return. Zhang et al. (2008) attributes this behavior (leverage) mainly to the conduct exhibited by participants in the energy markets as well as to non-renewable properties of oil (energy). In a similar study, Eydeland and Wolyniec (2003) observed leverage in energy markets. Kanamura (2009) examined the same in natural gas prices, while Knittel and Roberts (2005) found leverage in North California hourly electricity prices. Kristoufek (2014) examines Brent, WTI, natural gas and heating oil and finds that the leverage effect is exhibited by natural gas while it is not significant for heating oil; Brent crude oil and WTI exhibit standard leverage effect. Nomikos and Andriosopoulos (2012) employed the EGARCH model to incorporate time varying volatility and found asymmetric leverage effect for gasoline, natural gas, propane, and Interconnection Electricity Firm On Peak Price Index, whereas standard leverage effect was observed for WTI, heating oil and heating oil-WTI crack spread. Sadorsky (2006) estimates the volatility of petroleum futures through alternative models and states that the TGARCH fits better for heating oil and natural gas volatility and the GARCH fits well for crude oil and unleaded gasoline volatility. Kang et al. (2009), argues that CGARCH and FIGARCH models are better able to capture persistence in the volatility of crude oil prices than other GARCH family models. They further argue that CGARCH has a superior ability to forecasting accuracy in crude markets than other models. Agnolucci (2009), on the other hand, argues that among GARCH-type models, implied volatility models do not improve the predictive power of the model, while models assuming errors follow ged distribution and perform better. Mohammadi and Su (2010) argue that APARCH models largely outperform others when modeling and forecasting the conditional mean and volatility of crude markets. Wei et al. (2010) demonstrated that GARCH models, which can control long-memory effects as well as asymmetry, will certainly provide accurate estimates of volatility.
Some studies have analyzed the volatility persistency in energy markets. Cheong (2009) examines WTI and Brent crude oil markets for the evidence of persistence using GARCH specification models and finds that the persistence in WTI is more than Brent crude oil. Kristoufek (2014) examines Brent, WTI, natural gas and heating oil and observes that the volatility is long-term dependent. Alvarez-Ramirez et al. (2002); Serletis and Andreadis (2004) also observed persistence in crude oil volatility. Feng and Shi (2017) show that tempered stable distribution outperforms normal, Student-t and GED distributions. They recommend using this distribution especially in high frequency data and in FIGARCH-type long memory models. Aloui and Mabrouk (2010) confirms fat tails and long memory in energy markets while examining WTI, Brent, New York Harbor conventional gasoline regular (NYHCGR) and Rotterdam conventional gasoline regular (RCGR). They further show that the FIAPARCH model with skewed Student-t distribution provides accurate estimates in predicting one-day ahead VaR. Cartea and Villaplana (2008) demonstrated that electricity markets exhibit a clustering behavior and thus models with time-varying volatility perform better than constant volatility models; Chkili et al. (2014) examines crude oil, natural gas, gold and silver markets and finds that non-linear GARCH specifications over-perform linear GARCH model. Hung et al. (2008) show that the choice of the distribution of the innovation process is crucial for accurately estimating volatility. They demonstrate that heavy tail GARCH models outperform GARCH with t distributions as well as GARCH normal distribution. Gunay (2015) also found that Markov regime switching GARCH (MRS-GARCH) model outperforms all other GARCH models in estimating volatility by accurately capturing the most stylized facts of financial time series. Billio and Pelizzon (2000) found that MRS-GARCH model is superior to GARCH-t models in estimating volatility. Aloui and Jammazi (2009) show that MRS-GARCH can more precisely estimate the volatility in energy markets by substantially capturing volatility clustering and the leverage effect. Giot and Laurent (2003), in their out-of-sample investigation on crude oil, heating oil, propane and conventional gasoline regular show that the skew student APARCH models provides the best estimates of volatility. Huang and Lin (2004) suggests the student-t-APARCH model to consider the nature of distribution in modeling financial time series data. They show that normal APARCH performs better at lower confidence levels, while as student-t-APARCH provides better out-of-sample estimates at a higher confidence level. Bali (2003Bali ( , 2007 used generalized extreme values (gev) distribution for the measurement of risk; the author shows that GJR-GARCH-skewed-t performs better than GARCH-normal in estimating volatility while the gev provides a more robust and precise estimate of risk than models using the skewed-t distribution and normal distribution."

GARCH and APARCH Model
A convenient starting point to talk about volatility modeling is ARCH and GARCH processes. In an ARCH process, Engle (1982) assumes that the variance of the current error term is related to the error terms (weighted average of residuals) of the previous periods. Bollerslev (1986) generalized Engle's (1982) model independently to make it more realistic to fit the real situations, which are presented below as GARCH (p, q) model: where σ 2 t is the conditional variance of y k up to time t -1 and has an autoregressive structure, positively correlated to its past values of squared returns y 2 t . This model provides a systematic way to adjust the parameters ω, α, β to give the best fit, where ω > 0, α i ≥ 0, β j ≥ 0, and α + β < 1 and the innovation sequence {ε t } ∞ i=−∞ is independent and identically distributed with E(ε 0 ) = 0 and E ε 2 0 = 1. Usually in financial time series, large negative returns appear to increase volatility in larger magnitudes than positive returns do; this behavior is known as the leverage effect. Standard ARCH and GARCH models have failed to accommodate this empirical irregularity, because they model σ t as a function of past values of α 2 t without considering whether the past values of α t are positive or negative. A variety of models have been introduced to account for this asymmetry in the financial time series, which include EGARCH (Nelson 1991), TGARCH (Zakoian 1994;Glosten et al. 1993) and APARCH (Ding et al. 1993). The APARCH (p, q) model (Ding et al. 1993) can be defined as: where δ > 0, and −1 < γ < 1. Positive γ i values would indicate that negative shocks have a deeper impact on conditional volatility, whereas a negative γ i would indicate that positive shocks have a larger impact on the current conditional variance (Black 1976). δ plays the role of a Box-Cox transformation of the conditional standard deviation σ t .

Alpha-Stable Distribution
Nolan (2005) describes the stable process as processes that retain their shape under addition; i.e., if random variables X, X 1 , X 2 , . . . X n are independent, symmetrically stable variable, then for every n: Risk Financial Manag. 2018, 11, 30 5 of 19 for some constants C n > 0 and d n . The symbol d = means equality in distribution and the law is called strictly stable if d n = 0 for all n. The class of distributions that satisfy (equation above) are said to be stable and are described by four parametersα, β, γ and δ, which are described as follows: α must be in the range 0 < α ≤ 2 and is known as the index of the law or the index of stability. It determines the weight of the tails and smaller values means greater frequency and size of extreme events. β is known as the skewness parameter and must be in the range −1 ≤ β ≤ 1, β = 0 implies a symmetric distribution, β < 0 implies skewed to the left and β > 0 means the distribution is skewed to the right. γ is a scale parameter or dispersion parameter and must be a positive number that measures dispersion. The parameter δ is a location parameter, δ > 0 implies the right shift in the distribution and δ < 0 implies a left shift in the distribution. The two possible alpha stable parameterizations are described as follows: A random variable X is S(α, β, γ, δ 0 ; 0) if it has the following characteristic function: and a random variable X is S(α, β, γ, δ 1 ; 1) if it has the following characteristic function:

Generalized Extreme Values (GEV) Distribution
The growing interest in modeling shape and thickness of the tails in the financial time series data has led to the development of models that characterize models for extreme events in financial markets. Extreme value theory is one such framework to analyse the behavior of tails in distributions. According to Extreme Value Theory, the asymptotic distribution can be modeled to converge to the Gumble, Frechet, or Wiebull distributions, which could be standardized to generalized extreme values (gev) distribution. It was Jenkinson (1955) who proposed a gev distribution, which included a three-limit distribution-Fréchet, Weibull, and Gumbel distinguished by Gnedenko (1943). Generalized Extreme Values (GEV) is represented below: where ξ is the tail shape parameter and tail index (α) = ξ −1 . ξ = 0 yields thin tail distribution with tail index equal to infinity. At ξ ≥ 0, all moments upto the integer part of the tail index exist and when ξ ≤ 0, all moments are finite or zero. The distribution associated with ξ = 0 belongs to the Gumbel Class, ξ > 0 represents fat tailed Fréchet distributions such as Pareto, Cauchy and Student-t distributions, and ξ < 0 represents Weibull class of distributions. Reiss and Thomas (2001) shows that even for a small positive value of ξ rate of growth of skewness and kurtosis of the distribution results in concentration of probability density at the right tail for the Fréchet distribution. Since the probability of extreme economic losses are higher than for extreme gains, Fréchet distribution is more appropriate to model economic losses (Markose and Alentorn 2011).

Gat Distribution
The Student's t-distribution has been widely used in modelling risk in financial markets. After the introduction of skewed Student's t-distribution for risk modelling by Hansen (1994), several researchers have proposed various extensions of skewed Student's t-distribution. The generalized asymmetric Student-t (AST) skew extension proposed by Zhu and Galbraith (2010) allows for modeling asymmetry by assigning two separate parameters to control the thickness of the right and the left tail, and a third parameter to allow for the skewness to change independently of tail parameters.
The AST distribution is defined by: , and α * =

Empirical Analysis
In this section of the study, we investigate the performance of different conditional distributions (gev, gat and alpha-stable) by means of GARCH and APARCH models. We also implement Gaussian

Unit Root Test
Using stationary time series is one of the requirements for empirical analyses in the domain of robustness and validation of results. Broadly speaking, stationarity is a desirable statistical feature of a time series with regards to its constant mean and variance over time. Using non-stationary time series might bring on spurious results in analysis. Although many tests are introduced in literature, the model proposed by Kapetanios (2005) is quite robust and takes structural breaks into account. This test investigates stationarity of time series against the occurrence of an unspecified number of breaks up to 5. Results of the Kapetanios (2005) multiple break unit root test is presented in Table 1. According to the test statistics, we reject the null hypothesis of the model under all confidence levels, which indicates that each variable is stationary within the consideration of structural breaks.

Normality Tests
Normality in return distributions is one of the primary assumptions of many financial models. However, many studies have already showed that financial time series exhibit departures from normal distribution and display fat tails and high kurtosis features. In addition, as stated by Mandelbrot (1972), financial time series displays long memory features. These features are named in literature as stylized facts of financial time series. Taleb (2007) states that normal distribution underestimates the probability of extreme events in financial markets. Therefore, any financial model using normal distribution offers very low chance to events in tail probabilities and thus might cause failures in risk management. Table 2 shows the normality test results of NGF, BOF and HOF returns. To test normality, we use Anderson-Darling, Cramer-VonMises, Lilliefors and Pearson Chi-Square methods, which test if the time series comes from a certain probability distribution (such as normal distribution). The p values in each test statistic indicate that null hypothesis for all return series is rejected. Accordingly, it is can be argued that return series exhibits a substantial departure from the normal distribution. In addition to the normality tests, we plot autocorrelation graphs for the squared returns of each variable. High autocorrelation in squared returns imply long memory in returns of the time series. It means that there is a persistent relationship even between distant observations. This persistence, which is called long memory feature in returns, can be linked to alpha-stable distributions that consider the fat tails and heavy kurtosis in probability distributions of returns. As stated by Panas (2001) alpha-stable distributions do not only correspond to the distribution of returns, but is also a useful instrument in detection of long memory. The relationship between the Hurst exponent, which is the measure of long memory, and alpha parameter in alpha-stable distributions can be addressed as follows: α = 1/H. Any H value between 0.5 and 1, which corresponds to alpha value between 1 and 2, is the signature of long memory. On the other hand, H value between 0 and 0.5 and alpha value greater than 2 is the indicator of short memory or mean reverting in time series. As a preliminary evaluation for the long memory feature, we plot the squared return autocorrelations together with two-sided 5% critical values in Figure 1. According to the test statistics, we reject the null hypothesis of the model under all confidence levels, which indicates that each variable is stationary within the consideration of structural breaks.

Normality Tests
Normality in return distributions is one of the primary assumptions of many financial models. However, many studies have already showed that financial time series exhibit departures from normal distribution and display fat tails and high kurtosis features. In addition, as stated by Mandelbrot (1972), financial time series displays long memory features. These features are named in literature as stylized facts of financial time series. Taleb (2007) states that normal distribution underestimates the probability of extreme events in financial markets. Therefore, any financial model using normal distribution offers very low chance to events in tail probabilities and thus might cause failures in risk management. Table 2 shows the normality test results of NGF, BOF and HOF returns. To test normality, we use Anderson-Darling, Cramer-VonMises, Lilliefors and Pearson Chi-Square methods, which test if the time series comes from a certain probability distribution (such as normal distribution). The values in each test statistic indicate that null hypothesis for all return series is rejected. Accordingly, it is can be argued that return series exhibits a substantial departure from the normal distribution. In addition to the normality tests, we plot autocorrelation graphs for the squared returns of each variable. High autocorrelation in squared returns imply long memory in returns of the time series. It means that there is a persistent relationship even between distant observations. This persistence, which is called long memory feature in returns, can be linked to alpha-stable distributions that consider the fat tails and heavy kurtosis in probability distributions of returns. As stated by Panas (2001) alpha-stable distributions do not only correspond to the distribution of returns, but is also a useful instrument in detection of long memory. The relationship between the Hurst exponent, which is the measure of long memory, and alpha parameter in alpha-stable distributions can be addressed as follows: = 1/ . Any value between 0.5 and 1, which corresponds to alpha value between 1 and 2, is the signature of long memory. On the other hand, value between 0 and 0.5 and alpha value greater than 2 is the indicator of short memory or mean reverting in time series. As a preliminary evaluation for the long memory feature, we plot the squared return autocorrelations together with two-sided 5% critical values in Figure 1. As can be seen from the autocorrelation graphs, there is significant autocorrelations for each squared returns series up to 100 lags at 95% confidence level and especially in the NGF autocorrelations decay at a hyperbolic rate. Even though the level of the autocorrelations is not high, they remain positive in a substantial number of lags. Autocorrelations are more pronounced and persistence in the BOF. These findings demonstrate that there might be a long memory feature in squared return series. Therefore, it can be concluded that using alpha-stable distributions in volatility modeling would be useful since it considers the long memory and fat tails in return distributions. As stated by Nolan (2005), stable distributions consider skewness and fat tails and also have very interesting mathematical features; they can also fit financial data better because they consider stylized facts of financial time series. In Table 3, we present the estimated alpha-stable distribution parameters ( , , and ) of three returns series. The first parameter of the alpha-stable distributions is the stationary index (characteristic exponent) and theoretically has values between 0 and 2. Results show that for all returns series, alpha parameter is lower than 2, which indicates that distribution of return series follow a power law. As stated by Weron (2001), a smaller alpha value indicates thicker tails in the distribution. This situation reveals that the BOF returns are subject to more extreme returns in comparison to HOF and NGF returns. The second parameter denotes skewness and has values between −1 and 1; for = 0, the distribution is symmetrical. In cases where > 0, the distribution is skewed right and < 0 indicates that the distribution is skewed left. While in the NGF and HOF returns value is larger than 0, BOF has values lower than 0. Thus, while NGF and HOF returns are skewed right, BOF returns are skewed left. Thus, it can be said that for NGF and HOF returns the right tail of the As can be seen from the autocorrelation graphs, there is significant autocorrelations for each squared returns series up to 100 lags at 95% confidence level and especially in the NGF autocorrelations decay at a hyperbolic rate. Even though the level of the autocorrelations is not high, they remain positive in a substantial number of lags. Autocorrelations are more pronounced and persistence in the BOF. These findings demonstrate that there might be a long memory feature in squared return series. Therefore, it can be concluded that using alpha-stable distributions in volatility modeling would be useful since it considers the long memory and fat tails in return distributions. As stated by Nolan (2005), stable distributions consider skewness and fat tails and also have very interesting mathematical features; they can also fit financial data better because they consider stylized facts of financial time series. In Table 3, we present the estimated alpha-stable distribution parameters (α, β, γ and δ) of three returns series. The first parameter of the alpha-stable distributions α is the stationary index (characteristic exponent) and theoretically has values between 0 and 2. Results show that for all returns series, alpha parameter is lower than 2, which indicates that distribution of return series follow a power law. As stated by Weron (2001), a smaller alpha value indicates thicker tails in the distribution. This situation reveals that the BOF returns are subject to more extreme returns in comparison to HOF and NGF returns. The second parameter β denotes skewness and has values between −1 and 1; for β = 0, the distribution is symmetrical. In cases where β > 0, the distribution is skewed right and β < 0 indicates that the distribution is skewed left. While in the NGF and HOF returns β value is larger than 0, BOF has β values lower than 0. Thus, while NGF and HOF returns are skewed right, BOF returns are skewed left. Thus, it can be said that for NGF and HOF returns the right tail of the distribution is longer than left side and the frequency of negative returns is higher than the positive returns. On the contrary, in BOF, positive returns have a higher frequency than negative returns.

Volatility Modeling-GARCH (1.1) and APARCH (1.1)
In this section of the study, we model the variance through GARCH and APARCH for the return series. As normal distribution fails to capture most of the stylized facts of time series such as fat tails and excess kurtosis in returns, based on the results of Table 2, alternative return distributions (gev, gat and alpha-stable distributions) have been used to model volatility. The most commonly used method in volatility modeling literature is the GARCH model of Bollerslev (1986) which allows for a much more flexible lag structure than the ARCH model of Engle (1982). Unlike the original ARCH model, GARCH processes use lagged conditional variances in modeling besides the lagged squared error terms. Although several versions of GARCH family models exist in literature, in this study we prefer to use APARCH with original GARCH, as it models volatility by taking into account long memory and leverage features of asset returns. Results of the GARCH model can be found in Table 4. Results show that for all three return series, α and β parameters are statistically significant under all types of distributions. Besides, regarding the tail behavior of return distributions, both skewness (λ) and shape (γ 1 and γ 2 ) parameters of gev, gat and alpha-stable distributions are statistically significant under all return distributions, except for the λ parameter in the HOF. According to AIC, BIC and AICc, for the NGF and the HOF returns the most successful results are obtained under alpha-stable distribution. BOF, however, performs better under gat distribution. The sum of α and β coefficients is close to unity for all type of distributions, which, according to Bollerslev (1986) implies extreme persistence in volatility. However, when alpha-stable distribution is incorporated into the GARCH model, persistence in variance comparatively decreases to 0.98. Persistence in conditional volatility reveals that besides the GARCH model, utilization of APARCH model would be significant in modeling volatility.    Figure 2 below presents conditional volatility graphs of return distributions. The results clearly indicate that for all return series, volatility is significantly higher in alpha-stable distribution comparing to other distributions. Provided that the alpha-stable distribution is the accurate distribution, the highest volatility is obtained for NGF returns. Second and third highest distributions are exhibited by HOF and BOF returns (0.0858, 0.0591, and 0.0586 respectively). While the NGF returns have high volatility throughout the period we analyze, BOF and HOF returns is very volatile, especially in the following periods: 1990, 1998-2003-2016. J. Risk Financial Manag. 2018 Figure 2 below presents conditional volatility graphs of return distributions. The results clearly indicate that for all return series, volatility is significantly higher in alpha-stable distribution comparing to other distributions. Provided that the alpha-stable distribution is the accurate distribution, the highest volatility is obtained for NGF returns. Second and third highest distributions are exhibited by HOF and BOF returns (0.0858, 0.0591, and 0.0586 respectively). While the NGF returns have high volatility throughout the period we analyze, BOF and HOF returns is very volatile, especially in the following periods: 1990, 1998-2003, 2008-2009 and 2015-2016 Figure 2. Conditional Volatility of Return Series.
As shown by Cont (2001), stylized facts of financial time series are vital to modeling. Therefore, deficiencies in the first step of risk management and misreading the data might cause obstacles for further policies. Failure of the LTCM is one such example. Most discussed stylized facts in literature are fat tails in return distribution along with leverage effect and long memory in returns. In order to take into account these facts, in this study we use APARCH and GARCH models with alternative fat tail distributions. As stated before, the APARCH model introduced by Ding et al. (1993) provides some advantages to Bollerslev's (1986) GARCH model through its flexibility. Ding et al. (1993) shows that rather than quadratic returns in the GARCH model, other powers of conditional As shown by Cont (2001), stylized facts of financial time series are vital to modeling. Therefore, deficiencies in the first step of risk management and misreading the data might cause obstacles for further policies. Failure of the LTCM is one such example. Most discussed stylized facts in literature are fat tails in return distribution along with leverage effect and long memory in returns. In order to take into account these facts, in this study we use APARCH and GARCH models with alternative fat tail distributions. As stated before, the APARCH model introduced by Ding et al. (1993) provides some advantages to Bollerslev's (1986) GARCH model through its flexibility. Ding et al. (1993) shows that rather than quadratic returns in the GARCH model, other powers of conditional residuals display significant temporal dependencies in the APARCH model (Sáfadi and Alencar 2014).
Results obtained in Table 5 demonstrate that all α and β parameters are statistically significant in the APARCH model at minimum 95% confidence level. When we look at the alternative distributions, we see that the best-fitting distribution type for the data is the gat distribution based on the AIC, BIC and AICc. According to the results of the gat distribution, both skewness (λ) and shape (γ 1 and γ 2 ) parameters are statistically significant for all three variables, while as for BOF returns, only leverage effect (γ) is significant with 0.1505 value. As stated by Laurent (2004), a positive leverage means that past negative returns affect today's volatility more than past positive returns. That is, past bad news has more prominent effect on volatility than good news. The second important parameter in the APARCH estimations is δ; it represents a Box-Cox power transformation of the conditional standard deviation process. According to the results, δ parameter, which is equal to 2 as default in the GARCH model, is statistically significant in gat distribution for all return series with values 0.5787 for the NGF, 1.2334 for the BOF and 1.1864 for the HOF. AIC, BIC and AICc indicate that comparing with the previous results, the APARCH model is superior to the GARCH model in all return series, except for alpha-stable distribution in NGF and HOF returns. residuals display significant temporal dependencies in the APARCH model (Sáfadi and Alencar 2014). Results obtained in Table 5 demonstrate that all and parameters are statistically significant in the APARCH model at minimum 95% confidence level. When we look at the alternative distributions, we see that the best-fitting distribution type for the data is the gat distribution based on the AIC, BIC and AICc. According to the results of the gat distribution, both skewness ( ) and shape ( and ) parameters are statistically significant for all three variables, while as for BOF returns, only leverage effect ( ) is significant with 0.1505 value. As stated by Laurent (2004), a positive leverage means that past negative returns affect today's volatility more than past positive returns. That is, past bad news has more prominent effect on volatility than good news. The second important parameter in the APARCH estimations is ; it represents a Box-Cox power transformation of the conditional standard deviation process. According to the results, parameter, which is equal to 2 as default in the GARCH model, is statistically significant in gat distribution for all return series with values 0.5787 for the NGF, 1.2334 for the BOF and 1.1864 for the HOF. AIC, BIC and AICc indicate that comparing with the previous results, the APARCH model is superior to the GARCH model in all return series, except for alpha-stable distribution in NGF and HOF returns.      According to Figure 3, conditional variance graphs for the APARCH model demonstrate that in all return series the highest volatility values are obtained under gat distribution. As discussed before, the best fitting distribution for returns was gat distribution in the APARCH model. The lowest variance values, on the other hand, are attained for the gev distribution. In terms of the average values when we compare the volatilities in the gat distribution, the highest volatility in the decreasing order is exhibited by NGF, HOF, and BOF returns. Recalling GARCH model results, NGF exhibited largest volatilities among all commodities. One other important information provided in conditional variance figures is with regard to the fluctuation in specific time periods. Although the NGF returns have largest fluctuations on an average, the BOF returns present severe responses to the developments occurring in 1990, 2008-2009, and 2015. The analyses so far conducted demonstrate that Natural Gas Futures, Brent Oil Futures and Heating Oil Futures are quite risky and Natural Gas Futures display wild behavior compared to others. Besides, normal distribution in volatility modeling underestimates the level of risk. In order to further investigate these findings, and provide more evidence on this context, we calculate the Value at Risk for each variable. In conjunction with our current variables, we also use S&P 500 index returns as a benchmark. VaR analyses are performed through three different models: Gaussian VaR, which is based on the normal distribution, Historical VaR and Modified (Cornish-Fisher) VaR, which takes fat tails into account in return distributions. Results obtained are presented in Table 6.  Table 6 support the previous findings in GARCH and APARCH analysis. Accordingly, it is seen that NGF is still the most-risky variable among all. Additionally, all these energy futures display riskier behavior than S&P 500 returns. For example, VaR NGF is almost three times greater than the Value at Risk of S&P 500 returns. Similarly, VaR BOF and VaR HOF are approximately two times greater than VaR S&P500 . These findings indicate that investors should exercise extreme care and caution in their hedging strategies, while taking positions in energy commodities.

Results in
The graphical representation of VaR results for each variable is presented in Figure 4 above. In these graphs, the horizontal axis displays the confidence level and the vertical axis represents the Value at Risk. A similar pattern is observed in each variable; while Gaussian and Historical VaR track very close to each other under respective confidence levels, Modified VaR scatters in a wider range against changes in the confidence level. Figure 3, conditional variance graphs for the APARCH model demonstrate that in all return series the highest volatility values are obtained under gat distribution. As discussed before, the best fitting distribution for returns was gat distribution in the APARCH model. The lowest variance values, on the other hand, are attained for the gev distribution. In terms of the average values when we compare the volatilities in the gat distribution, the highest volatility in the decreasing order is exhibited by NGF, HOF, and BOF returns. Recalling GARCH model results, NGF exhibited largest volatilities among all commodities. One other important information provided in conditional variance figures is with regard to the fluctuation in specific time periods. Although the NGF returns have largest fluctuations on an average, the BOF returns present severe responses to the developments occurring in 1990, 2008-2009, and 2015.

According to
The analyses so far conducted demonstrate that Natural Gas Futures, Brent Oil Futures and Heating Oil Futures are quite risky and Natural Gas Futures display wild behavior compared to others. Besides, normal distribution in volatility modeling underestimates the level of risk. In order to further investigate these findings, and provide more evidence on this context, we calculate the Value at Risk for each variable. In conjunction with our current variables, we also use S&P 500 index returns as a benchmark. VaR analyses are performed through three different models: Gaussian VaR, which is based on the normal distribution, Historical VaR and Modified (Cornish-Fisher) VaR, which takes fat tails into account in return distributions. Results obtained are presented in Table 6.  Table 6 support the previous findings in GARCH and APARCH analysis. Accordingly, it is seen that NGF is still the most-risky variable among all. Additionally, all these energy futures display riskier behavior than S&P 500 returns. For example, VaRNGF is almost three times greater than the Value at Risk of S&P 500 returns. Similarly, VaRBOF and VaRHOF are approximately two times greater than VaRS&P500. These findings indicate that investors should exercise extreme care and caution in their hedging strategies, while taking positions in energy commodities.
The graphical representation of VaR results for each variable is presented in Figure 4 above. In these graphs, the horizontal axis displays the confidence level and the vertical axis represents the Value at Risk. A similar pattern is observed in each variable; while Gaussian and Historical VaR track very close to each other under respective confidence levels, Modified VaR scatters in a wider range against changes in the confidence level.

Conclusions
In this study, volatility modeling was performed for some energy futures: Natural Gas Futures (NGF), Brent Oil Futures (BOF) and Heating Oil Futures (HOF). In order to test whether return series obey normality assumption or not, we conducted Anderson-Darling, Cramer-Von Mises, Lilliefors and Pearson Chi-Square test and found that all return series exhibited a significant departure from normal distribution. Therefore, we estimate the parameters of alpha-stable distribution and observe that fat tail distributions for the returns series would be more suitable in empirical modeling of volatility. Accordingly, in conjunction with the normal distribution, we employed gev, gat and alpha-stable distribution as well. In addition, as a preliminary test through autocorrelation graphs, we checked the persistence of volatility of returns, and since the results exhibit high long memory, we implemented the APARCH model in conjunction with the conventional GARCH model. Our results show that the APARCH model fits the returns series better than the GARCH model in most cases. Specifically, the APARCH model outperforms the GARCH model in BOF returns and performs better under gat distribution. According to the AIC, BIC and AICc information criteria, the APARCH model fits better in the NGF and HOF returns (normal, gev and gat distribution, and the GARCH model fits the data slightly better for alpha-stable distribution.

Conclusions
In this study, volatility modeling was performed for some energy futures: Natural Gas Futures (NGF), Brent Oil Futures (BOF) and Heating Oil Futures (HOF). In order to test whether return series obey normality assumption or not, we conducted Anderson-Darling, Cramer-Von Mises, Lilliefors and Pearson Chi-Square test and found that all return series exhibited a significant departure from normal distribution. Therefore, we estimate the parameters of alpha-stable distribution and observe that fat tail distributions for the returns series would be more suitable in empirical modeling of volatility. Accordingly, in conjunction with the normal distribution, we employed gev, gat and alpha-stable distribution as well. In addition, as a preliminary test through autocorrelation graphs, we checked the persistence of volatility of returns, and since the results exhibit high long memory, we implemented the APARCH model in conjunction with the conventional GARCH model. Our results show that the APARCH model fits the returns series better than the GARCH model in most cases. Specifically, the APARCH model outperforms the GARCH model in BOF returns and performs better under gat distribution. According to the AIC, BIC and AICc information criteria, the APARCH model fits better in the NGF and HOF returns (normal, gev and gat distribution, and the GARCH model fits the data slightly better for alpha-stable distribution. The highest volatility values are obtained in the NGF returns for gat distribution, while the BOF displays the lowest level of fluctuation. Additionally, when we look at the conditional variance graphs for the APARCH model, it is seen that the gat distribution, which is the best-fitting tail distribution based on AIC, BIC and AICs information criteria, exhibits higher volatility than other models. With regard to models with normal distribution, it is observed that conditional variance levels are lower than the volatility in the gat distribution model on an average. For the volatilities of NGF, BOF and HOF returns, normal distribution presents 56%, 45% and 67% lower levels, respectively, than the best-fitting distribution, which is gat. It can thus be stated that normal distribution underestimates the true risk level. Furthermore, Gaussian, Historical and Modified (Cornish-Fisher) VaR analyses indicate that energy commodities exhibit very high risk, and are even riskier than the stock market, based on our benchmarking analysis in S&P 500 returns.
These findings demonstrate the importance of accurate methods in the measurement of risk. When normal distribution is assumed instead of gat distribution, the level of risk might be underestimated, and thus lead to valuation issues in risk management instruments, which by and large still continue to operate on normal distribution. Therefore, it can be said that researchers and practitioners are relying more on normal distribution, which has a tendency to underestimate the level of risk. To conclude, our results indicate that the level of risk in energy markets is higher than what is predicted by the estimates of conventional models, and energy markets are riskier than the stock markets as indicated by the VaR analysis. In this regard, investors should be more careful in formulating their hedging strategies. In addition, by considering this wild volatility, policy makers might take actions to protect customers against exposures to fluctuations. Beyond investment strategies, high fluctuation in energy prices also brings instability in capital budgeting plans of institutions. Unforeseen price changes could lead to negative net present value in project valuations, and if left unaddressed, the consequences, of course, could be drastic to different participants in the economy.