ARMA–GARCH model with fractional generalized hyperbolic innovations

In this study, a multivariate ARMA–GARCH model with fractional generalized hyperbolic innovations exhibiting fat-tail, volatility clustering, and long-range dependence properties is introduced. To define the fractional generalized hyperbolic process, the non-fractional variant is derived by subordinating time-changed Brownian motion to the generalized inverse Gaussian process, and thereafter, the fractional generalized hyperbolic process is obtained using the Volterra kernel. Based on the ARMA–GARCH model with standard normal innovations, the parameters are estimated for the high-frequency returns of six U.S. stocks. Subsequently, the residuals extracted from the estimated ARMA–GARCH parameters are fitted to the fractional and non-fractional generalized hyperbolic processes. The results show that the fractional generalized hyperbolic process performs better in describing the behavior of the residual process of high-frequency returns than the non-fractional processes considered in this study.

successfully capture volatility clustering but sometimes fail to adequately account for conditional leptokurtosis. Hence, a need for more flexible conditional distributions to be explored has arisen. In this sense, much research has tried to account for the fat-tail property found in most financial return data: for example, the Student's t-distribution (Bollerslev 1987), non-central t-distribution (Harvey and Siddique 1999), generalized hyperbolic skew Student's t-distribution (Aas and Haff 2006), Johnson's SU-distribution (Johnson 1949;Yan 2005;Rigby and Stasinopoulos 2005;Choi and Nam 2008;Naguez and Prigent 2017;Naguez 2018), generalized error distribution (GED) (Nelson 1991;Baillie and Bollerslev 2002), and Laplace distribution (Granger and Ding 1995). However, these distributions generally have numerous issues in the applications, even though they may be capable in some cases, and the issues include that a GED distribution does not have a sufficient fat tail to account for extreme events and the density of a non-central t-distribution has the sum of an infinite series.
In addition, during the past decades, there has been growing interest in other types of fat-tail distributions in finance, including a generalized hyperbolic distribution (Barndorff-Nielsen 1977), a normal inverse Gaussian distribution (Barndorff-Nielsen 1995), a variance gamma distribution (Madan and Seneta 1990), a stable distribution (Mandelbrot 1963a;Liu and Brorsen 1995;Panorska et al. 1995;Mittnik et al. 1998;Mittnik and Paolella 2003), and a tempered stable distribution (Rosiński 2007;Kim et al. 2009). In particular, the generalized hyperbolic distribution is a very general form that involves other types of distributions as special cases, such as a Student's t-distribution, hyperbolic distribution, normal inverse Gaussian distribution, variance gamma distribution, and Laplace distribution (Anderson 2001;Jensen and Lunde 2001;Venter and De Jongh 2002;Forsberg and Bollerslev 2002;Stentoft 2008). Although the aforementioned fat-tail distributions have been studied for empirical innovations in ARMA-GARCH models to capture the fat-tail property on financial return data, Sun et al. (2008) and Beck et al. (2013) evoke that these distributions are unable to fully describe the stylized facts, especially in high-frequency time series analysis.
However, since Mandelbrot and Van Ness (1968) have introduced fractional Brownian motion 2 to explain short-and long-range dependence, several empirical studies on asset prices (Lo 1991;Cutland et al. 1995;Willinger et al. 1999;Robinson 2003;Casas and Gao 2008;Chronopoulou and Viens 2012;Cont 2005;Mariani et al. 2009;Cheridito 2003;Comte and Renault 1998;Rosenbaum 2008) have investigated the long-range dependence properties of asset returns. In the sense that the fractional Brownian motion can capture the long-range dependence properties, the ARCH and GARCH models with innovations utilizing the fractional Brownian motion can capture not only the long-range dependence but also the volatility clustering effect on high-frequency return data. However, the fat-tail properties remain unconsidered. Hence, in an attempt to better describe high-frequency returns, Sun et al. (2008) offer a univariate model with three stylized facts by taking the ARMA-GARCH model with fractional stable distributed residuals, and they suggest that the model performs effectively in high-frequency returns. In addition, Kim (2012) introduced the fractional multivariate normal tempered stable process as a subclass of the fractional stable process 3 using time-changed fractional Brownian motion with the fractional tempered stable subordinator. The fractional tempered stable process was redefined by Kim (2015) with the stochastic integral for the Volterra kernel presented in Houdre and Kawai (2006), and it was applied to the innovations on the multivariate ARMA-GARCH model 4 . In addition, Kozubowski et al. (2006) applied fractional Laplace motion by subordinating fractional Brownian motion to a gamma process to model a financial time series. The fractional normal inverse Gaussian process was proposed in Kumar et al. (2011) and Kumar and Vellaisamy (2012).
Much research has addressed wide issues, especially the simulation and parameter estimations in processes with long-range dependence, as well as the detection of longrange dependence. This study is motivated by the availability of high-frequency time series return data and the need for a new market model accommodating three stylized facts. The main contribution of this study is the development of a general fractional model that incorporates three stylized facts into the volatility process. This is the first study to explore the fractional generalized hyperbolic process and apply it to a multivariate ARMA-GARCH model. The generalized hyperbolic process is derived using time-changed Brownian motion with the generalized inverse Gaussian subordinator. Thereafter, the univariate fractional generalized hyperbolic process is defined by the stochastic integral for the Volterra kernel. The basic properties and covariance structure of the multivariate fractional generalized hyperbolic process are discussed, and a numerical method for generating the sample paths of the processes is presented. Using highfrequency stock return data, parameter estimations on the ARMA-GARCH model with fractional generalized hyperbolic innovations are provided, and goodness-of-fit tests are performed for the estimated parameters to validate the model.
The purpose of this study is threefold. First, a new ARMA-GARCH model that accounts for three stylized facts is presented to provide the broad framework of an alternative model that can be used in financial economic applications. Second, we examine whether considering three stylized facts, especially with a fractional generalized hyperbolic innovation, produces superior forecast estimates and thereby motivates their use in effectively managing financial risk as well as optimizing portfolios in high-frequency time series return data. Finally, this study calls for attention to an assumption underlying the mean-variance model (Markowitz 1952) applied to high-frequency returns, that is, the assumption that asset returns are normally distributed. By being the first to derive a fractional generalized hyperbolic process using time-changed Brownian motion, we contribute to the empirical literature on modeling for financial risk and portfolio management. In the portfolio optimization based on Markowitz's mean-variance model, the Gaussian assumption can be replaced by the ARMA-GARCH model with fractional generalized hyperbolic innovations, and the portfolio value-at-risk (VaR) and average value-at-risk (AVaR) based on the model can supersede the variance risk measure.
The remainder of this paper is organized as follows: in the next section, Models and Methodology, the notion of long-range dependence is discussed, time-changed 3 Refer to Samorodnitsky and Taqqu (1994) for stable process. 4 The fractional tempered stable process is also discussed on the option pricing by Kim et al. (2019).
Brownian motion is reviewed, the generalized hyperbolic process is derived from the time-changed Brownian motion, and the fractional generalized hyperbolic process as well as the covariance structure for two elements of the processes are discussed. Section Simulation presents the simulation of fractional generalized hyperbolic processes and an illustration of the sample paths for representative parameter values. The multivariate ARMA-GARCH model with long-range dependence is defined in Section ARMA-GARCH Model with fGH Innovations, with empirical illustration. The main findings and implications are summarized in Section Conclusion.

