Tail risk forecasting using Bayesian realized EGARCH models
TTail risk forecasting using Bayesian realized EGARCH models
Vica Tendenan, Richard Gerlach, and Chao WangDiscipline of Business Analytics, Business School,The University of Sydney, Australia.
Abstract
This paper develops a Bayesian framework for the realized exponential generalizedautoregressive conditional heteroskedasticity (realized EGARCH) model, which canincorporate multiple realized volatility measures for the modelling of a return series.The realized EGARCH model is extended by adopting a standardized Student- t and astandardized skewed Student- t distribution for the return equation. Different types ofrealized measures, such as sub-sampled realized variance, sub-sampled realized range,and realized kernel, are considered in the paper. The Bayesian Markov chain MonteCarlo (MCMC) estimation employs the robust adaptive Metropolis algorithm (RAM)in the burn in period and the standard random walk Metropolis in the sample period.The Bayesian estimators show more favourable results than maximum likelihood es-timators in a simulation study. We test the proposed models with several indices toforecast one-step-ahead Value at Risk (VaR) and Expected Shortfall (ES) over a periodof 1000 days. Rigorous tail risk forecast evaluations show that the realized EGARCHmodels employing the standardized skewed Student- t distribution and incorporatingsub-sampled realized range are favored, compared to a range of models. Keywords : realized EGARCH; Markov chain Monte Carlo; realized volatility; Valueat Risk; Expected Shortfall.
1. Introduction
Value at Risk (VaR) has become a standard measure in assessing market risks since theintroduction of the RiskMetric model by J.P. Morgan in 1993. It is a measure of themaximum potential loss of certain portfolio over a specified holding period at a givenconfidence level and has been widely employed by corporations and financial institutionsin risk management due to its simplicity and convenience in computation. However, theVaR suffers from several weaknesses. First, it provides no information about the size ofthe losses beyond a particular quantile. Second, it is considered as an incoherent measureof market risks since it violates subadditivity property (McNeil et al., 2005).Expected Shortfall (ES) or Conditional VaR was proposed by Artzner et al. (1997,1999) as a coherent risk measure to overcome the shortcomings of VaR. It measures theconditional expectation of the return that is higher than a VaR threshold. It can beviewed as an extension of VaR by exploring further details of the tail distribution and1 a r X i v : . [ q -f i n . R M ] A ug herefore can provide a measure about the magnitude of the loss of tailed events, whichthe VaR cannot provide (Yamai and Yoshiba, 2005). Since May 2012, Basel Committeeon Banking Supervision (2012) has recommended the use of ES to better capture the risksof fat-tailed events. The Basel III Accord, which was implemented in 2019, places newemphasis on ES: “ES must be computed on a daily basis for the bank-wide internal modelsto determine market risk capital requirements. ES must also be computed on a daily basisfor each trading desk that uses the internal models approach (IMA).” (Basel Committeeon Banking Supervision, 2019).The calculation of parametric VaR and ES requires an accurate volatility estimation.The parametric approach includes the use of volatility models that estimate conditionalquantiles by using the resulted forecasts of conditional volatility. The popular models tocapture volatility dynamics are the Autoregressive Conditional Heteroskedasticity (ARCH)model (Engle, 1982) and the Generalized ARCH model (Bollerslev, 1986). The next de-velopment of GARCH-class models includes the EGARCH model (Nelson, 1991) and theGJR-GARCH model (Glosten et al., 1993) that capture the asymmetry in volatility inresponse to positive and negative returns. These early models utilize daily returns inmodelling the conditional variance, which provides noisy signals about the current level ofvolatility (Hansen et al., 2012). In addition, these models have a slow adjustment to volatil-ity movements, since they rely on long and slowly decaying weighted moving averages ofpast squared returns (Andersen et al., 2003). With the availability of high-frequency data,more efficient realized measures of volatility have been proposed, such as realized variance(Andersen and Bollerslev, 1998; Andersen et al., 2003), realized range (Martens and vanDijk, 2007; Christensen and Podolskij, 2007), and realized kernel (Barndorff –Nielsen andShephard, 2002; Barndorff –Nielsen et al., 2008).The realized GARCH model was proposed by Hansen et al. (2012) as an extension ofthe GARCH-type models, which takes advantage of the availability of the high-frequencydata by jointly model the returns and realized measures of volatility. The key featureof the realized GARCH model is an additional equation, called the measurement equa-tion, which contemporaneously relates the conditional variance with a realized measure.Recently, Hansen and Huang (2016) proposed the realized EGARCH model to improvethe realized GARCH model by allowing the use of multiple realized volatility measuresand facilitating a more flexible modeling of the dependence between returns and volatil-ity. Another extension of the realized GARCH model, called threshold realized GARCHmodel, was proposed by Chen and Watanabe (2019) by incorporating high-frequency real-ized volatility measures into a threshold framework to model different volatility behaviorsacross regimes.Error distribution specification in the volatility models is also an important aspect ofthe estimation of the parametric tail risks. Early financial time series modelling based onARCH and GARCH type models were largely based on the assumption of normal errors,while vast evidence have shown that financial time series are highly leptokurtic and often-times skewed. For instance, Bollerslev (1987) allowed for conditional t -distributed errorsin ARCH and GARCH models; Hansen (1994) and Fernndez and Steel (1998) proposed2 skewed Student- t distribution; Chen et al. (2012) combined an asymmetric Laplace dis-tribution (ALD) with a GJR-GARCH model; Chen and Gerlach (2013) developed a moreflexible extension of the ALD, a two-sided Weibull distribution, as the conditional returndistribution in four GARCH-type models.The original realized GARCH model (Hansen et al., 2012) and the realized EGARCHmodel (Hansen and Huang, 2016) both assume Gausiannity of both observation and mea-surement equation errors. While it is adequate to adopt the normal distribution for mea-surement equation errors, the study of Watanabe (2012) and Contino and Gerlach (2017)on the realized GARCH model suggests that the choice of observation error distribution isfound to be important, which in most cases the skewed Student- t distribution is favored.Wang et al. (2019) also extended the realized GARCH model by proposing the two-sidedWeibull distribution developed by Chen and Gerlach (2013) for the observation equationerrors.The common method to estimate ARCH and GARCH type models is by using theMaximum Likelihood (ML) estimation. However, the ML estimation faces several chal-lenges. First, since the optimization procedure is usually encumbered by the positivityconstraints and stationarity condition, the algorithm can be less stable and produces lessconsistent results, which will then create problems in standard error calculation (Silva-pulle and Sen, 2011). This problem might be more severe for complex models such as therealized EGARCH with multiple realized measures. Second, for the realized EGARCHmodels, the asymptotic analysis and properties of the ML estimators are complicated andhave not been well studied. For simplification, Hansen and Huang (2016) used a scorefunction and the numerical Hessian matrix of the log likelihood function to compute thestandard errors for the estimators.It is argued that the application of the Bayesian approach in the GARCH-type mod-elling will lead to a better inference and prediction compared to the classical ML ap-proaches. First, the parameter constraints and the the covariance stationarity conditionin the GARCH-type models can be incorporated by specifying appropriate priors. Second,the inference of the GARCH-type model parameters is straightforward in the Bayesian set-ting. Appropriate Markov chain Monte Carlo (MCMC) procedures will allow us to obtainthe posterior distribution of any function of the GARCH-type model parameters. Theissue of local maxima convergence or convergence to incorrect values which is usually en-countered in the ML estimation of complex GARCH-type models can be avoided by usingthe MCMC techniques. Third, the Bayesian estimation could provide more reliable andconsistent results even for small finite samples (Ardia and Hoogerheide, 2010; Virbick-aite et al., 2015). Fourth, Bayesian approach presents a natural way to assess parameteruncertainty in volatility estimation that can be used in decision making, for instance fore-casting purposes. The predictive distributions of the future volatilities and the measureof precision for risk measures can be obtained via predictive interval (Paap, 2002; Aus´ınand Galeano, 2007). Relevant to this, the predictive distribution of risk measures such asthe VaR is of paramount importance in risk management (Jorion, 2001).This paper focuses on forecasting the VaR and ES based on the extension of the realized3GARCH model of Hansen and Huang (2016). The extension of the realized EGARCHis conducted in several ways. First, two types of return distributions, i.e. a standardizedStudent- t and a standardized skewed Student- t distribution, are proposed, to take intoaccount the non-Gaussianity in the error distribution. The measurement equation errorsare kept to be Normal distribution. The two frameworks are called RE-SkN and RE-tN.Second, three types of realized measures are utilized in the realized EGARCH model,either individually or jointly. Those are the sub-sampled realized variance (RVSS), sub-sampled realized range (RRSS), and realized kernel (RK). The sub-sampling process isbased on Zhang et al. (2005) to deal with micro-structure effects. The choice of theRVSS and RRSS is motivated by the recent works of Gerlach and Wang (2016) andWang et al. (2019) who found that the RVSS and RRSS could improve the out-of-sampleforecasting of their proposed models. Meanwhile, the RK is chosen since it has somerobustness to market microstructure effects (Barndorff –Nielsen et al., 2008) and has beenconsistently used in the original realized GARCH and realized EGARCH models. Third,an adaptive MCMC procedures is adopted and applied to the proposed models. To assessthe performance of the proposed RE-SkN and RE-tN models, their 2.5% and 1% (followingthe recommendation of Basel Committee (2012; 2016; 2019)) VaR and ES forecasts arecompared to those of various competing models such as the original realized GARCH andrealized EGARCH models, the Conditional Autoregressive Expectile (CARE) model byTaylor (2008), and several extensions of GARCH-type models.The structure of the paper is as follows: Sections 2 reviews the employed realizedmeasures, Section 3 discusses the proposed extension to the realized EGARCH models, andSection 4 discusses the Bayesian approach to estimate the proposed models. In Section 5,we present the result of a simulation study using the proposed method and models. Section6 provides an overview of tail risks forecasting methods and relevant model evaluations.Section 7 presents empirical study results, followed by the conclusion of the paper andpossible extensions in Section 8.
2. Realized Measures
This section reviews three types of realized measures employed in this study, i.e. therealized variance (RV), realized range (RR), and realized kernel (RK). For day t , let H t , L t and C t denote the intra-day high, low, and closing prices, respectively. The dailyclose-to-close log return is calculated as: r t = log( C t ) − log( C t − ) , where r t is the associated volatility estimator.In a high frequency intra-day framework, each day t from open to close can be dividedinto N intervals of equal size with length (cid:52) , where each intra-day time is subscripted as i = 0 , , , ..., N . Let P t − i (cid:52) denotes the closing price at the i - th interval of day t . RV is4alculated as the sum of the N intra-day squared returns, at frequency (cid:52) , for day t : RV (cid:52) t = N (cid:88) i =1 (cid:2) log( P t − i (cid:52) ) − log( P t − i − (cid:52) ) (cid:3) . (1)The RR, proposed by Martens and van Dijk (2007) and Christensen and Podolskij(2007), may contain more information about volatility than the RV, in the same way thatthe intra-day range contains more information than squared returns, i.e. it uses all theprice movements in a time period, not just the price at the start and end of each sub-period. Let H t,i = sup ( i − (cid:52) RV with T /N sub-sampling isgiven by: RV T/M,T/N = (cid:80) n k − i =0 RV i n k . (5)Let H t,i = sup ( i + n k ( m − (cid:52) 1) and i + n k m , respectively.Wang et al. (2019) proposed the T /M RR with T /N sub-sampling that takes the followingform: RR i = M (cid:88) m =1 ( H t,i − L t,i ) ; where i = 0 , , ..., n k − . (6) RR T/M,T/N = (cid:80) n k − i =0 RR i n k . (7)For details regarding the properties of the sub-sampled RR, please refer to Wang et al.(2019).The selection of sampling frequency needs to take into account a trade-off betweenbias and variance, e.g. see Bandi and Russell (2008) and Zhang et al. (2005). The RV andRR are usually computed from intra-day returns sampled at a moderate frequency, such6s 5-minute frequency, to balance the bias and variance. This study employs 5-minute RVand RR with 1 min sub-sampling proposed by Wang et al. (2019), as follows: RV , , = (log C t − log C t ) + (log C t − log C t ) + ...RV , , = (log C t − log C t ) + (log C t − log C t ) + ...RV , = (cid:80) i =0 RV , ,i ,RR , , = (log H t ≤ t ≤ t − log L t ≤ t ≤ t ) + (log H t ≤ t ≤ t − log L t ≤ t ≤ t ) + ...RR , , = (log H t ≤ t ≤ t − log L t ≤ t ≤ t ) + (log H t ≤ t ≤ t − log L t ≤ t ≤ t ) + ...RR , = (cid:80) i =0 RR , ,i . The third realized measure proposed to be employed in this study is the realized kernel(RK), which is somewhat robust to noise. The RK was developed by Barndorff –Nielsenet al. (2008) and was also employed the original realized GARCH and realized EGARCHmodels. The estimator takes the following form: K ( X ) = H (cid:88) h = − H k (cid:16) HH + 1 (cid:17) γ h , where γ h = n (cid:88) j = | h | +1 x j,t x j −| h | ,t . (8)The kernel weight function, k ( x ), in equation (8), is chosen to be the Parzen weightfunction, since it satisfies the smoothness conditions, k (cid:48) (0) = k (cid:48) (1) = 0. This will ensurea non-negative estimate. The parzen kernel function is defined as follows: k ( x ) = − x + 6 x ≥ x ≥ / , − x ) / ≥ x ≥ , x > , (9)where x j is the high-frequency return in equation (8) and (9), H is a standard bandwidthspelt out in Barndorff –Nielsen et al. (2009), given by: H ∗ = c ∗ ξ / n / , c ∗ = ((12) / . / = 3 . , ξ = ω T (cid:82) T σ u du . For details, please see Barndorff –Nielsen et al. (2009). 3. Proposed model The realized EGARCH model of Hansen and Huang (2016) with K realized measurestakes the following form: r t = µ + (cid:112) h t (cid:15) t , log h t = ω + β log h t − + τ ( (cid:15) t − ) + γ (cid:48) u t − , log x k,t = ξ k + ϕ k log h t + δ ( k ) ( (cid:15) t ) + u k,t , k = 1 , ..., K, (10)7here r t is the daily return, (cid:15) t iid ∼ D (0 , 1) and u t iid ∼ D (0 , Σ), (cid:15) t and u t = ( u ,t , ..., u K,t ) (cid:48) are mutually and serially independent, x t is the vector of realized measures, and τ ( (cid:15) ) and δ ( k ) ( (cid:15) ) are leverage functions given by τ ( (cid:15) ) = τ (cid:15) + τ ( (cid:15) − 1) and δ ( k ) ( (cid:15) ) = δ k, (cid:15) + δ k, ( (cid:15) − , k = 1 ..., K. , respectively.The realized EGARCH model in equation (10) consists of three types of equations.The first equation is the standard return equation. The second equation is the GARCHequation and the third type is the measurement equation(s). There are two main featuresof the GARCH equation in the realized EGARCH model. The first feature is the pres-ence of a leverage function, τ ( z t − ) that distinguishes the realized EGARCH model fromthe realized GARCH model, which only includes a leverage function in the measurementequation. The second key feature is the last term, γ (cid:48) u t − , which facilitates the chan-nelling of the impact of realized measure on the expectations of future volatility. Since u t is K-dimensional, it allows us to incorporate multiple realized measures of volatility.The parameter β represents volatility persistence, whereas γ reflects the information therealized measures contain about future volatility. The measurement equation(s) show therelationship between the ( ex-post ) realized measures of volatility and the ( ex-ante ) condi-tional variance. In their paper, Hansen and Huang (2016) used the RK, the daily range,and six types of realized variances and adopt a Gaussian specification, D (0 , 1) = N (0 , D (0 , Σ) = N (0 , Σ).This study extends the realized EGARCH model by allowing D (0 , 1) to be a stan-dardized Student- t , and a standardized skewed Student- t of Hansen (1994), followingWatanabe (2012) and Contino and Gerlach (2017) for the case of the realized GARCH.These proposed models are denoted hereafter as RE-tN and RE-SkN, respectively, whilethe original realized EGARCH model employing normal distribution for both return andmeasurement equations is denoted as RE-NN.The Student- t distribution could be potentially restrictive, since it only allows a varia-tion in the location, scale, and tail fatness. Meanwhile, the skewed Student- t distributionis more flexible and allows for a richer behaviors. It is able to model not only leptokur-tosis, but also asymmetry (Hansen, 1994). The pdf of a standardized skewed Student- t distribution is:( ε | ν, λ ) = bc (cid:18) ν − (cid:16) bε + a − λ (cid:17) (cid:19) − ( ν +1) / if ε < − ab ,bc (cid:18) ν − (cid:16) bε + a λ (cid:17) (cid:19) − ( ν +1) / if ε ≥ − ab , (11)where ν is the degrees of freedom parameter that controls the kurtosis of the distributionwith range 2 < ν < ∞ , λ is the skewness parameter with range − < λ < 1, while a, b, and c are some constants given by: a = 4 λc (cid:18) ν − ν − (cid:19) , b = (cid:112) λ − a , c = Γ (cid:16) ν +12 (cid:17)(cid:112) π ( ν − (cid:16) ν (cid:17) . Hansen (1994) showed that this is a proper density function that has zero mean and8 unit variance and specializes to a Student- t density by setting λ = 0. When λ > 0, themode of the density is to the left of zero and the variable is right-skewed, and vice versawhen λ < h t = ω − γ (cid:48) ξ k + ( β − γ (cid:48) ϕ k ) log h t − + a t , (12)where a t = τ ( (cid:15) t − ) + γ (cid:48) (log x k,t − − δ ( k ) ( (cid:15) t )).Taking the expectation of both sides of equation (12), we derive the long run value oflog h as given by: E (log h ) = ω − γ (cid:48) ξ k − ( β − γ (cid:48) ϕ k ) . (13)The stationarity of the realized EGARCH model in equation (10) is then achieved byimposing the following restriction: β − γ (cid:48) ϕ k < . (14) 4. Bayesian Estimation This section discusses the Bayesian estimation approach and Markov chain Monte Carlo(MCMC) sampling procedures employed in estimating the proposed models and conduct-ing the tail-risk forecasting. Basically, the Bayesian model is a probability model thatconsists of a likelihood function p ( y | θ ) and a prior distribution π ( θ ). The product of thelikelihood and the prior probability is equivalent to the posterior distribution p ( θ | y ). The log-likelihood function for the RE-SkN type models, where D ∼ Skt ν,λ (0 , 1) and D = N (0 , Σ) is given by: (cid:96) ( r ; θ ) = T [ (cid:96) ( b ) + (cid:96) ( c )] − (cid:80) Tt =1 (cid:34) v +12 log (cid:32) v − (cid:32) b(cid:15) t + a − λ (cid:33) (cid:33) + 0 . h t )] if (cid:15) t < − ab ,T [ (cid:96) ( b ) + (cid:96) ( c )] − (cid:80) Tt =1 (cid:34) v +12 log (cid:32) v − (cid:32) b(cid:15) t + a λ (cid:33) (cid:33) + 0 . h t )] if (cid:15) t ≥ − ab .(cid:96) ( x ; θ, Σ) = − n (cid:88) t =1 K log(2 π ) + log( | Σ | ) + u (cid:48) t Σ − u t ] ,(cid:96) ( r, x ; θ, Σ) = (cid:96) ( r ; θ ) + (cid:96) ( x ; θ, Σ) . (15)9eanwhile, the log-likelihood function for the RE-tN type model, where D ∼ t ν (0 , D = N (0 , Σ) takes the following form: (cid:96) ( r, x ; θ, Σ) = T (cid:104) log Γ (cid:16) v + 12 (cid:17) − log Γ (cid:16) v (cid:17) − . v − π ] (cid:105) − T (cid:88) t =1 (cid:34) v + 12 log (cid:32) r t − µ ) ( v − h t (cid:17) + 0 . h t (cid:105) − n (cid:88) t =1 (cid:104) K log(2 π ) + log( | Σ | ) + u (cid:48) t Σ − u t (cid:105) . (16)The parameters in the realized EGARCH models are: θ = ( µ, ω, β, τ (cid:48) , γ (cid:48) , ψ (cid:48) , ..., ψ (cid:48) K , ν, λ ) (cid:48) and Σ , (17)where ψ k = ( ξ k , ϕ k , δ (cid:48) k ) (cid:48) , k = 1 , ..., K. Joint Conditional Likelihood is given by: L ( θ, Σ , r, x ) = T (cid:89) t =1 p ( r t , x t |F t − , θ, Σ) (18)= T (cid:89) t =1 p ( r t |F t − , θ ) p ( x t | r t , F t − , θ, Σ) (19)where F t = ( r , ..., r t ; x , ..., x t ) (cid:48) ; the predictive density p ( r t | F t − , θ ) is determined by thedistributional choice of (cid:15) t ; and the conditional density p ( x t | r t , F t − , θ, Σ) is determined bythe distribution of u t . A combination of mostly uninformative and Jeffreys-type priors is chosen over the possibleregion for the parameters: π ( θ ) ∝ σ u,k =1 σ u,k =2 .. σ u,k = K ν I ( A ) , (20)where A is the region described by the parameter restriction in Table 1. The indicatorfunction I ( A ) = 1 over the region A and zero otherwise.10able 1. Parameter RestrictionParameter Bounds µ, ω, τ , τ , ξ, δ , δ ±∞ σ u > ν (4 , λ ( − , β − γ (cid:48) ϕ < u t in Equation (10). The priorof the degrees of freedom parameter ( ν ) is ν − , which is an improper prior obtained bybeing flat on ν − (Bauwens et al., 2006). This is equivalent to U (0 , ) so that ν > 4. Thisensures that the first 4 moments of the distribution of (cid:15) t are finite (Chen et al., 2006).Other research works have adopted this type of prior on σ u and ν on Bayesian estimationof the previous model, the realized GARCH models, e.g. see Gerlach and Wang (2016),Contino and Gerlach (2017), and Wang et al. (2019). In order to numerically perform an inference on the model parameters, a MCMC procedureis performed to draw samples from the posterior distribution of the model parameters.An adaptive MCMC method of Contino and Gerlach (2017) is adopted and extended,by considering the balancing between the target acceptance rate and convergence, thescaling factor for the tuning process, and the selection of the block of the parameters inthe updating process. The adaptive MCMC approach consists of two steps, burn-in periodand post burn-in period.In the burn-in period, a Robust Adaptive Metropolis (RAM) algorithm of Vihola(2012) is employed. The RAM algorithm estimates the shape of the target distributionand simultaneously allow us to attain a given mean acceptance rate. The RAM is definedrecursively through the following process:1. Compute the proposal Y n := X n − + S n − U n , where S n − is a lower diagonal ma-trix with positive diagonal elements and U n is a vector of random variable from d -dimensional Gaussian distribution.2. The proposal is accepted with probability α n := min (cid:110) , π ( Y n ) π ( X n − ) (cid:111) .3. If the proposal Y n is accepted, set X n := Y n . Otherwise, set X n := X n − .4. Compute the lower diagonal matrix, S n , satisfying the equation: S n S Tn = S n − (cid:16) I + η n ( α n − α ∗ ) U n U Tn || U n || S Tn − (cid:17) , (21)where (cid:8) η n (cid:9) n ≥ ⊂ (0 , 1] is a step size sequence decaying to zero, often defined as η n = n − γ with an exponent γ ∈ (1 / , α ∗ ∈ (0 , 1) is the target mean acceptance11robability, and I ∈ R dxd is the identity matrix.To obtain S n in step (4), we use rank-one Cholesky update of S n − when ( α n − α ∗ ) > S n − when ( α n − α ∗ ) < 0. This approach is computationally moreefficient than computing the right hand side of equation (21) and performing new Choleskydecomposition. Folowing Vihola (2012), η n = min { , d · n − / } , where d is the dimensionof the block of parameters. The target mean acceptance rate α ∗ = 0 . 234 (if d > ≤ d ≤ d = 1), as recommended by Roberts et al. (1997). Fortheoretical details of the RAM algorithm, see Vihola (2012).The algorithm is run for 20,000 iterations for each model parameter to achieve conver-gence. The first half of burn-in iterates (10,000 iterates) is discarded, then the remainingsecond half is used to calculate the sample mean ( X n ) and covariance matrix (Σ) for eachparameter as the proposal distribution in the post burn-in period. In the post burn-inperiod, a standard RWM algorithm is employed and a mixture of three Gaussian proposaldistribution is used as follows: Y n ∼ N (cid:0) X n − , C i Σ (cid:1) , (22)where C = 1; C = 100 , C = 0 . 01, following Wang et al. (2019). The algorithm is runfor 10,000 iterations for each model parameter, which are then used to calculate the tailrisk forecasts. The posterior mean of the tail risk forecasts is used as the final forecast.The convergence of the adaptive MCMC method is assessed via the comparison of theestimated variances between-chains and within-chains for each model parameter by usingthe Gelman-Rubin diagnostic. The convergence is indicated by small differences betweenthese variances (Gelman et al., 2013). The efficiency of the MCMC method is evaluatedby using an effective sample size testing (Gelman et al., 2013) and the autocorrelationtime (Chib and Jeliazkov, 2001).To determine the setting of the parameter blocks, two settings of parameter blocks arecompared: 1 block (as in Contino and Gerlach (2017)) and 4 blocks. The latter settingis adapted from Wang et al. (2019) and adjusted to the proposed models in this paper asfollows: θ = ( µ ), θ = ( ω, β, τ (cid:48) , γ (cid:48) , ϕ k ), θ = ( ξ k , δ (cid:48) k , Σ), and θ = ( ν, λ ). The parameterswithin the same parameter block are more likely to be correlated, in the posterior (orlikelihood), than those between blocks. For example, the stationarity condition causescorrelation between iterates of β, γ (cid:48) , ϕ k . Therefore, those parameters are are grouped inthe same block. 5. Simulation study A simulation study is conducted to compare the performance of the MCMC estimators andmaximum likelihood (ML) estimators of the proposed RE-SkN model. The purpose is todemonstrate the superior performance of the MCMC estimators in terms of unbiasednessand precision properties.A total of 1000 replicated datasets are simulated from the proposed RE-SkN model12sing 1 realized measure, with a sample size of 2000. r t = 0 . (cid:112) h t (cid:15) t (23)log h t = − . 12 + 0 . 98 log h t − − . (cid:15) t − + 0 . (cid:15) t − − 1) + 0 . u t − log x t = − . 17 + 0 . 94 log h t − . (cid:15) t + 0 . (cid:15) t − 1) + u t , where (cid:15) t ∼ Skt ν =4 . ,λ =0 . , u t ∼ N (0 , . h = 0 . ν = 5, for bothestimation methods. The tail risks are considered for the quantile level α = 2 . 5% and α = 1%, following the recommendation of Basel Committee (2012; 2016; 2019).First, to assess the convergence of the MCMC method, m = 10 independent MCMCchains are run for two different parameter block settings. Then, we analyze the variancebetween and within the chains according to the Gelman-Rubin diagnostic (Gelman et al.,2013). Table 2 presents the summary of the the Gelman-Rubin diagnostic ˆ R , the numberof effective sample size n eff (Gelman et al., 2013), and the autocorrelation time ˆ t (Chib andJeliazkov, 2001) for both 1 parameter block and 4 parameter blocks settings. Favorableresults are indicated by the values printed in bold. The ˆ R for each parameter is very closeto 1 and all ˆ R are < . 1, the recommended threshold by Gelman et al. (2013), exceptfor ν under the 1 block setting. The ˆ R under the 4 blocks settings is closer to 1 for 11out of 13 parameters, suggesting an improved convergence property under this parameterblock setting. The convergence can also be assessed visually by using the trace plots ofthe MCMC iterations of those independent chains, as depicted for γ and σ u parameters inFigure 1. As can be seen, the MCMC chains show a rapid convergence and chain mixing.The result also suggests that the number of effective sample size for most parameters, n eff , is higher than the suggested benchmark (5 m ) by Gelman et al. (2013). The numberof effective sample size for 10 parameters under the 4 blocks setting is higher than thatunder the 1 block setting. As regards the correlation time ˆ t , the parameters under the4 blocks setting also have a more favorable property, i.e. have a smaller autocorrelationtime. Overall, the MCMC method under the 4 parameter blocks setting provides a muchmore improved convergence property, is more efficient, and has a smaller autocorrelationtime than the 1 block setting. 13able 2. Summary statistics of the Gelman-Rubin diagnostic, effective sample size andautocorrelation time of the RE-SkN model estimated using Bayesian method forsimulated dataset Parameters 1 Block 4 Blocks ˆ R n eff ˆ t ˆ R n eff ˆ tµ ω β γ τ τ ξ ϕ δ δ σ u ν λ Figure 1. MCMC convergence and effficiency using 4 parameter blocks settingSecond, the estimation results of the proposed RE-SKN model based on the ML methodand the MCMC method using the 4 parameter block setting are summarized in Table 3.The optimal method is printed in bold. For each dataset, the true one step-ahead VaR14nd ES forecasts are calculated and averaged over the 1000 data sets. The average is givenin the True column. Both the ML estimators and MCMC estimators generate accurateestimates for the parameters and the VaR and ES forecasts, being very close to the truevalues. However, the MCMC method produce more favorable estimators, both in termsof unbiasedness and precision. All parameters and tail risk forecasts under the MCMCframework, except for ν parameter, have a smaller bias than those under the ML method.Regarding the estimation precision, all MCMC estimators have smaller root-mean-squareerrors (RMSE) than their respective ML counterpart.Table 3. Summary statistics of the proposed RE-SkN model estimators based on theMaximum Likelihood and MCMC Methods n=2000 True ML MCMCParameter Mean RMSE Mean RMSE µ ω -0.12 -0.1407 0.1208 -0.1342 0.0329 β γ τ -0.12 -0.1133 0.0665 -0.1210 0.0125 τ ξ -0.17 -0.0896 0.6982 -0.2002 0.1884 ϕ δ -0.09 -0.0826 0.0587 -0.0901 0.0115 δ σ u ν λ t +1 -0.0763 -0.0857 0.1169 -0.0768 0.0037 t +1 -0.0948 -0.1062 0.1459 -0.0952 0.0055 1% VaR t +1 -0.0916 -0.1047 0.1472 -0.0921 0.0049 1% ES t +1 -0.1130 -0.1291 0.1817 -0.1133 0.0073 The performance of the MCMC and ML estimators in terms of the unbiasedness andprecision of the risk forecasts is also illustrated by using the distribution of 1000 2.5%and 1% VaR and ES forecasts in Figure 2 and Figure 3, respectively. Although boththe MCMC and ML methods generate VaR and ES forecasts relatively closed to the truevalue of VaR and ES (represented by the red horizontal lines), it can be seen that theML method produces more outlying forecasts than the MCMC method. The mean valuesof both VaR and ES forecasts generated using the MCMC method are also very close totheir respective true values, compared to those of the ML method. This might be dueto the convergence issue in the ML method, which results in a higher RMSE of the riskforecasts. This suggests that the MCMC method is more preferable than the ML methodsince it produces more accurate risk forecasts.15igure 2. Histogram of 1000 2.5% and 1% VaR true and forecasts: MCMC and MLestimationFigure 3. Histogram of 1000 2.5% and 1% ES true and forecasts: MCMC and MLestimation16 . Tail risk forecasting method and evaluation The realized EGARCH models are estimated using the proposed MCMC procedures witha fixed size in-sample data combined with a rolling window approach to produce 1000one-step ahead forecasts. In a Bayesian setting, each MCMC iterates of each parameteris used to form MCMC iterates for VaR and ES forecasts at 2.5% and 1% quantile levels.The posterior mean estimates of the VaR and ES are then used as the final forecast. The original realized EGARCH model estimated using the proposed MCMC method (RE-NN) is used to calculate VaR forecasts by using the inverse cumulative distribution func-tion of the distribution (CDF), Φ − , as a benchmark to the proposed realized EGARCHmodel using the standardized Student- t distribution (RE-tN) and the standardized skewedStudent- t distribution (RE-SkN) for the return equation errors. V aR α = µ + (cid:112) h t +1 Φ − ( α ) (24)For the standardized Student- t distribution, the VaR forecasts are calculated by usingthe inverse CDF, t − ν , given by: V aR α = µ + (cid:112) h t +1 t − v ( α ) (cid:114) v − v , (25)while for the skewed Student- t distribution, the VaR estimation is conducted using theinverse CDF, skt − ν,λ as specified in Contino and Gerlach (2017) as follows: V aR α = µ + (cid:112) h t +1 skt − v,λ , (26)where skt − v,λ is skt − ( α | v, λ ) = − λb (cid:113) v − v t − v (cid:16) α − λ , v (cid:17) − ab , if α < − λ λb (cid:113) v − v t − v (cid:16) . α λ (cid:16) α − − λ (cid:17) , v (cid:17) − ab , if α ≥ − λ . Similar to the VaR, three types of ES forecasts are also calculated. For a normally dis-tributed loss, McNeil et al. (2005) formulated ES as follows: ES α = µ + (cid:112) h t +1 φ (Φ − ( α ))1 − α , (27)where φ denotes the density of the standard normal distribution and Φ − is inverse of thenormal CDF.For the standardized Student t -distribution, the ES takes the following form (McNeil17t al., 2005): ES α = µ + (cid:112) h t +1 g v ( t − v ( α ))1 − α (cid:32) v + ( t − v ( α )) v − (cid:33)(cid:114) v − v , (28)where t v is the density function, g v is the density of the standard Student- t , t − v is theinverse CDF, and v is the estimated degrees of freedom with a restriction v > t distribution, Contino and Gerlach (2017) de-rived the ES, given the standardized VaR value Θ, as follows: ES α = µ + (cid:112) h t +1 E [ (cid:15) | (cid:15) < Θ] , (29)where E [ (cid:15) | < Θ] = 1 skt v,λ (Θ) (cid:90) Θ −∞ bc (cid:32) v − (cid:32) b(cid:15) + a − λ (cid:33) (cid:33) − v +12 (cid:15)d(cid:15) = c (1 − λ ) b · skt v,λ (Θ) v − − v (cid:32) v − (cid:32) b Θ + a − λ (cid:33) (cid:33) − v − ab , (30)where skt v,λ is the CDF of the Skewed Student- t distribution: skt ( (cid:15) | v, λ ) = (1 − λ ) t v (cid:16)(cid:113) vv − (cid:16) b(cid:15) + a − λ (cid:17) , v (cid:17) if (cid:15) < − ab − λ (1 + λ ) (cid:104) t v (cid:16)(cid:113) vv − (cid:16) b(cid:15) + a − λ (cid:17) , v (cid:17) − . (cid:105) if (cid:15) ≥ − ab . Several tests are conducted to assess the accuracy of VaR and ES forecasts. As an initialassessment of VaR forecasting accuracy, the VaR violation rate (VRate) is examined.VRate is calculated as the proportion of returns exceeding the forecasted VaR level in theforecasting period, as given by: V Rate = 1 m n + m (cid:88) t = n +1 I ( r t < V aR t ) , (31)where n is the in-sample size and m is the forecasting sample size.To formally evaluate the VaR forecast accuracy, we employ three standard VaR back-test procedures that typically focus on coverage tests. These are based on the comparisonbetween the number of times that losses exceed the VaR and the expected number usingstatistical tests, such the unconditional coverage (UC) test of Kupiec (1995), the condi-tional coverage (CC) test of Christoffersen (1998), and the dynamic quantile (DQ) test ofEngle and Manganelli (2004). The UC test is employed to test whether the VaR violationrate is statistically different from the confidence level α . The CC test improves the UCtest by not only testing the joint assumption of the unconditional coverage, but also test-ing whether the probability of VaR violation is independent over time. The DQ test is a18egression-based model of hits variable on a constant, the lagged values of the hit variable,and any function of the past information. In this study, we use 4 lags of the hit variable asrecommended by Engle and Manganelli (2004). For details regarding those tests, pleaserefer to Kupiec (1995), Christoffersen (1998), and Engle and Manganelli (2004).Further, to assess whether the ES forecast is correctly specified, two recently introducedES backtests are employed. The first ES backtest utilized in this study is the bivariateExpected Shortfall Regression (ESR) backtest developed by Bayer and Dimitriadis (2018).The ESR backtest is based on a joint regression model for the quantile and the ES, wherethe returns r t is regressed on the ES forecast ˆ e t and an intercept: r t = α + β ˆ e t + u et , (32)where ES τ ( u et |F t − ) = 0.The null hypothesis of this test is H : ( α, β ) = (0 , H : ( α, β ) (cid:54) = (0 , e t = ES τ ( r t |F t − ). For details on ESR backtest, see Bayer andDimitriadis (2018). The ESR backtest is conducted based on the R code downloaded from’https://rdrr.io/github/BayerSe/esback/’.The second ES backtest employed is the test developed by Couperier and Leymarie(2019). This method is based on the theory of multi-quantile regression that exploits therelationship between the VaR and ES by using the definition of a Riemann sum, wherethe α -level ES is approximated as a finite sum of several VaR series as given by: ES α ≈ p p (cid:88) j =1 V aR t ( u j ) , (33)where the risk level u j is defined by u j = α + ( j − − αp for j = 1 , , ..., p , and p representsthe number of quantiles applied in the approximation.Couperier and Leymarie (2019) generalized the idea of Gaglianone et al. (2011) whointroduced the VaR as a regressor in a quantile regression model into the assessment ofmultiple VaRs. This is conducted by estimating a multi-quantile regression model wherethe ex-post losses { L t , t = 1 , ..., T } is regressed on the p VaR forecasts { VaR t ( u j ) , t =1 , .., T } j =1 , ,..,p L t = β ( u j ) + β ( u j ) V aR ( u j ) + (cid:15) j,t ∀ j = 1 , , ..., p, (34)where β ( u j ) and β ( u j ) denote the intercept and the slope parameters at level u j , re-spectively, and (cid:15) j,t denotes the error term at risk level u j and time t such that the u j -thconditional quantile of (cid:15) j,t satisfies Q (cid:15)j,t ( u j ; F t − ) = 0.Given the multi-quantile regression model in equation (34), the u j -th conditional quan-tile of L t is defined as: Q L t ( u j ; F t − ) = β ( u j ) + β ( u j ) V aR ( u j ) ∀ j = 1 , , ..., p. (35)19he procedure verifies whether VaR( u j ) perfectly matches Q L t ( u j ; F t − ). Under thenull hypothesis that a sequence of VaR is valid, the parameters satisfy β ( u j ) = 0 and β ( u j ) = 1, for j = 1 , , ..., p . These parameter values will produce valid tail distributionof the ES model and ES forecasts. There are four null hypotheses of this ES backtests,which assess different settings that the regression coefficients should fulfill for valid ESforecasts. These hypotheses, H ,J , H ,J , H ,I , H ,S are given by: H ,J : p (cid:88) j =1 (cid:0) β ( u j ) + β ( u j ) (cid:1) = p, (36) H ,J : p (cid:88) j =1 β ( u j ) = 0 , and p (cid:88) j =1 β ( u j ) = p, (37) H ,I : p (cid:88) j =1 β ( u j ) = 0 , (38) H ,S : p (cid:88) j =1 β ( u j ) = p, (39)where J and J denotes the joint backtests, I refers to the intercept backtest, and S S α ( Q t , r t ) = n + m (cid:88) t = n +1 ( r t − Q t )( α − I ( r t < Q t )) (40)where Q t +1 , ..., Q n + m is a series of quantile forecasts. True VaR α,t minimizes E (cid:2) S α (VaR α,t , r t ) |F t − (cid:3) . Therefore, the most accurate VaR forecasting model should minimize the scoringfunction in equation (40).The ES, however, is not an elicitable risk measure (Gneiting, 2011). Therefore, ingeneral there is no loss function minimized by the true ES. Nevertheless, Fissler and Ziegel(2016) showed that the ES is elicitable of higher order in the sense that the VaR and ESare jointly elicitable. They developed a strictly consistent function that is minimized by20he true VaR and ES: S t ( r t , V aR t , ES t ) = ( I t − α ) G ( V aR t ) − I t G ( r t )+ G ( ES t ) (cid:16) ES t − V aR t + I t α ( V aR t − r t ) (cid:17) − H ( ES t ) + a ( r t ) , (41)where I t = 1 if r t < V aR t and 0 otherwise for t = 1 , ..., T , G () is increasing, G () isstrictly increasing and strictly convex, G = H (cid:48) and lim x →−∞ G ( x ) = 0 and a( · ) is a realvalued integrable function.This study employs two special joint loss functions. The first one is based on therecommendation of Fissler and Ziegel (2016) which has been used in the study of Wanget al. (2019), while the second one is based on Taylor (2019). The first joint loss functionconsiders G ( x ) = x, G ( x ) = exp( x ) , H ( x ) = exp( x ) and a ( r t ) = 1 − log(1 − α ) that givesthe following scoring function: S t ( r t , V aR t , ES t ) = ( I t − α )( V aR t ) − I t ( r t )+ exp( ES t ) (cid:16) ES t − V aR t + I t α ( V aR t − r t ) (cid:17) − exp( ES t ) + 1 − log(1 − α ) , (42)where the loss function S = (cid:80) Tt − S t is strictly consistent and jointly minimized by thetrue VaR and ES series.The second joint loss function is a method for predicting the ES corresponding tothe VaR forecasts by using the equivalence between quantile regression and maximumlikelihood based on an asymmetric Laplace (AL) density (Taylor, 2019). In this framework,a joint model of conditional VaR and ES is estimated by maximizing the AL log-likelihood.The scoring function proposed by Taylor (2019) is based on the scoring function of Fisslerand Ziegel (2016) in equation (41) by considering G = 0 , G = − /x, H ( x ) = − log( − x )and a ( r t ) = 1 − log(1 − α ). The function is the negative of the AL log-likelihood, referredto as the AL log score, which takes the following form: S ( Q t , ES t , y t ) = − log (cid:32) α − ES t (cid:33) − ( y t − Q t )( α − I ( y t ≤ Q t )) αES t . 7. Empirical Study In this section, we present the empirical results of the proposed models, the RE-tN andRE-SkN type models, estimated using the MCMC approach using returns and 3 typesof realized measures (individual or combined) specified in Section 2. The risk forecastperformance of the models are evaluated based on various tests discussed in Section 6 andcompared with the original realized EGARCH model estimated also within a Bayesian set-ting (the RE-NN type models) . Several GARCH-family type models, all with Student- t distribution, are also included for comparison purposes, such as the GARCH, EGARCH,and GJR-GARCH models. We also include a GARCH-t-EVT model to produce condi-tional quantile and ES estimates by applying the peaks over threshold of extreme valuetheory (EVT) method to the standardized residuals of the GARCH model with Student- t distribution. We follow Taylor (2008) and set the threshold as the 10% unconditionalquantile for the lower tail. For details regarding GARCH-t-EVT models, see McNeil andFrey (2000) and Section 7.2.3 in McNeil et al. (2005). The comparison models also in-clude a filtered GARCH approach (GARCH-t-HS) used in Wang et al. (2019), where astandardized VaR and ES are estimated based on a GARCH with with Student- t distri-bution (GARCH-t). The returns ( r , .., r n ) are standardized by using by their estimatedconditional standard deviation, rt/ (cid:112) ˆ h t . The standardized VaR and ES are then multi-plied by the forecast (cid:113) ˆ h n +1 from the GARCH-t model. We also generate VaR and ESforecasts based on symmetric absolute value conditional autoregressive expectiles (CARE-SAV) model of Taylor (2008) and the realized GARCH models of Hansen et al. (2011)with a Student- t distribution for the observation equation errors for the return (RG-tN).The GARCH-type models are estimated using the Econometrics toolbox in Matlab, whilethe GARCH-t-EVT, CARE-SAV, realized GARCH, and realized EGARCH models areestimated using Matlab code developed by authors. The parameters of the models areused to generate 1000 one-step-ahead VaR and ES forecasts. This study utilizes the data of daily, including open, high, low and closing prices, as wellas high frequency data, observed at 5-minute and 1-minute intervals within trading hours.The data are downloaded from Thomson Reuters Tick History. The data are collectedfor 7 market indices: ASX200 (Australia), DAX (Germany), FTSE 100 (UK), Hang Seng(Hong Kong), NASDAQ (US), SMI (Swiss) and S&P500 (US), covering the period ofJanuary 2000 to June 2016.The daily return is calculated based on the close-to-close prices, including overnight22rice movements. The sub-sampled RV (RVSS) and sub-sampled RR (RRSS) used in thisstudy are based on the dataset used in Wang et al. (2019), in which the daily RV and RRare calculated based on the 5-minute data, and the 1-minute data are used to generatedaily RVSS and RRSS. Meanwhile, the data for the realized kernel using a Parzen weightfunction is obtained from the Oxford-Man Institute’s realized library, Library Version:0.3, which also originates from Thomson Reuters Tick History. The realized measuresmight be biased downward measure for the true volatility, since they are calculated duringthe market is open, ignoring the overnight price movement. How the realized EGARCHmodels adjust for such downward bias is discussed in section 7.2. Before discussing tail risk forecast evaluations, the estimated parameters for all forecastingsteps of the realized EGARCH models are presented. In this section, we take an exampleof the estimated parameters for the S&P500 dataset presented in Table 4 and Table 5.In general, the parameter estimates are in line with with those in Table 3 of Hansenand Huang (2016). First, γ parameter, which provides the channel for the impact of therealized measures on future volatility, ranges from 0.141 to 0.357 in the models with singlerealized measure. Furthermore, estimated γ in the models using the RVSS and RRSSare much greater (0.237 and 0.222 and on average, respectively) than those using the RK(0.046 on average), indicating that the RVSS and RRSS provide a stronger signal thanthe RK. The estimated γ in models employing the RVSS and RRSS is several times largerthan the typical value for α (about 0.05) in a conventional GARCH model, which measuresthe coefficient associated with squared returns. This suggests that the realized measuresoffers a stronger signal about future volatility than the squared return does.Second, the estimates of β are close to 1 in all cases and β − γ (cid:48) ϕ k are also relatively highalthough they are less than unity, suggesting volatility persistence. The volatility is lesspersistent in the models employing the RVSS and RRSS, compared to those employing theRK. For example, the average estimates of β − γ (cid:48) ϕ k of the RE-RV-SkN and RE-RR-SkNmodels are 0.648 and 0.622, respectively, while that of the RE-RK-SkN model is 0.829.Third, the leverage functions ( τ , τ , δ , δ ) are of the expected sign: τ and δ arenegative and τ and δ positive. The degree of asymmetry in the leverage effect is consistentacross models, where the average estimate of τ is -0.176 and τ δ and δ tend to be smaller (more negative for δ ) when using the RVSS andRRSS than the RK. The fact that δ < x t will be larger when (cid:15) t < (cid:15) t > 0, which will make h t +1 larger when (cid:15) t < γ > Model µ ω β γ γ γ τ τ ξ ξ ξ ϕ ϕ ϕ RE-RV-NN 0.0001 -0.322 0.965 0.319 -0.179 0.037 -0.934 0.967RE-RV-SkN 0.0000 -0.324 0.965 0.336 -0.186 0.039 -1.165 0.943RE-RV-tN 0.0003 -0.305 0.968 0.335 -0.183 0.038 -1.158 0.944RE-RR-NN 0.0000 -0.326 0.965 0.344 -0.177 0.039 -1.066 0.977RE-RR-SkN 0.0000 -0.327 0.965 0.357 -0.183 0.041 -1.227 0.961RE-RR-tN 0.0003 -0.313 0.967 0.356 -0.179 0.040 -1.210 0.963RE-RK-NN 0.0001 -0.264 0.972 0.141 -0.174 0.031 -0.538 1.003RE-RK-SkN 0.0001 -0.272 0.971 0.144 -0.181 0.032 -0.694 0.987RE-RK-tN 0.0003 -0.254 0.973 0.144 -0.178 0.032 -0.662 0.990RE-RV-RR-NN 0.0001 -0.355 0.962 0.185 0.125 -0.162 0.042 -0.239 -0.354 1.041 1.054RE-RV-RR-SkN 0.0001 -0.348 0.963 0.198 0.133 -0.172 0.044 -0.707 -0.817 0.991 1.004RE-RV-RR-tN 0.0003 -0.334 0.965 0.199 0.128 -0.168 0.043 -0.587 -0.702 1.005 1.018RE-RV-RK-NN 0.0002 -0.330 0.965 0.256 0.011 -0.173 0.033 -0.220 -0.150 1.043 1.043RE-RV-RK-SkN 0.0001 -0.339 0.964 0.271 0.013 -0.182 0.034 -0.539 -0.422 1.009 1.015RE-RV-RK-tN 0.0004 -0.318 0.967 0.266 0.013 -0.177 0.034 -0.425 -0.337 1.022 1.024RE-RR-RK-NN 0.0001 -0.339 0.964 0.282 0.011 -0.174 0.034 -0.408 -0.205 1.048 1.038RE-RR-RK-SkN 0.0001 -0.347 0.963 0.295 0.012 -0.182 0.036 -0.700 -0.457 1.016 1.011RE-RR-RK-tN 0.0003 -0.333 0.965 0.291 0.013 -0.178 0.036 -0.592 -0.360 1.028 1.021RE-RV-RR-RK-NN 0.0002 -0.349 0.963 0.158 0.120 0.015 -0.167 0.039 -0.244 -0.360 -0.117 1.041 1.053 1.047RE-RV-RR-RK-SkN 0.0002 -0.353 0.962 0.163 0.118 0.015 -0.169 0.039 -0.275 -0.407 -0.162 1.037 1.048 1.042RE-RV-RR-RK-tN 0.0004 -0.337 0.965 0.161 0.117 0.016 -0.165 0.040 -0.246 -0.368 -0.138 1.042 1.053 1.045 Average 0.0002 -0.323 0.966 0.237 0.222 0.046 -0.176 0.037 -0.562 -0.684 -0.353 1.007 1.019 1.022 Table 5. Estimated parameter mean for all forecasting steps for each parameter:S&P500-Part 2 Model δ δ δ δ δ δ σ u σ u σ u ν ν skw λ RE-RV-NN -0.159 0.064 0.194RE-RV-SkN -0.160 0.065 0.194 10.266 -0.167RE-RV-tN -0.156 0.064 0.194 9.594RE-RR-NN -0.153 0.065 0.175RE-RR-SkN -0.154 0.066 0.175 10.688 -0.165RE-RR-tN -0.150 0.066 0.175 10.187RE-RK-NN -0.053 0.205 0.413RE-RK-SkN -0.056 0.204 0.413 10.599 -0.166RE-RK-tN -0.042 0.205 0.413 9.750RE-RV-RR-NN -0.158 -0.151 0.063 0.066 0.194 0.176RE-RV-RR-SkN -0.159 -0.153 0.064 0.067 0.195 0.176 10.568 -0.162RE-RV-RR-tN -0.155 -0.147 0.064 0.067 0.195 0.176 9.766RE-RV-RK-NN -0.158 -0.050 0.062 0.202 0.194 0.402RE-RV-RK-SkN -0.159 -0.053 0.063 0.202 0.194 0.402 10.797 -0.161RE-RV-RK-tN -0.156 -0.041 0.063 0.203 0.194 0.401 10.093RE-RR-RK-NN -0.152 -0.054 0.064 0.204 0.176 0.402RE-RR-RK-SkN -0.153 -0.056 0.065 0.203 0.176 0.402 11.078 -0.160RE-RR-RK-tN -0.150 -0.043 0.065 0.204 0.176 0.401 10.370RE-RV-RR-RK-NN -0.157 -0.150 -0.050 0.063 0.066 0.203 0.195 0.176 0.403RE-RV-RR-RK-SkN -0.158 -0.151 -0.051 0.063 0.065 0.203 0.195 0.176 0.403 11.047 -0.157RE-RV-RR-RK-tN -0.154 -0.147 -0.038 0.063 0.066 0.204 0.195 0.175 0.403 9.762 Average -0.157 -0.151 -0.049 0.064 0.066 0.203 0.194 0.176 0.405 9.932 10.720 -0.162 ξ and ϕ would be 0 and 1, respectively. Therealized measures are measured from market open to close so they may slightly underesti-mate volatility on average. Nevertheless, the parameters in the realized EGARCH modelscan adjust for the bias. The estimated ϕ < ξ .Fifth, parameter σ u partially indicates how much error contained in the realized mea-sures, compared with log( h ), but the log-scale makes this difficult to interpret. The modelsusing the RRSS generates a smaller σ u compared to those using the RVSS and RK, mean-ing that the RRSS can potentially provide extra efficiency. This is in line with the findingsof Martens and van Dijk (2007), Christensen and Podolskij (2007), and Wang et al. (2019)that conclude that the realized range has lower mean squared errors than the realized vari-ance, which may allow models employing the realized range to produce a higher accuracyand efficiency in volatility estimation and forecasting.Sixth, the estimates of ν are around 9 to 12, indicating that the choices of the Student- t and skewed Student- t return equation errors are justified, compared with the usualGaussian distribution. Lastly, the estimated λ is negative, suggesting that the return isskewed to the left. A prudent risk management requires accurate risk forecasts. On the one hand, financialregulators are usually concerned about how many times losses exceed the VaR (the VaRviolations) and the uncovered loss size. On the other hand, risk managers must maintainthe balance between the safety goal and the profit maximization goal. An excessively highVaR will result in an inefficiency since financial institutions face large opportunity costsof capital due to reserving too much capital to buffer againts the market risks.The VaR violation rates (VRates) for each model for all market indices are presentedin Table 6 for 2.5% forecasting and Table 7 for 1% forecasting. Overall, the proposedRE-tN and RE-SkN type models generate more optimal VaR forecast series comparedto other competing models, except for the CARE-SAV model for the 2.5% forecasting.The proposed models are consistently conservative, in the sense that they the VRatesare closer to the nominal VRates of 2.5% and 1%. Other competing models tend to beanti-conservative, generating VaR forecasts much higher than 2.5% and 1%, except for theCARE-SAV and GARCH-t-HS models. The GARCH-t-EVT in particular produces toolow VaR forecasts for the 2.5% level.For the 2.5% VaR forecasting, the model with the average VRate closest to the nom-inal VRate of 2.5% is the CARE-SAV model (2.54%), followed by the RE-RR-tN model(2.59%), RE-RV-RR-tN model (2.60%), and GARCH-t-HS model (2.60%). Individually,the CARE-SAV model produces the most accurate VRates for 3 out of 7 indices, whilethe RE-RV-SkN, RE-RR-tN, and RE-RK-SkN models produce the most accurate VRatesfor 2 out of 7 indices. 25he favorable performance of the proposed RE-SkN and RE-tN type models is moredistinctive for 1% VaR forecasting level. The best models generating average VRate closestto the nominal VRate of 1% are the RE-RR-tN and RE-RK-tN models (both 1%), followedby the RE-RR-RK-tN model (0.99%) and RE-RV-RR-RK-tN model (1.04%). Individually,the model that generates the most accurate VRates is the RE-RV-RR-RK-SkN model (4out of 7 indices), followed by the RE-RR-SkN and RE-RV-RR-RK-tN models (3 out of 7indices).With respect to the number of realized measure employed, there is no clear patternwhether using multiple realized measures improves the 2.5% VaR forecast accuracy. How-ever, a clear pattern can be observed for the 1% forecasting level, where models usingmultiple realized measures produce VRate closer to 1% more frequently than those usingsingle realized measure. The result also suggests that the models utilizing the RRSS (indi-vidually or jointly with other realized measures) with the Student- t and Skewed Student- t distributions for observation equation errors tend to be more conservative and producerelatively accurate VaR violation rates. 26able 6. 2.5% VaR forecasting violation rate Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Average GARCH-t Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Average GARCH-t The result of three types of VaR backtests (the UC, CC, and DQ tests) conducted at 5%significance level is presented in Table 8. The proposed RE-SkN and RE-tN type modelshave a lower number of rejections than the competing models, as highlighted in green. Thestandard GARCH-family models have relatively high rejection rates, as expected. The RE-NN type models tend to also have a relatively high number of rejections. Meanwhile, mostof the RE-SkN and RE-tN type models are not rejected in the backtests for all indices.Specifically, the 2.5% VaR forecasts generated by the RE-RK-SkN, RE-RV-RK-tN, andRE-RV-RR-RK-SkN models are not rejected for all three types of tests. For the 1% VaRforecasts, the RE-RK-SkN and RE-RR-SkN models have no rejection for all three types oftests. This finding is consistent across all types of realized measures, either using a singleor multiple realized measures. 28able 8. VaR Backtest: UC, CC, and DQ tests Model α = 2.5% α = 1%UC CC DQ4 UC CC DQ4GARCH-t 7 4 6 4 4 5EGARCH-t 3 1 3 2 1 4GJR-t 4 3 2 2 1 3GARCH-t-HS RE-RR-TN 1 1 1 RE-RK-TN 1 1 Table 9 shows the result of the bivariate ESR backtests (Bayer and Dimitriadis, 2018)in equation (32) across all market indices at the 5% significance level. There is a smallnumber of rejection for the null hypothesis of correct ES forecasts by almost all ESR back-tests, even for the traditional GARCH-family models employing the Student- t distribution,except for the GARCH-t-EVT model. This finding is consistent with the empirical stud-ies conducted by Bayer and Dimitriadis (2018) that GARCH and GJR-GARCH modelsemploying Student- t and skewed Student- t are not rejected by the ESR backtests. FromTable 9, we can see the difference in performance across models at the 1% forecasting level,where the proposed RE-tN and RE-SkN class models could generate more accurate ES29orecasts compared to the competing models and the RE-NN class models. All proposedRE-tN and RE-SkN models have 0 ESR backtest rejection, except for the RE-RR-SkNand RE-RK-tN models. The best performing models with only 1 rejection for the 2.5% ESforecasts and 0 rejection for the 1% ES forecasts are the RE-RV-SkN, RE-RV-RK-SkN,RE-RR-RK-SkN, RE-RV-RR-RK-SkN, and RE-RV-RR-RK-tN models.Table 9. ES regression-based backtest (ESR) Model Total Rejections α =2.5% α =1% GARCH-t 1 1EGARCH-t 2 1GJR-t 1 1GARCH-t-HS 1 1GARCH-t-EVT 7 5CARE-SAV 4 3RG-RV-tN 3 1RG-RR-tN 1 1RG-RK-tN 2 1RE-RV-NN 3 1RE-RV-SkN 1 0RE-RV-tN 2 0RE-RR-NN 2 0RE-RR-SkN 1 1RE-RR-tN 2 0RE-RK-NN 2 2RE-RK-SkN 2 0RE-RK-tN 2 1RE-RV-RR-NN 2 0RE-RV-RR-SkN 2 0RE-RV-RR-tN 2 0RE-RV-RK-NN 3 0RE-RV-RK-SkN 1 0RE-RV-RK-tN 2 0RE-RR-RK-NN 2 0RE-RR-RK-SkN 1 0RE-RR-RK-tN 2 0RE-RV-RR-RK-NN 3 0RE-RV-RR-RK-SkN 1 0RE-RV-RR-RK-tN 1 0Note: Under the null hypothesis, the ES forecast are cor-rectly specified. The tests are conducted at 5% level ofsignificance. Less rejection is favored. In Table 10, the result of the ES backtests based on the multi-quantile regression (ES-MQR) approach of Couperier and Leymarie (2019) in equation (34) across models and allmarket indices is presented. The ES-MQR is conducted by using the number of quantile30 = 6 at coverage level τ = 2 . Model Hypothesis Total Rejections J J I S GARCH-t 4 4 4 4 16EGARCH-t 3 5 3 4 15GJR-t 4 3 4 5 16GARCH-t-HS 4 6 4 4 GARCH-t-EVT 7 7 7 7 CARE-SAV 2 7 RE-RK-SkN 1 6 1 3 11RE-RK-tN 2 7 1 3 13RE-RV-RR-NN 1 2 1 1 5RE-RV-RR-SkN 1 4 1 3 9RE-RV-RR-tN 2 6 1 4 13RE-RV-RK-NN RE-RV-RK-SkN RE-RR-RK-SkN RE-RV-RR-RK-SkN All four hypothesis testings of the ES-MQR backtests are conducted at the 5% signif-icance level based on bootstrap critical values. With respect to the total rejections of 4hypotheses of valid ES forecasts for 7 market indices, the ES-MQR backtests of the real-31zed EGARCH models estimated using the adaptive MCMC approach (the RE-NN,RE-tN,RE-SkN class models) are less likely to be rejected than the standard GARCH-type modelsusing the Student- t distribution. However, some of those proposed models only performslightly better than the CARE-SAV and RG-tN models, or even slightly worse than theRG-tN models. For instance, the RE-RV-tN model has a higher total number of rejectionscompared to the CARE-SAV and RG-tN models.In terms of the individual hypotheses testings, the first joint backtests that sums thethe intercept and the slope coefficients together ( H ,J ) and the intercept backtests ( H ,I )for the realized EGARCH models are less likely to reject the validity of ES forecasts. Thenon-rejection of H ,I implies that the forecasting errors are not constant across time. Theother joint backtests that sum the intercept and the slope coefficients separately ( H ,J )tend to reject the validity of ES forecasts in most models. This is due to the rejectionof the slope backtests ( H ,S ), which indicates that the errors are time-varying since theyfluctuate with respect to the VaR forecasts.While the RE-SkN type models are relatively preferable than the RE-tN type models,the ES-MQR backtests surprisingly are less likely to reject the tail risk validity producedby the RE-NN class models compared to the RE-SkN and RE-tN class models. This israther contradictory to the result of the ESR backtests which provide more supports for theRE-SkN and RE-tN type models. However, the rejection of the ES-MQR backtests for theRE-SkN and RE-tN models are mainly due to the rejection of H ,J and H ,S hypotheses.For the other two ES backtests hypothesis tests ( H ,J and H ,I ), the RE-SkN and RE-tNtype models perform as good as the RE-NN type models. Lastly, the realized EGARCHmodels with multiple realized measures have a lower ES-MQR backtest rejections thanthose with single realized measure.Figure 4. 6 quantiles VaR forecasts of the RE-RV-RK-SkN model for S&P50032 .4. Quantile loss functions In Table 11 and 12, we present the quantile loss function values at 2.5% and 1% forecastinglevels, respectively. The proposed RE-SkN and RE-tN type models in general generatesmaller quantile losses and rank better compared to other models and the RE-NN classmodels.For the 2.5% VaR forecast, the best model generating the lowest average quantileloss value is the RE-RV-RK-SkN model (63.215), followed by the RE-RV-RR-RK-SkNmodel (63.220). In terms of average rank, the RE-RV-RK-SkN model also ranks thebest, followed by the RE-RV-RK-tN model. On the contrary, the GARCH-t-EVT andGARCH-t models produce the highest average quantile loss values, which is consistentwith the VRate of these models in Table 6. On the one hand, the GARCH-t-EVT modeltends to underestimate risks with an average VRate of 1.43%, which is much lower thanthe nominal VRate of 2.5%. On the other hand, the GARCH-t model overestimate risks,producing VRate much higher than 2.5%, i.e. 3.76%.For the 1% forecast, the model that has the lowest average quantile loss value (29.541)and ranks the best (5.4 on average) is the RE-RV-RR-RK-SkN model. As has beendiscussed, the VRate of this model is also the best for 4 out of 7 market indices in thisstudy, although its average VRate across 7 market indices is not the lowest. The secondbest model is the RE-RV-RK-SkN model with the average quantile loss value of 29.643,while the RE-RV-RR-SkN model has the second best average rank (7.3). It is evidentfor both 2.5% and 1% forecasting levels that the proposed RE-SkN type models are morepreferred in VaR forecasting. 33able 11. Quantile loss function values; α = 2 . Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Mean Loss Mean Rank GARCH-t EGARCH-t 54.925 78.204 59.180 73.726 63.206 68.774 53.840 64.550 19.3GJR-t 56.707 79.343 59.779 73.114 63.810 68.708 54.000 65.066 22.3GARCH-t-HS 57.756 80.133 62.454 75.875 66.424 69.844 56.700 67.027 26.4GARCH-t-EVT CARE-SAV 55.334 78.659 59.965 RE-RV-RK-TN 53.345 78.397 57.780 72.922 61.642 RE-RR-RK-NN 53.396 80.086 58.403 72.002 61.636 67.816 52.544 63.698 12.6RE-RR-RK-SkN α = 1% Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Mean Loss Mean Rank GARCH-t CARE-SAV 25.362 RE-RV-NN 24.943 38.711 28.378 35.116 29.267 36.626 25.437 31.211 21.6RE-RV-SkN 24.876 35.926 26.976 33.231 28.452 34.398 23.735 29.656 8.0RE-RV-TN RE-RV-RR-TN 24.917 36.237 26.842 34.658 28.164 34.382 23.396 29.799 8.4RE-RV-RK-NN 25.873 38.653 28.166 34.859 29.185 36.423 24.772 31.133 22.1RE-RV-RK-SkN 25.104 36.011 26.711 33.571 28.292 34.405 23.404 RE-RV-RR-RK-TN 24.977 35.986 26.896 33.935 The values of the special case of the VaR and ES joint loss function of Fissler and Ziegel(2016) based on equation (42) are presented in Table 13 for 2.5% forecasting level and Table14 for 1% forecasting level. The proposed RE-SkN type models show a more desirableproperty compared to other models, generating lower VaR and ES joint loss values for 5out of 7 market indices for either 2.5% or 1% forecasting level. The RE-RR-SkN modelconsistently produces the lowest average of the VaR and ES joint losses for both forecastinglevels, followed by the RE-RR-RK-SkN model for the 2.5% forecasting level and the RE-RV-RR-SkN model for the 1% forecasting level. On average, the RE-RR-SkN model alsoranks first for the 1% forecasting level and ranks second for the 2.5% forecasting level (theRE-RR-RK-SkN model ranks first for the 2.5% forecasting level). As regards the type ofrealized measures, the proposed RE-SkN type models employing RRSS, both individuallyor jointly with RVSS and/or RK, are more likely to have lower VaR and ES joint losses. We35an conclude that the choice of the RRSS variable and the skewed Student- t distributionas the return equation errors improves the tail risk forecasting efficiency.Table 13. VaR and ES joint loss function values based on equation (42); α = 2 . Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Mean Loss Mean Rank GARCH-t EGARCH-t 970.7 1073.3 995.5 1042.6 1002.8 1036.1 948.2 1009.9 24.0GJR-t 980.2 1076.4 986.3 1037.2 1003.6 1030.2 951.1 1009.3 23.7GARCH-t-HS 965.2 1058.7 988.8 1045.8 1007.1 1012.1 962.9 1005.8 20.4GARCH-t-EVT CARE-SAV 974.2 RE-RR-tN 938.7 1061.1 968.6 1037.7 985.9 1003.4 927.9 989.0 9.6RE-RK-NN 957.2 1077.0 988.2 1033.1 1004.0 1016.4 931.2 1001.0 19.9RE-RK-SkN 951.6 1061.9 974.0 1029.8 992.5 1012.9 RE-RV-RR-RK-tN α = 1% Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Mean Loss Mean Rank GARCH-t 981.2 1042.2 982.0 1015.4 983.0 999.3 957.1 994.3 22.29EGARCH-t CARE-SAV 963.5 1026.1 982.0 RE-RV-SkN 940.4 1026.9 955.3 1005.3 969.3 985.5 916.8 971.4 9.86RE-RV-tN 941.5 1036.5 959.4 1004.8 967.5 993.9 911.8 973.6 12.14RE-RR-NN 940.7 1052.3 973.7 1006.3 979.4 RE-RR-tN 938.1 1025.3 954.2 1004.5 962.0 989.0 913.6 969.5 7.71RE-RK-NN 965.4 RE-RV-RR-tN 939.1 1034.1 955.6 1004.5 964.1 992.2 910.4 971.4 9.14RE-RV-RK-NN 960.4 1067.9 981.3 1015.4 995.4 1036.3 941.8 999.8 24.14RE-RV-RK-SkN 940.6 1027.6 954.7 1000.5 970.6 RE-RR-RK-tN 942.2 1033.6 RE-RV-RR-RK-tN 934.2 1033.7 956.8 1001.0 961.7 996.9 910.8 970.7 7.43Note: For individual models, green highlight indicates the favoured model, bold indicates the 2nd ranked model, bold italicsindicates the least favoured model, italics indicates the 2nd lowest ranked model, in each column. The other special case of the joint loss function of VaR and ES forecasts based on theasymmetric Laplace (AL) log score by Taylor (2019) in equation (43) across models andmarket indices is presented in Table 15 for 2.5% forecast level and Table 16 for 1% forecastlevel, respectively. The proposed RE-SkN type models also produce the smallest averageof the AL log scores (highlighted in green), i.e. for 5 out of 7 market indices, either forthe 2.5% or 1% forecast levels. Some models in this class type also have the second lowestAL log score average.The RE-RV-RR-RK-SkN model consistently generates the lowest joint loss and ranksfirst, on average, for both forecasting levels. For the 2.5% forecasting level, the secondbest models with the smallest joint loss values are the RE-RR-SkN, RE-RR-RK-SkN, andRE-RV-RK-SkN models, whereas for 1% forecasting level, the second best model is theRE-RR-SkN model.The comparison of VaR and ES joint loss functions using two different score functions,one by equation (42) and another one by equation (43), shows the superiority of theRE-SkN type models. It is also evident that the choice of the RRSS variable, either37ndividually or jointly with the other two realized measures, can improve the tail-riskforecast performance of the RE-SkN type models.Table 15. Joint Loss for VaR and ES forecasts based on AL log score in equation (43); α = 2 . Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Mean loss Mean Rank GARCH-t EGARCH-t 1.837 2.207 1.898 2.099 1.961 2.057 1.775 1.976 24.00GJR-t 1.873 2.222 1.882 2.081 1.972 2.045 1.785 1.980 24.43GARCH-t-HS 1.823 2.172 1.892 2.119 1.974 2.009 1.817 1.972 21.57GARCH-t-EVT CARE-SAV 1.831 2.170 1.892 RE-RV-NN 1.762 2.216 1.870 2.067 1.931 2.020 1.769 1.948 19.14RE-RV-SkN 1.749 2.169 1.830 2.057 1.911 RE-RV-RK-tN 1.762 2.182 1.834 2.071 1.915 1.969 1.709 1.920 12.00RE-RR-RK-NN 1.762 2.213 1.842 2.058 1.910 2.018 1.731 1.933 15.00RE-RR-RK-SkN RE-RR-RK-tN 1.755 2.182 RE-RV-RR-RK-tN α = 1% Model ASX200 DAX FTSE HK NASDAQ SMI S&P500 Mean loss Mean Rank GARCH-t EGARCH-t 2.027 2.441 2.094 2.261 2.051 2.347 1.903 2.161 GJR-t 2.029 2.439 2.025 2.231 2.063 2.335 1.900 2.146 21.71GARCH-t-HS 1.960 2.347 2.056 2.266 2.049 2.230 1.968 2.125 20.29GARCH-t-EVT CARE-SAV 1.962 2.321 2.070 RE-RV-SkN 1.900 2.303 1.965 2.202 2.032 RE-RR-RK-tN 1.918 2.320 RE-RV-RR-RK-tN 1.891 2.316 1.966 2.205 The joint loss function based on the AL log score function (equation (43)) for both 2.5%and 1% forecasting levels are further incorporated as the loss function in the calculation ofthe MCS (Hansen et al., 2011). The tests are conducted by using both R and SQ methodsat 90% and 75% confidence levels. The result of the tests are presented in Table 17. Themodels that consistently enter the MCS for both MCS methods and across market indicesare highligted in green. Overall, the proposed extensions of the realized EGARCH models(the RE-SkN and RE-tN type models), demonstrate a superior performance than thecompeting models. The RE-SkN and RE-tN type models are more likely to be includedin the R and SQ MCS methods, both at the 90% and 75% confidence levels. They arealso more favored than the original realized EGARCH models that are estimated usingthe proposed adaptive MCMC in this paper (the RE-NN type models).39able 17. 75% and 90% Model confidence set (MCS) with R and SQ methods Model α =2.5% α =1%90% MCS 75% MCS 90% MCS 75% MCSR SQ R SQ R SQ R SQ GARCH-t 4 2 1 2 6 5 4 3EGARCH-t 6 5 5 3 6 5 5 4GJR-t 5 3 4 2 6 5 6 5GARCH-t-HS 4 4 3 3 6 5 5 4GARCH-t-EVT 2 2 1 1 4 2 3 2CARE-SAV 5 4 4 2 4 4 4 2RG-RV-tN 7 3 4 3 6 4 4 2RG-RR-tN 6 5 6 3 6 5 6 4RG-RK-tN 6 2 3 2 6 6 5 3RE-RV-NN 7 4 6 3 6 5 4 2RE-RV-SkN 7 7 7 7 7 6 7 6RE-RV-tN 7 5 7 5 7 6 7 6RE-RR-NN 7 6 7 5 7 6 7 5RE-RR-SkN 7 7 7 7 6 6 6 6RE-RR-tN 7 6 6 6 6 6 5 6RE-RK-NN 7 5 6 3 6 5 4 2RE-RK-SkN 7 7 7 7 7 6 7 6RE-RK-tN 6 5 6 4 7 5 7 4RE-RV-RR-NN 7 5 6 4 6 5 5 4RE-RV-RR-SkN RE-RV-RR-tN 7 6 7 6 7 6 6 6RE-RV-RK-NN 7 4 6 3 6 5 4 3RE-RV-RK-SkN RE-RV-RK-tN RE-RR-RK-NN 7 5 6 5 6 5 5 5RE-RR-RK-SkN 7 6 7 6 6 6 6 6RE-RR-RK-tN 7 6 6 6 7 6 7 6RE-RV-RR-RK-NN 7 5 6 4 6 5 5 4RE-RV-RR-RK-SkN RE-RV-RR-RK-tN Note: For each MCS method, green indicates favored models for all marketindices. The result in Table 17 also indicates that there is a tendency that the realized EGARCHmodels employing multiple realized measures have a more favorable performance thanthose using single realized measure. In particular, there are 3 RE-SkN type models and2 RE-tN type models using multiple realized measures that are consistently included inthe MCS using both the R and SQ methods at both confidence levels, i.e. the RE-RV-RR-SkN, RE-RV-RK-SkN, RE-RV-RR-RK-SkN, RE-RV-RK-tN, and RE-RV-RR-RK-tNmodels. The RE-RV-RR-RK-SkN is also the most preferred model based on the VaR andES joint loss function values using the AL log score function (equation (43)). Meanwhile,the RE-RR-SkN model, which is the second best model based on the joint loss functionvalues using the loss functions in equation (42) and equation (43), is all included in the40CS for 7 market indices when the MCS incorporates the 2.5% VaR and ES levels.Although the RE-RR-SkN model is included in the MCS for only 6 market indices whenusing the 1% tail risk forecasting, the performance of this model overall is still highlyfavorable. 8. Conclusion In this paper, we propose a Bayesian approach to estimate the realized EGARCH model forthe purpose of VaR and ES forecasting. The original realized EGARCH model is extendedby employing Student- t and skewed Student- t distributions for the observation equationerrors (RE-tN and RE-SkN type models). The result of the simulation study demonstratesthe advantage of the MCMC estimators over maximum likelihood estimators in terms ofunbiasedness and efficiency.In the empirical study using 7 market indices, the sub-sampled realized variance, thesub-sampled realized range, as well as the realized kernel are chosen to improve the out-of-sample forecasting performance of the models. The empirical study result shows that thechoice of the return equation error distribution is found to be consequential. The proposedRE-SkN class models have a superior tail risk forecasting performance in most cases thanthe standard GARCH-type models, CARE models, realized GARCH models employingthe Student- t distribution for the observation equation errors. The proposed RE-SkN classmodels consistently produce conservative and accurate violation rates and sufficient riskcoverage, are generally less likely to be rejected in both VaR and ES backtests, produce thelowest joint loss function values, and most of the time are included in the model confidenceset.This paper also finds that both RE-SkN and and RE-tN type models are more favorablethan than the original realized EGARCH model using the Gaussian distribution estimatedunder our proposed Bayesian framework (RE-NN type models). With respect to thechoice of realized measures, the sub-sampled realized range (either employed individuallyor jointly with sub-sampled realized variance and/or realized kernel) is found to be crucialin improving the accuracy of the tail risk forecasts.We conclude that the RE-SkN and RE-tN type models incorporating sub-sampled real-ized range should be considered in tail risk forecasting. This will help financial institutionsto improve capital allocation efficiency and allow them to to have larger profit maximiza-tion opportunities while at the same time be able to maintain safe goals in accordance tothe Basel Capital Accords.This paper can be extended by: allowing other distributions for both the return andmeasurement equation errors, e.g. by using asymmetric Laplace or two-sided Weibulldistributions; incorporating a higher number of different realized measures; using otherdata frequencies and sub-sampling frequency for the realized measures; and employingother MCMC algorithms. Another area to explore is to develop a variable selection andshrinkage method for selecting realized measures in the realized EGARCH models undera Bayesian framework. 41 eferences Andersen, T. G. and T. Bollerslev (1998). Answering the skeptics: Yes, standard volatility modelsdo provide accurate forecasts. International economic review , 885–905.Andersen, T. G., T. Bollerslev, F. X. Diebold, and P. Labys (2000). Great realizations. Risk ,105–108.Andersen, T. G., T. Bollerslev, F. X. Diebold, and P. Labys (2003). Modeling and forecastingrealized volatility. Econometrica 71 (2), 579–625.Ardia, D. and L. F. Hoogerheide (2010). Efficient Bayesian estimation and combination ofGARCH-type models , Volume II, Book section 1. London: RiskBooks.Artzner, P., F. Delbaen, J. Eber, and D. Heath (1997). Thinking coherently. risk, 10. November,68 71 .Artzner, P., F. Delbaen, J. Eber, and D. Heath (1999). Coherent measures of risk. Mathematicalfinance 9 (3), 203–228.Aus´ın, M. C. and P. Galeano (2007). Bayesian estimation of the gaussian mixture garch model. Computational Statistics & Data Analysis 51 (5), 2636–2652.Bandi, F. M. and J. R. Russell (2008). Microstructure noise, realized variance, and optimalsampling. The Review of Economic Studies 75 (2), 339–369.Barndorff –Nielsen, O. E., P. R. Hansen, A. Lunde, and N. Shephard (2008). Designing realizedkernels to measure the ex post variation of equity prices in the presence of noise. Economet-rica 76 (6), 1481–1536.Barndorff –Nielsen, O. E., P. R. Hansen, A. Lunde, and N. Shephard (2009). Realized kernels inpractice: trades and quotes. The Econometrics Journal 12 (3), C1–C32.Barndorff –Nielsen, O. E. and N. Shephard (2002). Econometric analysis of realized volatility andits use in estimating stochastic volatility models. Journal of the Royal Statistical Society: SeriesB (Statistical Methodology) 64 (2), 253–280.Basel Committee on Banking Supervision (2012). Fundamental review of the trading bookconsul-tative document.Basel Committee on Banking Supervision (2016). Minimum capital requirements for market risk-consultative document.Basel Committee on Banking Supervision (2019). Minimum capital requirements for market risk-bank for international settlements.Bauwens, L., S. Laurent, and J. V. Rombouts (2006). Multivariate garch models: a survey. Journalof applied econometrics 21 (1), 79–109.Bayer, S. and T. Dimitriadis (2018). Regression based expected shortfall backtesting. arXivpreprint arXiv:1801.04112 .Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of econo-metrics 31 (3), 307–327.Bollerslev, T. (1987). A conditionally heteroskedastic time series model for speculative prices andrates of return. The review of economics and statistics , 542–547.Chen, C. W., R. Gerlach, and M. K. So (2006). Comparison of nonnested asymmetric heteroskedas-tic models. Computational Statistics & Data Analysis 51 (4), 2164–2178.Chen, C. W. and T. Watanabe (2019). Bayesian modeling and forecasting of value-at-risk viathreshold realized volatility. Applied Stochastic Models in Business and Industry 35 (3), 747–765. hen, Q., R. Gerlach, and Z. Lu (2012). Bayesian value-at-risk and expected shortfall forecastingvia the asymmetric laplace distribution. Computational Statistics & Data Analysis 56 (11),3498–3516.Chen, Q. and R. H. Gerlach (2013). The two-sided weibull distribution and forecasting financialtail risk. International Journal of Forecasting 29 (4), 527–540.Chib, S. and I. Jeliazkov (2001). Marginal likelihood from the metropolishastings output. Journalof the American Statistical Association 96 (453), 270–281.Christensen, K. and M. Podolskij (2007). Realized range-based estimation of integrated variance. Journal of Econometrics 141 (2), 323–349.Christoffersen, P. F. (1998). Evaluating interval forecasts. International economic review , 841–862.Contino, C. and R. H. Gerlach (2017). Bayesian tail-risk forecasting using realized garch. AppliedStochastic Models in Business and Industry .Couperier, O. and J. Leymarie (2019). Backtesting expected shortfall via multi-quantile regression.(halshs-01909375).Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the varianceof united kingdom inflation. Econometrica: Journal of the Econometric Society , 987–1007.Engle, R. F. and S. Manganelli (2004). Caviar: Conditional autoregressive value at risk by regres-sion quantiles. Journal of Business & Economic Statistics 22 (4), 367–381.Fang, Y. (1996). Volatility modeling and estimation of high-frequency data with gaussian noise .Thesis.Fernndez, C. and M. F. Steel (1998). On bayesian modeling of fat tails and skewness. Journal ofthe American Statistical Association 93 (441), 359–371.Fissler, T. and J. F. Ziegel (2016). Higher order elicitability and osbands principle. The Annals ofStatistics 44 (4), 1680–1707.Gaglianone, W. P., L. R. Lima, O. Linton, and D. R. Smith (2011). Evaluating value-at-risk modelsvia quantile regression. Journal of Business & Economic Statistics 29 (1), 150–160.Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin (2013). BayesianData Analysis (3rd ed.). Hoboken: CRC Press.Gerlach, R. and C. Wang (2016). Forecasting risk via realized garch, incorporating the realizedrange. Quantitative Finance 16 (4), 501–511.Giacomini, R. and I. Komunjer (2005). Evaluation and combination of conditional quantile fore-casts. Journal of Business & Economic Statistics 23 (4), 416–431.Glosten, L. R., R. Jagannathan, and D. E. Runkle (1993). On the relation between the expectedvalue and the volatility of the nominal excess return on stocks. The journal of finance 48 (5),1779–1801.Gneiting, T. (2011). Making and evaluating point forecasts. Journal of the American StatisticalAssociation 106 (494), 746–762.Gneiting, T. and A. E. Raftery (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association 102 (477), 359–378.Hansen, B. E. (1994). Autoregressive conditional density estimation. International EconomicReview , 705–730.Hansen, P. R. and Z. Huang (2016). Exponential garch modeling with realized measures of volatil-ity. Journal of Business & Economic Statistics 34 (2), 269–287.Hansen, P. R., Z. Huang, and H. H. Shek (2012). Realized garch: a joint model for returns andrealized measures of volatility. Journal of Applied Econometrics 27 (6), 877–906.Hansen, P. R. and A. Lunde (2006). Realized variance and market microstructure noise. Journalof Business & Economic Statistics 24 (2), 127–161. ansen, P. R., A. Lunde, and J. M. Nason (2011). The modal confidence set. Econometrica 79 (2),453–497.Jorion, P. (2001). Value at risk : the new benchmark for managing financial risk (2nd. ed. ed.).New York: McGraw-Hill. Includes bibliographical references (p. 521-530) and index.Kupiec, P. (1995). Techniques for verifying the accuracy of risk measurement models. The Journalof Derivatives 3 (2), 78–84.Lopez, J. A. (1999). Methods for evaluating value-at-risk estimates. Economic review 2 , 3–17.Martens, M. and D. van Dijk (2007). Measuring volatility with the realized range. Journal ofEconometrics 138 (1), 181–207.McNeil, A. J. and R. Frey (2000). Estimation of tail-related risk measures for heteroscedasticfinancial time series: an extreme value approach. Journal of empirical finance 7 (3-4), 271–300.McNeil, A. J., R. Frey, and P. Embrechts (2005). Quantitative risk management: Concepts, tech-niques and tools . Princeton University Press.Nelson, D. B. (1991). Conditional heteroskedasticity in asset returns: A new approach. Econo-metrica: Journal of the Econometric Society , 347–370.Paap, R. (2002). What are the advantages of mcmc based inference in latent variable models? Statistica Neerlandica 56 (1), 2–22.Roberts, G. O., A. Gelman, and W. R. Gilks (1997). Weak convergence and optimal scaling ofrandom walk metropolis algorithms. The Annals of Applied Probability 7 (1), 110–120.Silvapulle, M. J. and P. K. Sen (2011). Constrained statistical inference: Order, inequality, andshape constraints , Volume 912. John Wiley & Sons.Taylor, J. W. (2008). Estimating value at risk and expected shortfall using expectiles. Journal ofFinancial Econometrics 6 (2), 231–252.Taylor, J. W. (2019). Forecasting value at risk and expected shortfall using a semiparametricapproach based on the asymmetric laplace distribution. Journal of Business & Economic Statis-tics 37 (1), 121–133.Vihola, M. (2012). Robust adaptive metropolis algorithm with coerced acceptance rate. Statisticsand Computing 22 (5), 997–1008.Virbickaite, A., M. C. Ausin, and P. Galeano (2015). Bayesian inference methods for univariateand multivariate garch models: A survey. Journal of Economic Surveys 29 (1), 76–96.Wang, C., Q. Chen, and R. Gerlach (2019). Bayesian realized-garch models for financial tailrisk forecasting incorporating the two-sided weibull distribution. Quantitative Finance 19 (6),1017–1042.Watanabe, T. (2012). Quantile forecasts of financial returns using realized garch models. TheJapanese Economic Review 63 (1), 68–80.Yamai, Y. and T. Yoshiba (2005). Value-at-risk versus expected shortfall: A practical perspective. Journal of Banking & Finance 29 (4), 997–1015.Zhang, L., P. A. Mykland, and Y. At-Sahalia (2005). A tale of two time scales: Determiningintegrated volatility with noisy high-frequency data. Journal of the American Statistical Asso-ciation 100 (472), 1394–1411.(472), 1394–1411.