Models and methodology
This section reviews the background knowledge regarding long-range dependence and the generalized hyperbolic process. Subsequently, a fractional generalized hyperbolic process is derived.

Long-range dependence
Long-range dependence in a stochastic process is typically defined based on its autocorrelation functions. While a short memory stationary process usually refers to a stochastic process whose autocorrelation functions decay fast and whose spectral density is bounded everywhere, a long memory process 5 is associated with the slow decay of autocorrelation functions as well as a certain type of scaling linked to self-similar processes. Although the notion of long-range dependence varies from author to author, the commonly used notion is recalled in terms of the autocorrelation function of a process. For lag k, a stationary process with autocorrelation function ρ(k) is referred to as a process with long-range dependence if does not converge. According to this notion, absolute autocorrelations are not integrable. Thus, a weakly stationary process is said to have long-range dependence if ρ(·) hyperbolically decays, that is, ρ(k) ∼ L(k)k 2d−1 as k → ∞ , 0 < d < 1 2 , where L(k) is a slowly varying function as x → ∞ , if for every positive constant c, exists and is equal to 1. Hence, the memory length in a stochastic process is measured by the rate at which its correlations decay with a lag. To illustrate, suppose an autocorrelation function ρ(·) of an ARMA process at lag k. The function converges rapidly to zero as k → ∞ in the sense that there exists r > 1 such that In this case, the correlations for the ARMA processes decay exponentially fast with k. Conversely, consider the autocorrelation function ρ(k) of the FARIMA(p, d, q) process with 0 < |d| < 1 2 . This autocorrelation function has the following property: indicating that ρ(k) converges to zero as k → ∞ at a much slower rate than that of an ARMA process. Therefore, ARMA processes are regarded as having short memory, whereas FARIMA processes are regarded as exhibiting long memory. Likewise, the longrange dependence property depends on the behavior of the autocorrelation function at large lags. However, it may be difficult to empirically estimate the long-range dependence property (Beran 1994). Thus, self-similar processes are often used to formulate models with the long-range dependence property. A stochastic process {X(t)} t≥0 is said to be self-similar if there exists H > 0 such that for any scaling factor c > 0 , the processes {X(t)} t≥0 and {c H X(t)} t≥0 have the same law: where H is referred to as the self-similarity exponent of the process X. For a continuous stochastic process x t , we assume that Thereafter, we have: By the self-similarity property, the correlation is Based on L'Hopital's rule, Thus, lim n→∞ k−n k=−n ρ(k) exist if and only if H < 1 2 , and H > 1 2 implies that the series does not converge; thus, X(t) is a long-memory process.

Time-changed Brownian motion
By substituting the usual deterministic time t as subordinator T = {T (t)} t≥0 in a stochastic process Y = {Y (t)} t≥0 , we obtain a new process X = {X(t)} t≥0 , and Y is said to be subordinated by process T. Namely, Hence, we change the time to the stochastic process Y to run on a "new clock" whose stochastic time is dominated by the subordinator T. Intuitively, it can be considered as the business time (or trading time) that can be viewed as stochastic. If we employ an arithmetic Brownian motion subordinated by a subordinator, we obtain the time-changed Brownian motion. More specifically, considering that a subordinator T = {T (t)} t≥0 is independent of the standard Brownian motion {B(t)} t≥0 , we have a new process X = {X(t)} t≥0 , which is called the time-changed Brownian motion, by substituting physical time t as the subordinator T as follows: Let F T (t) be the filtration of T. We can derive the characteristic function φ X(t) via subordination theorem (Cont and Tankov 2004;Barndorff-Nielsen and Levendorskii 2001;Sato 1999). That is to say, where φ T (t) is the characteristic function of the subordinator T(t). A Lévy process can be constructed by subordinating Brownian motion with a particular subordinator (Clark 1973). Therefore, subordination is used as a constructive tool to define a particular class of Lévy processes. In the next section, a generalized hyperbolic process is presented by applying time-changed Brownian motion.

Generalized hyperbolic process
The density function of the generalized inverse Gaussian distribution described by the three parameters ( , δ, γ ) is given by where Y is the modified Bessel function of the second kind with index as follows: The parameter domain of the generalized inverse Gaussian distribution is

Proposition 1 The characteristic function of a generalized inverse Gaussian random variable G is given by
Proof See Appendix.
Using the same approach as in the proof of Proposition 1, the rth moment can be derived as follows, given a generalized inverse Gaussian random variable G with parameters ( , δ, γ ), where q( , δ, γ ) represents the norming constant of the corresponding generalized inverse Gaussian density. In this study, we set the expected value of the generalized inverse Gaussian random variable G to one, that is, E[G] = 1 . Using the moment formula in Equation (2), the expected value is defined by and set Then, we have a generalized inverse Gaussian random variable G under parametrization ( , α) . It follows identity Hence, the variance of G is given by Owing to the infinite divisibility of a generalized inverse Gaussian distribution (Barndorff-Nielson and Halgreen 1977), the characteristic function of a generalized inverse Gaussian process In the remainder of this paper, we assume that N is a positive integer that represents the dimension, and β = (β 1 , β 2 , · · · , β N ) ⊤ , θ = (θ 1 , θ 2 , · · · , θ N ) ⊤ , while � = [σ m,n ] m,n∈{1,2,··· ,N } is a correlation matrix such that σ n,n = 1 for n ∈ {1, 2, · · · , N } . Thereafter, we consider an N-dimensional Brownian motion B = {B(t)} t≥0 , such that B(t) = (B 1 (t), B 2 (t), · · · , B N (t)) ⊤ , and suppose that cov(B m (t), B n (t)) = σ m,n t for all m, n ∈ {1, 2, · · · , N } and that B and the generalized inverse Gaussian process {G(t)} t≥0 with parameters ( , α) are independent. In addition, consider the N-dimensional process , with X(t) = (X 1 (t), X 2 (t), · · · , X N (t)) ⊤ . For n ∈ {1, 2, · · · , N } , we define {X(t)} t≥0 by the time-changed Brownian motion as Thereafter, process X is referred to as the N-dimensional generalized hyperbolic process with parameters ( , α, β, θ, �) . We define the characteristic function of a generalized hyperbolic process X n (t) in the following proposition.

Proposition 2
The characteristic function of a generalized hyperbolic process X n (t) is given by Furthermore, the covariance between X m (t) and X n (t) is The variance of X n (t) is equal to
Proposition 4 For m, n ∈ {1, 2, · · · , N } , the covariance between Z m (t) and Z n (t) is equal to Proof See Appendix.
Let L be the lower triangular matrix obtained by the Cholesky decomposition for with = L L T , where is the correlation matrix in Equation (3). Thus, we obtain a mutually independent vector of Brownian motion B (t) = (B 1 (t),B 2 (t), · · · ,B N (t)) ⊤ .
com, Inc., and CVS Health Co. for 7 trading days from February 5 to February 13, 2020, are obtained 7 . The estimated ARMA(1,1)-GARCH(1,1) parameters µ n , a n , b n , ω n , ξ n , and    Table 3 Estimated parameters of non-fractional GH and fractional GH process  ζ n with standard normal innovations are presented in Table 1. The residuals, which are obtained from the estimated parameters, are used to compute the Hurst index H n for n = 1, 2, · · · , N , and they are reported in Table 2, with the estimated H = 0.5387 as the mean of H n . The estimated fractional generalized hyperbolic parameters H, , α , β , and θ are listed in Table 3 along with the estimated non-fractional generalized hyperbolic parameters. Note that the ARMA-GARCH model with fractional generalized hyperbolic innovations becomes the ARMA-GARCH model with non-fractional generalized hyperbolic innovations when parameter H = 1/2 . The estimated matrix is given in Table 4. Three criteria are used to compare the goodness-of-fit of the innovation processes under consideration: the Kolmogorov-Smirnov (KS) test statistics, log-likelihood (LL), and Akaike information criterion (AIC). The KS test is performed based on the statistic: where F (x) and F(x) are the empirical sample and estimated theoretical distributions, respectively. The KS statistics are calculated for the fractional generalized hyperbolic distribution with estimated parameters ( , α, β n , θ n ) and the empirical sample distribution of {X n (t k ) − X n (t k−1 )} for k = 1, 2, · · · , M , where X n (t k ) is the extracted process using Equation (8). The p-values of the KS statistic 8 represent the amount of evidence available to invalidate the null hypothesis that the data are drawn from the concerned theoretical distribution. Hence, the smaller the p-values, the more the evidence against the null hypothesis. The LL is an overall measure of the goodness-of-fit, implying that higher values are more likely to be the distribution candidates for modeling the data. The AIC is a measure of the goodness-of-fit that estimates the relative support for a model, which is where w is the number of parameters in the calibrated model and L is the maximized value of the likelihood function of the estimated model. Table 5 presents the goodness-of-fit tests of the innovation processes of each model. The ARMA-GARCH model with fractional generalized hyperbolic innovations is not rejected by KS statistics at the 10% significance level for all the six stock returns, which means the empirical sample distribution is not statistically different from the estimated theoretical distribution. However, the ARMA-GARCH models with nonfractional and normal innovations are rejected at the 10% significance level for all considered stocks. Indeed, the fractional generalized hyperbolic models have the highest LL as well as the smallest AIC for all the innovation processes analyzed. Furthermore, to assess the power of the goodness-of-fit tests for a small sample size, 390 observed 1-minute stock prices of the six companies are taken, and these are the prices on February 5, and the same estimation/test procedures are followed. The inferences from Table 5 are consistent as shown in Table 6. Consequently, these results suggest that the ARMA-GARCH model with fractional generalized hyperbolic innovations can be used as an approximation to the empirical sample distribution. Figure 4 exhibits the distributions of the residuals from the models for each stock. Therefore, it is noteworthy that the ARMA-GARCH model with fractional generalized hyperbolic innovations better describes the behavior of the high-frequency stock returns considered in this study than the ARMA-GARCH models with non-fractional generalized hyperbolic and normal innovations. In this regard, portfolio VaR and AVaR based on the ARMA-GARCH models with fractional generalized hyperbolic innovations can be applied to the portfolio optimization. As alternative risk measures, the VaR and AVaR can overcome the limitations of the use of portfolio variance as a measure of risk. In portfolio optimization that portfolio managers and investors use for portfolio rebalancing, a crucial measure is the marginal risk contribution to the portfolio AIC = 2w − 2 log L,  (Gourieroux et al. 2000), which is the risk measure with respect to a given portfolio holding. Hence, by employing the VaR and AVaR calculated by the ARMA-GARCH model with fractional generalized hyperbolic innovations, we can find the optimal portfolio. Moreover, we obtain encouraging results from the goodness-of-fit test that it can be applied to an option pricing model equipped with volatility clustering, leverage effect and conditional skewness, leptokurtosis, and long-range dependence.

Conclusion
In this study, the non-fractional generalized hyperbolic process is derived by subordinating the time-changed Brownian motion to the generalized inverse Gaussian process, and the fractional generalized hyperbolic process is obtained using the Volterra kernel. Thereafter, we introduce the multivariate ARMA-GARCH model with fractional generalized hyperbolic innovations that exhibit three-stylized facts: fat tail, volatility clustering, and long-range dependence properties. This model is compared with ARMA-GARCH models with non-fractional generalized hyperbolic and normal innovations. The fractional generalized hyperbolic process performs better in describing the behavior of the residual process of high-frequency returns than the non-fractional generalized hyperbolic process or Gaussian process. Although the results reveal some important insights into high-frequency return data, the models By the property of the generalized hyperbolic process X n , we have Hence, we obtain From Equation (5) and (6), Refer to Prause (1999), Mercado (2011), andEberlein andHammerstein (2004) in more details on the properties of the generalized hyperbolic process.