Efficient Bayesian estimation for GARCH-type models via Sequential Monte Carlo

TL;DR

This paper introduces an efficient Bayesian estimation method for GARCH models using Sequential Monte Carlo, enabling better uncertainty quantification, model selection, and exact inference for complex models.

Contribution

It develops novel SMC-based algorithms for GARCH parameter estimation, model selection, and likelihood estimation, addressing limitations of classical methods.

Findings

Posterior distributions are non-normal, justifying Bayesian methods.

Proposed methods outperform classical inference in long time series.

Unbiased likelihood estimator enables exact Bayesian inference for complex models.

Abstract

The advantages of sequential Monte Carlo (SMC) are exploited to develop parameter estimation and model selection methods for GARCH (Generalized AutoRegressive Conditional Heteroskedasticity) style models. It provides an alternative method for quantifying estimation uncertainty relative to classical inference. Even with long time series, it is demonstrated that the posterior distribution of model parameters are non-normal, highlighting the need for a Bayesian approach and an efficient posterior sampling method. Efficient approaches for both constructing the sequence of distributions in SMC, and leave-one-out cross-validation, for long time series data are also proposed. Finally, an unbiased estimator of the likelihood is developed for the Bad Environment-Good Environment model, a complex GARCH-type model, which permits exact Bayesian inference not previously available in the literature.

Tables

Table 1. Table 1: Comparison of unbiased Monte Carlo integration estimator and importance sampling estimator

	Number of Monte Carlo Iterations	Standard Deviation
		Monte Carlo Integration	Importance Sampling
Parameter Set 1	100	3.59	1.57
	500	1.61	1.04
	1000	1.24	0.91
	2000	0.94	0.78
	5000	0.79	0.70
Parameter Set 2	100	3.13	1.51
	500	1.48	0.87
	1000	1.13	0.80
	2000	0.90	0.64
	5000	0.74	0.51
Parameter Set 3	100	3.06	1.43
	500	1.33	0.92
	1000	1.00	0.79
	2000	0.76	0.67
	5000	0.64	0.49

Table 2. Table 3: Approximated elpd LFO subscript elpd LFO \mathrm{elpd}_{\mathrm{LFO}} for one-step ahead forecasts of simulated conditional variances from t = 201 𝑡 201 t=201 to t = 1099 𝑡 1099 t=1099

	Simulated Data from GARCH$(1,1)$	Simulated Data from GJR-GARCH$(1,1)$	Simulated Data from BEGE
GARCH$(1,1)$	6572.2	3739.8	3069.1
GJR-GARCH$(1,1)$	6352.1	6505.7	5834.7
BEGE	6447.1	6528.9	6188.1

Table 3. Table 4: Approximated elpd LFO subscript elpd LFO \mathrm{elpd}_{\mathrm{LFO}} for one-step ahead forecasts of real return from t = 201 𝑡 201 t=201 (March 1943) to t = 1099 𝑡 1099 t=1099 (January 2018)

	Real Time Series Data
GARCH$(1,1)$	1597.3
GJR-GARCH$(1,1)$	1611.9
BEGE	1613.6

Table 4. Table 5: Estimation of Log Evidence - Simulated Data

	Simulated Data from GARCH$(1,1)$
	1	2	3	4	5
GARCH$(1,1)$	1875.5	1898.6	1939.2	1947.4	1786.4
GJR-GARCH$(1,1)$	1874.4	1898.4	1938.8	1945.2	1785.7
BEGE	1870.7	1895.2	1931.0	1940.3	1776.9
	Simulated Data from GJR-GARCH$(1,1)$
	1	2	3	4	5
GARCH$(1,1)$	1889.0	2079.5	2002.3	2012.7	1929.6
GJR-GARCH$(1,1)$	1895.1	2080.9	2009.0	2018.2	1936.2
BEGE	1892.4	2080.3	2008.2	2015.1	1931.8
	Simulated Data from BEGE
	1	2	3	4	5
GARCH$(1,1)$	1945.4	1862.4	1845.4	1892.2	1826.4
GJR-GARCH$(1,1)$	1954.7	1869.0	1855.6	1899.1	1837.0
BEGE	1978.8	1881.4	1873.1	1915.7	1862.6

Table 5. Table 6: Estimation of Log Evidence - Real Data

	Real Time Series Data
	Data Annealing SMC	Likelihood Annealing SMC
GARCH$(1,1)$	1829.9	1830.1
GJR-GARCH$(1,1)$	1839.4	1838.8
BEGE	1845.4	1844.9

Equations

$$r_t=\mu +u_t,\mathrm{and}$$

$$(u_t|{\sigma }_t)\sim f_{\nu }(u_t|{\sigma }_t)=\frac{\Gamma \left[\frac{1}{2}(\nu +1)\right]}{{\pi }^{\frac{1}{2}}\Gamma \left[\frac{1}{2}\nu \right]}{\left[\left(\nu -2\right){\sigma }_t^2\right]}^{-\frac{1}{2}}{\left[1+\frac{u_t^2}{\left(\nu -2\right){\sigma }_t^2}\right]}^{-\frac{1}{2}\left(\nu +1\right)},$$

$${\sigma }_t^2={\alpha }_0+\sum _{j=1}^q{\alpha }_j{u}_{t-j}^2+\sum _{i=1}^p{\beta }_i{\sigma }_{t-i}^2,t\in \mathbb{Z},$$

$${\sigma }_t^2={\alpha }_0+{\beta }_1{\sigma }_{t-1}^2+{\alpha }_1{u}_{t-1}^2,$$

$$L(\mathbf{r}|{\mathbf{\theta }}_{\mathrm{Garch}})=$$

$$-\frac{1}{2}\left[\left(\nu +1\right)\log \left(1+\frac{u_t^2}{\left(\nu -2\right){\sigma }_t^2}\right)\right],$$

$${\sigma }_t^2={\alpha }_0+{\beta }_{\sigma }{\sigma }_{t-1}^2+\varphi u_{t-1}^2+{\varphi }^{-}{u}_{t-1}^2(I_{u_{t-1}<0}),$$

$$u_t={\sigma }_p{\omega }_{p,t}-{\sigma }_n{\omega }_{n,t},\mathrm{where}$$

$${\omega }_{p,t}\sim \tilde{\Gamma }(p_t,1),\mathrm{and}$$

$${\omega }_{n,t}\sim \tilde{\Gamma }(n_t,1),$$

$$p_t=p_0+{\rho }_p{p}_{t-1}+\frac{{\varphi }_p^{+}}{2{\sigma }_p^2}{u}_{t-1}^2{I}_{u_{t-1}\ge 0}+\frac{{\varphi }_p^{-}}{2{\sigma }_p^2}{u}_{t-1}^2(1-I_{u_{t-1}\ge 0}),\mathrm{and}$$

$$n_t=n_0+{\rho }_n{n}_{t-1}+\frac{{\varphi }_n^{+}}{2{\sigma }_n^2}{u}_{t-1}^2{I}_{u_{t-1}\ge 0}+\frac{{\varphi }_n^{-}}{2{\sigma }_n^2}{u}_{t-1}^2(1-I_{u_{t-1}\ge 0}).$$

$${\sigma }_t^2\equiv va{r}_t(r_t)={\sigma }_p^2{p}_t+{\sigma }_n^2{n}_t,$$

$$L(\mathbf{r}|{\mathbf{\theta }}_{\mathrm{BEGE}})=\sum _{t=1}^T\log f_{\mathrm{BEGE}}(r_t|r_{t-1},\dots ,r_1;{\mathbf{\theta }}_{\mathrm{BEGE}}),$$

$$f_{\mathrm{BEGE}}(u_t)$$

$$={\int }_{{\omega }_{p,t}}{f}_{{\omega }_{n,t}}({\omega }_{p,t}-u_t)f({\omega }_{p,t})\mathrm{d}{\omega }_{p,t}.$$

$${\hat{f}}_{\mathrm{BEGE}}(u_t)=\frac{1}{M}\sum _{i=1}^M{f}_{{\omega }_{n,t}}({\omega }_{p,t}^i-u_t),\mathrm{where}{\omega }_{p,t}^i\stackrel{\text{i.i.d.}}{\sim }f({\omega }_{p,t}),i=1,\dots ,M.$$

$${\mathbb{E}}_g\left[\frac{h(x)}{g(x)}\right]$$

$$\approx \frac{1}{M}\sum _{i=1}^M\frac{h(x_i)}{g(x_i)},x_i\stackrel{\text{i.i.d.}}{\sim }g(\cdot ).$$

$${\hat{f}}_{\mathrm{IS}}(u_t)=\frac{1}{M}\sum _{i=1}^M\frac{h({\omega }_{p,t}^i)}{g({\omega }_{p,t}^i)},{\omega }_{p,t}^i\stackrel{\text{i.i.d.}}{\sim }g({\omega }_{p,t}),$$

$$\pi (\mathbf{\theta }|\mathbf{y})=\frac{f(\mathbf{y}|\mathbf{\theta })\pi (\mathbf{\theta })}{Z},$$

$$\alpha (\mathbf{\theta }\to \mathbf{\theta }{}^{\ast })=\min \left(1,\frac{f(\mathbf{y}|{\mathbf{\theta }}^{\ast })\pi ({\mathbf{\theta }}^{\ast })q(\mathbf{\theta }|{\mathbf{\theta }}^{\ast })}{f(\mathbf{y}|\mathbf{\theta })\pi (\mathbf{\theta })q({\mathbf{\theta }}^{\ast }|\mathbf{\theta })}\right).$$

$${\pi }_t(\mathbf{\theta }|\mathbf{y})=\frac{f(\mathbf{y}|\mathbf{\theta }{)}^{{\gamma }_t}\pi (\mathbf{\theta })}{Z_t}\propto f(\mathbf{y}|\mathbf{\theta }{)}^{{\gamma }_t}\pi (\mathbf{\theta }),$$

$$w_t^i=W_{t-1}^i{\tilde{w}}_t^i=W_{t-1}^i\frac{{\eta }_t({\mathbf{\theta }}_{t-1}^i|\mathbf{y})}{{\eta }_{t-1}({\mathbf{\theta }}_{t-1}^i|\mathbf{y})},$$

$${\mathrm{ESS}}_t=\frac{N}{1+\mathrm{CV}[w_t(\mathbf{\theta })]},$$

$$\text{ESJD}=\alpha ({\mathbf{\theta }}_t^i\to \mathbf{\theta }{}^{\ast }){\Lambda }_t^i=\alpha ({\mathbf{\theta }}_t^i\to \mathbf{\theta }{}^{\ast })({\mathbf{\theta }}_t^i-\mathbf{\theta }{}^{\ast }{)}^T{\hat{\Sigma }}_t^{-1}({\mathbf{\theta }}_t^i-\mathbf{\theta }{}^{\ast }),$$

$${\pi }_t(\mathbf{\theta }|{\mathbf{y}}_{1:t})\propto f({\mathbf{y}}_{1:t}|\mathbf{\theta })\pi (\mathbf{\theta }),$$

$$w_t^i=W_{t-1}^i\frac{f({\mathbf{y}}_{1:t}|{\mathbf{\theta }}_{t-1}^i)}{f({\mathbf{y}}_{1:t-1}|{\mathbf{\theta }}_{t-1}^i)}=W_{t-1}^{i}f({\mathbf{y}}_t|{\mathbf{y}}_{1:t-1},{\mathbf{\theta }}_{t-1}^i),$$

$$\mathrm{P}({\hat{y}}_{t+1}|y_{1:t})={\int }_{\mathbf{\theta }}\mathrm{P}({\hat{y}}_{t+1}|\mathbf{\theta },y_{1:t}){\pi }_t(\mathbf{\theta }|y_{1:t})\mathrm{d}\mathbf{\theta }.$$

$${\mathrm{elpd}}_{\mathrm{LFO}}=\sum _{t=\tau -1}^{T-1}\log p(y_{t+1}|y_{1:t}),$$

$${\mathrm{BF}}_{m_1,m_2}=\frac{P(\mathbf{y}|m_1)}{P(\mathbf{y}|m_2)}=\frac{\int P(\mathbf{y}|{\mathbf{\theta }}_{m_1},m_1)P({\mathbf{\theta }}_{m_1}|m_1)\mathrm{d}{\mathbf{\theta }}_{m_1}}{\int P(\mathbf{y}|{\mathbf{\theta }}_{m_2},m_2)P({\mathbf{\theta }}_{m_2}|m_2)\mathrm{d}{\mathbf{\theta }}_{m_2}},$$

$$\frac{Z_t}{Z_{t-1}}={\int }_{\mathbf{\theta }}f(\mathbf{y}|\mathbf{\theta }{)}^{{\gamma }_t-{\gamma }_{t-1}}{\pi }_{t-1}(\mathbf{\theta }|\mathbf{y})\mathrm{d}\mathbf{\theta }.$$

$$\log Z$$

$$\approx \sum _{t=1}^T\log \left(\sum _{i=1}^N{w}_t^i\right),$$

Full text

Keywords: Markov chain Monte Carlo; Time series analysis; Volatility distribution; Cross-validation; Data annealing

1 Introduction

In many financial applications, volatility, defined as the conditional standard deviation of asset returns, is an important quantity. Many financial applications are based on volatility forecasts, such as risk management, asset pricing and portfolio optimization (Engle,, 2001). Therefore modelling of volatility has attracted a lot of research attention. One of the most influential models that describes the dynamics of volatility is the AutoRegressive Conditional Heteroscedasticity (ARCH) model by Engle, (1982) and the Generalized ARCH (GARCH) model by Bollerslev, (1986), which is a more parsimonious extension of ARCH. The standard GARCH model of Bollerslev, (1986) has successfully captured many characteristics of financial asset returns, and forms the basis for a wide range of expressions used to capture various empirical features of returns data. For instance, in order to capture the stylised fact of the asymmetric effect (Black,, 1976), a range of models such as the exponential GARCH (EGARCH) model (Nelson,, 1991), the Glosten-Jagannathan-Runkle GARCH model (GJR-GARCH) of Glosten et al., (1993), the absolute value GARCH (AGARCH) model of Hentschel, (1995), and the Bad Environment-Good Environment (BEGE) model (Bekaert et al.,, 2015) have been developed. Given the wide range of GARCH type models, it is crucial to develop an efficient statistical inference framework to select between these models.

Frequentist inference based on maximum likelihood estimation (MLE) is commonly employed with the GARCH class of models. A limitation of MLE is that it only produces point estimates of parameters, with the uncertainty of these parameter estimates based on asymptotic theory where parameter estimates are assumed to be normally distributed (Hamilton,, 1994). With the development of new financial instruments based on volatility, such as volatility options, modelling the full predictive density of volatility is attracting attention (e.g., Corradi et al.,, 2009, 2011; and Griffin et al.,, 2016). To produce such distributions, here we develop a Bayesian approach to statistical inference for the GARCH class of models. Under this approach, the unknown parameters are regarded as random variables, whose posterior distributions may be far from normal. In contrast to the confidence interval under classical statistics, the credible interval of Bayesian statistics is easily interpreted and gives a more direct and intuitive probability statement about uncertainty in the parameter estimates, (O’Hagan,, 2004). Furthermore, our Bayesian approach can also provide full predictive distributions of the conditional volatility required for pricing volatility instruments.

An issue with Bayesian estimation is developing an efficient method for sampling from the posterior distribution based on Monte Carlo methods (Robert and Casella,, 2004). The most commonly used approach is the Markov Chain Monte Carlo (MCMC) method (Metropolis et al.,, 1953; Hastings,, 1970), which has been utilised for Bayesian analysis of the GARCH class of models (see, for example, Kim et al.,, 1998; Vrontos et al.,, 2000; and Takaishi,, 2009). Despite the popularity of MCMC, it can be costly and difficult to tune the associated parameters when using MCMC (Roberts and Rosenthal,, 2009) with conventional MCMC algorithms becoming less effective when the target distribution is far from normally distributed. In more recent times, the sequential Monte Carlo (SMC) method for static Bayesian models (Chopin,, 2002; Del Moral et al.,, 2006) has been devised as a complementary method to MCMC. SMC traverses a set of particles through a sequence of probability distributions, which starts with one that is easy to sample from and ends with the ultimate target posterior distribution. The sequence of probability distributions can be attained either using a data annealing strategy (Chopin,, 2002) or the method of likelihood annealing (Neal,, 2001). SMC iteratively applies three main steps to traverse the population of particles, which include re-weighting, re-sampling, and moving (or mutation). SMC is parallelisable and easy to adapt to become more efficient compared with the MCMC method. It is possible to use parallelisation to save computation costs especially when the likelihood function is computationally intensive, which we demonstrate is easily achievable for the GARCH class of models. In addition, SMC is efficient at exploring irregular posterior distributions, often found with GARCH type models. The intermediate samples created by SMC embedded with data annealing strategy can be reused to perform posterior predictions. Also, the estimate of evidence, which is a by-product from SMC methods, can be used in Bayesian model selection.

In this paper, we utilise SMC methods to develop a novel fully Bayesian approach for the widely used GARCH class of models. This Bayesian framework can be used for estimation, prediction and model selection, and generating a full predictive density for conditional volatility. Our Bayesian SMC approach provides an attractive alternative to the MCMC or frequentist approaches. This paper also contributes to the SMC literature as we demonstrate that SMC data annealing is useful for leave-one-out cross-validation when dealing with time series data. Finally, we develop an unbiased likelihood estimator for the BEGE model, which permits exact Bayesian inference for this model. Bekaert et al., (2015) use an approximation to the likelihood, which we show produces biased posterior distributions. South et al., (2019) use the BEGE model with biased likelihood as an example to illustrate a new SMC approach for parameter estimation. Our paper is a comprehensive demonstration of the advantages that SMC can produce for GARCH-type models.

In Section 2, we provide an introduction to GARCH models and present the three models that are investigated in the paper. The methods for deriving an unbiased estimator for the likelihood of the BEGE model are also described. Section 3 provides a brief explanation of Bayesian statistics and the general SMC framework is given in detail. An efficient strategy to choose the sequence of intermediate probability distributions is also discussed. Section 4 provides a comparison of one-step forward predictions of volatility that are derived from the fully Bayesian approach and the classical approach, respectively. A leave-one-out cross-validation method for time series data is also described for evaluating the predictive performance of the forecasted volatility. In Section 5 a Bayesian model selection tool that is used to choose between candidate models is provided. Finally, a summary of the results is provided in Section 6 and concluding discussion is given in Section 7.

2 GARCH models

The standard GARCH model introduced by Bollerslev, (1986) successfully established a dynamic model of time-varying volatility which captures some of the important characteristics of financial asset returns (Engle,, 1982). However, the classical GARCH model does not account for the stylised feature of the asymmetric effect (Black,, 1976). The GJR-GARCH (Glosten et al.,, 1993) model is an improvement over the classical GARCH model, and Engle and Ng, (1993) suggest that the model captures the asymmetric effect in stock returns best among other popular nonlinear models such as the EGARCH (Nelson,, 1991) model. Originally, the assumption of Gaussian innovations in the standard GARCH and GJR-GARCH model fails to generate a non-Gaussian distribution for returns, which is one of the most important properties of financial time series (Francq and Zakoïan,, 2011). To solve this problem to some extent, the Student-t distribution suggested by Bollerslev, (1987) is employed here. Besides, we also consider an alternative model that accommodates conditional non-Gaussianity and was based on the GJR-GARCH model, which is called bad environment-good environment (BEGE) model (Bekaert et al.,, 2015).

2.1 Classical GARCH(1,1) model

GARCH models are observation driven, where the conditional variance is a function of not only its past values, but also past innovations to returns. To begin, define the innovations in returns, $u_{t}$ as a shock to returns:

[TABLE]

where $r_{t}$ is the time series of asset returns with fixed mean $\mu$, with $u_{t}$ having a conditional Student-t distribution and a time varying conditional variance $\sigma_{t}^{2}$, with different GARCH models defining different processes for $\sigma_{t}^{2}$. In what follows, we present three models within the GARCH family considered in this paper.

Under the classical \rm{GARCH}$$(p,q) model, the conditional variance is given by:

[TABLE]

where $\alpha_{0}>0$, $\alpha_{j}>0$, $\beta_{i}>0$, $p>0$, $q>0$, $i=1,\dots,p$, $j=1,\dots,q$, and the stationary condition requires that $\sum_{j=1}^{q}\alpha_{j}+\sum_{i=1}^{p}\beta_{i}<1$ (Bollerslev,, 1986). The wide success of the parsimonious GARCH$(1,1)$ model in extensive empirical applications makes this model the workhorse of financial applications (see, for example, Engle,, 2001; Alexander and Lazar,, 2006; and Francq and Zakoïan,, 2011). The GARCH$(1,1)$ model considered here is given by:

[TABLE]

where $\alpha_{0}$, $\alpha_{1}$, and $\beta_{1}$ are subject to the constraints discussed above in the general \rm{GARCH}$$(p,q) model. The log-likelihood function of the GARCH$(1,1)$ model is given by:

[TABLE]

where $\bm{r}$ is the observed time series of asset returns, and $\bm{\theta}_{\mathrm{Garch}}$ is a parameter set of the GARCH$(1,1)$ model.

2.2 GJR-GARCH(1,1) model

The asymmetric GJR-GARCH$(1,1)$ model is a GARCH-type model with a threshold indicator function, and has a different expression on conditional volatility $\sigma_{t}^{2}$ from the classical GARCH model, which is defined as:

[TABLE]

where $I$ is a dummy variable which takes the value of $1$ if the innovation is negative, and [math] otherwise; and $\alpha_{0}>0$, $\beta_{\sigma}>0$, $\phi^{-}>0$ and $\phi>0$ to ensure the conditional variance is positive and $\beta_{\sigma}+\phi+0.5\phi^{-}<1$ to ensure the conditional variance is stationary. In this model, the negative shocks increase the conditional variance relative to positive shocks due to the asymmetric effect in stock returns. The log-likelihood function of the GJR-GARCH$(1,1)$ model has the same expression with the one of the GARCH$(1,1)$ model provided in Equation (5).

2.3 Bad Environment-Good Environment (BEGE) model

The BEGE model is based on non-Gaussian innovations and is defined as:

[TABLE]

where $\widetilde{\Gamma}(k,\theta)$111The probability density function of $\widetilde{\Gamma}(k,\theta)$ is $f(x)=\frac{1}{\Gamma(k)\theta^{k}}(x+k\theta)^{k-1}\exp\left(-\frac{1}{\theta}(x+k\theta)\right)$ for $x>-k\theta$. is the so-called centered gamma distributions (Bekaert et al.,, 2015) with a zero mean. The innovation $u_{t}$ is a linear combination of a bad environment component $\omega_{n,t}$ with a shape parameter of $n_{t}$, and a good environment component $\omega_{p,t}$ with a shape parameter of $p_{t}$. The conditional shock distribution of $u_{t}$ is generated from two gamma-distributed shocks. The time-varying shape parameters $n_{t}$ and $p_{t}$ are given by:

[TABLE]

The overall conditional variance $\sigma_{t}^{2}$ of $u_{t}$ in the BEGE model is:

[TABLE]

which is derived from the moment generating function of the centered gamma distribution (Bekaert et al.,, 2015). As with the conditional variance, moments of higher order such as conditional skewness and excess kurtosis can be easily derived with tractable expressions. The conditional skewness is given by $\mathrm{skw_{t}}(r_{t})=2\left(\sigma_{p}^{3}p_{t}-\sigma_{n}^{3}n_{t}\right)/\sigma_{t}^{3}$, and the conditional kurtosis is given by $\mathrm{kur_{t}}(r_{t})=6\left(\sigma_{p}^{4}p_{t}+\sigma_{n}^{4}n_{t}\right)/\sigma_{t}^{4}$. These higher-order conditional moments are time-varying compared with traditional asymmetric volatility GARCH models. The log-likelihood function of the BEGE model is written as:

[TABLE]

where $\bm{\theta}_{\mathrm{BEGE}}$ is a parameter set of the BEGE model, and the conditional likelihood of observation $r_{t}$ $f_{\mathrm{BEGE}}(r_{t}|r_{t-1},\dots,r_{1};\bm{\theta}_{\mathrm{BEGE}})$ is approximated numerically according to Bekaert et al., (2015). The cumulative distribution function (CDF) of $f_{\mathrm{BEGE}}(r_{t}|r_{t-1},\dots,r_{1};\bm{\theta}_{\mathrm{BEGE}})$ at two points just above and below the observation $r_{t}$ are numerically evaluated. Then a finite difference approximation method is used to find the derivative of the CDF, which approximates the probability density function of $r_{t}$. As a result of the nature of approximation in this deterministic numerical method, the approximated likelihood of $r_{t}$ is biased. Below, we propose an unbiased estimation for the likelihood of $r_{t}$ of the BEGE model with importance sampling.

2.3.1 Unbiased estimation for the BEGE likelihood

According to Bekaert et al., (2015), the BEGE-distributed variable $u_{t}$ can be rewritten as $u_{t}=\omega_{p,t}-\omega_{n,t}$, where $\omega_{p,t}\sim\widetilde{\Gamma}(p_{t},\sigma_{p})$, $\omega_{n,t}\sim\widetilde{\Gamma}(n_{t},\sigma_{n})$, and $u_{t}=r_{t}-\mu$. Then the BEGE density can take the form as

[TABLE]

An unbiased estimator of $f_{\mathrm{BEGE}}(u_{t})$ can be constructed using Monte Carlo methods. Firstly, we sample $M$ independent samples $\omega_{p,t}^{1},\dots,\omega_{p,t}^{M}$ from $f(\omega_{p,t})$, then $f_{\mathrm{BEGE}}(u_{t})$ can be estimated as

[TABLE]

One main issue of the Monte Carlo integration method is that a large number of Monte Carlo draws, $M$, may be required to increase the accuracy of the Monte Carlo estimates, which makes it more computationally intensive. Therefore, we adopt importance sampling to reduce variance of Monte Carlo estimates and improve computational efficiency. When it is hard to sample from our target distribution $h(\cdot)$ directly, we can construct an importance distribution, $g(\cdot)$, which is easy to sample from. The following identity is the basis for importance sampling

[TABLE]

In the case of the BEGE model, the target distribution is $h(\omega_{p,t})=f_{\omega_{n,t}}(\omega_{p,t}-u_{t})f(\omega_{p,t})$, and the importance sampling estimator of $f_{\mathrm{BEGE}}(u_{t})$ is given by

[TABLE]

where a sound proposal distribution $g(\omega_{p,t}^{i})$ needs to be selected. The process of constructing an efficient importance distribution for estimating $f_{\mathrm{BEGE}}(u_{t})$ is provided in the Appendix. The general importance sampling estimator in (14) is unbiased.

Andrieu and Roberts, (2009) show that using an unbiased likelihood estimator within MCMC produces a Markov chain that converges to the idealised posterior that assumes the likelihood can be computed exactly. Chopin et al., (2013) and Duan and Fulop, (2015) show that using an unbiased likelihood estimator in an SMC algorithm for data annealing and likelihood annealing, respectively, produces an algorithm that also targets the exact posterior. We adopt a Bayesian statistical framework for parameter estimation, model selection and prediction for GARCH-type models that is described in the next section.

3 Bayesian Inference

3.1 Introduction to Bayesian Statistics

Bayes’ theorem captures learning from observed data and experience. Here, $\bm{y}$ represents the observed data and $\bm{\theta}$ represents the set of parameters of a model believed to have generated $\bm{y}$. The posterior distribution of $\bm{\theta}$ is given by:

[TABLE]

where ${f(\bm{y}|\bm{\theta})}$ is the likelihood function, ${\pi(\bm{\theta})}$ is the prior, and ${Z}=f(\bm{y})=\int_{\bm{\theta}}f(\bm{y}|\bm{\theta})\pi(\bm{\theta})\mathrm{d}\bm{\theta}$ is the normalising constant of the posterior that is independent of $\bm{\theta}$.

When our interest is in generating samples of $\bm{\theta}$ from the posterior, the normalising constant ${Z}$ can be ignored and then the posterior function of $\bm{\theta}$ can be expressed as $\pi(\bm{\theta}|\bm{y})\propto f(\bm{y}|\bm{\theta})\pi(\bm{\theta})$. Given that directly sampling from the posterior distribution $\pi(\bm{\theta}|\bm{y})$ is generally infeasible, the Markov chain Monte Carlo (MCMC) method (Metropolis et al.,, 1953; Hastings,, 1970) is the most widely used sampling method. One standard MCMC algorithm is the Metropolis-Hastings (MH) algorithm, which is straightforward to implement. Under the MH algorithm, a proposal distribution $q(\bm{\theta^{*}}|\bm{\theta})$ shall be carefully chosen, where the new sample candidate value $\bm{\theta^{*}}$ that depends on the current state $\bm{\theta}$ is proposed with probability

[TABLE]

The reader is referred to Chib and Greenberg, (1995) for the theoretical details of MH algorithms. According to Andrieu and Roberts, (2009), the exact likelihood $f(\bm{y}|\bm{\theta})$ can be replaced by its unbiased estimator $\hat{f}(\bm{y}|\bm{\theta})$, with the corresponding MCMC algorithm still converging to the posterior distribution $\pi(\bm{\theta}|\bm{y})$.

However, the popular MCMC algorithm is not free from limitations. For the MH algorithm, an optimal proposal distribution $q(\bm{\theta^{*}}|\bm{\theta})$ can lead to better Markov chains that converge more quickly and efficiently explore the parameter space. But it may be costly to tune the proposal distribution. In order to improve the efficiency of MCMC algorithms, adaptive MCMC methods have been studied in the literature (e.g., Haario et al.,, 2001, 2006; Roberts and Rosenthal,, 2009). However, the Markov property of the process may be destroyed with many adaptive strategies and optimal acceptance rates targeted by some adaptive methods are based on strong posterior distribution assumptions, which are unlikely to hold in many applications. Moreover, when the target distribution is multimodal, the conventional MCMC algorithm may become stuck in a local mode for many iterations. An alternative method named sequential Monte Carlo (SMC) was designed, which could be viewed as a complementary method to MCMC. SMC has several advantages over MCMC and a more detailed description of this method is given in the next section.

3.2 Sequential Monte Carlo

3.2.1 Sequential Monte Carlo framework

Here, we are focusing on the use of SMC for static Bayesian models (Chopin,, 2002). In SMC, a set of $N$ weighted particles, $\left\{W_{t}^{i},\bm{\theta}_{t}^{i}\right\}_{i=1}^{N}$, are traversed through a sequence of posterior distributions. Therefore, we need to build a sequence of distributions $\pi_{t}(\bm{\theta}|\bm{y})$ for $t=0,\dots,T$, which explore the ultimate target distribution gradually and where $\pi_{T}(\bm{\theta}|\bm{y})$ equals $\pi(\bm{\theta}|\bm{y})$. One method to construct the sequence of distributions is the data annealing strategy (Chopin,, 2002). Under this method the target sequence is created as $\pi_{t}(\bm{\theta}|\bm{y}_{1:t})\propto f(\bm{y}_{1:t}|\bm{\theta})\pi(\bm{\theta})$, where $\bm{y}_{1:t}$ represents the observation data up to time $t$. Another method called likelihood annealing (Neal,, 2001) can also be used to generate the sequence when the whole observed dataset is already available. The sequence formed by the likelihood annealing method is given by:

[TABLE]

where $\gamma_{t}$ is referred to as the temperature and $0=\gamma_{0}\leq\dots\leq\gamma_{t}\leq\dots\leq\gamma_{T}=1$. We use both of the methods in this study, however, the data annealing approach is preferred when data are collected sequentially. A more detailed analysis and comparison between these two sequence constructions is demonstrated in Section 3.2.2. As with MCMC, unbiased likelihood estimators in SMC and exact Bayesian inference (up to Monte Carlo error) can be produced (Chopin et al.,, 2013; Duan and Fulop,, 2015).

In order to generate a properly weighted sample for target $t$, we need to reweight the particle set at target $t-1$, $\left\{W_{t-1}^{i},\bm{\theta}_{t-1}^{i}\right\}_{i=1}^{N}$. The reweighting step can be achieved by importance sampling using:

[TABLE]

(Chopin,, 2002), for $i=1,\dots,N$ where $w_{t}^{i}$ needs to be normalised, and $\eta_{t}(\cdot)$ is the unnormalised version of $\pi_{t}(\cdot)$. The so-called incremental weight $\tilde{w}_{t}^{i}$ is calculated by dividing the unnormalised target density at $t$, $\eta_{t}(\bm{\theta}_{t-1}^{i}|\bm{y})$, by the unnormalised target density at $t-1$, $\eta_{t-1}(\bm{\theta}_{t-1}^{i}|\bm{y})$. Unfortunately, the reweighting step generally leads to reductions in the weighted samples efficiency, which can be measured by the effective sample size (ESS; Liu,, 2008). ESS refers to the number of perfect samples among all the $N$ generated samples and is computed by:

[TABLE]

where CV[$w_{t}(\bm{\theta})$] is the coefficient of variation of the unnormalised weights of target $t$. The ESS can be approximated using the normalized weights as $1/\sum_{i=1}^{N}(W_{t}^{i})^{2}$. Therefore, when the ESS becomes too small, we resample the particle values proportional to their weights. This resets the approximate ESS back to $N$.

However, the resampling strategy will result in duplication, especially for particle values with relatively large weight. A moving step is required to increase the number of unique particles. The moving step is usually the most computationally expensive part of the SMC process, and the MCMC kernels of invariant distribution $\pi_{t}$ are commonly employed in the literature. Generally, in static Bayesian models $R_{t}$ iterations of a multivariate normal random walk MCMC kernel (Chopin,, 2002) are taken, the approach adopted here is designed to diversify the particles. We propose to use the MH MCMC algorithm. It is important to note that the population of particles can be used to automatically tune the covariance matrix of the multivariate normal random walk. As the proposals of particles can be rejected by the MCMC kernel, it is then necessary to repeat the MCMC kernel $R_{t}$ times to maintain a target level of diversity.

According to Salomone et al., (2018), $R_{t}$ can be determined adaptively by choosing an optimal tuning scale for the covariance matrix of the MCMC kernel. The optimal scale is selected based on expected square jumping distance (ESJD; Sherlock and Roberts,, 2009). We randomly assign each particle one of several tuning scales and we estimate the ESJD for each particle using a single MCMC iteration. The tuning scale which resulted in the highest median ESJD value is selected as optimal. The ESJD is given by:

[TABLE]

which is the product of the MH acceptance ratio and the Mahalanobis square jumping distance (Salomone et al.,, 2018), $i=1,\dots,N$, and $t=0,\dots,T$. After obtaining the optimal scale, we keep applying an MCMC kernel to all the particles until a target number of particles’ ESJD values are higher than a threshold. We set this threshold at the median value of the Mahalanobis distances between each pair of particles before the move step. As can be seen, the total number of MCMC iterations can be adaptively attained after meeting the threshold. The SMC sampling method is summarized in Algorithm 3.1.

3.2.2 Efficient sequence construction method

Under the SMC approach, the sequence of target distributions can be defined through both likelihood annealing and data annealing. These two annealing methods provide similar estimates of posterior distributions for the parameters of the three GARCH-type models, which are shown in Section 6.1. Under the data annealing strategy, the sequence of target distributions are obtained smoothly when the observations arrive one at a time. The sequence of posterior distributions constructed by using data annealing method are given by:

[TABLE]

where $t=0,\dots,T$, and $\bm{y}_{1:t}$ represents the data available up to current time $t$. The corresponding unnormalised weight of the particle $i$ at time $t$ is thus:

[TABLE]

which needs to be normalised. It is important to note that to update the weights only the conditional likelihood of the introduced data point needs to be computed. As more returns become available, the re-weighting step in data annealing does not need to re-process previously collected data. This provides extensive computational savings, especially when the time series is long. It is worth noting that all data still need to be traversed from the beginning when moving the particles with an MCMC kernel. Therefore, it is not an online method but some benefits can be provided from the re-weighting step.

The data annealing strategy is also useful for performing posterior prediction. One way to consider if a GARCH-type model is appropriate would be to focus on the model’s ability to forecast. Our Bayesian approach provides a clear method to deduce the posterior predictive distribution conditional on currently available data. The intermediate posterior distribution $\pi_{t-1}(\bm{\theta}|\bm{y}_{1:t-1})$ from data annealing is utilised when forecasting one-step ahead posterior prediction for data at time $t$. Meanwhile, the corresponding posterior predictive density that is used to measure the prediction accuracy can be estimated. Section 4 details the procedure of estimating the posterior predictive distribution.

4 Bayesian posterior predictive distribution

In this section, the fully Bayesian and traditional classical approaches are compared in the context of generating one-step ahead forecasts of conditional variances. The distribution of the one-step ahead conditional variances together with 95% prediction intervals are offered by the fully Bayesian approach, which cannot be provided by the classical approach. We also provide a fully Bayesian approach for assessing time series models’ predictive performance, which can be achieved easily with SMC data annealing.

4.1 Bayesian posterior predictive distribution

Under the fully Bayesian approach, the posterior predictive distribution for the unknown observations can be produced. New observations $\hat{y}_{t+1}$ can be drawn from the posterior predictive distribution given the original data $y_{1:t}$ up to date $t$:

[TABLE]

The posterior predictive distribution accounts for the uncertainty associated with the parameters without making asymptotic assumptions. Here, we consider one-step ahead forecasts for the posterior predictions of the conditional variance $\hat{\sigma}^{2}_{t+1}$ given its past values $\sigma_{1:t}^{2}$ and the real data $y_{1:t}$ up to time $t$.

The posterior predictive distribution is particularly easy to estimate via SMC with data annealing. Denote the SMC weighted sample from the posterior $\pi_{t}(\bm{\theta}|y_{1:t})$ as $\left\{W_{t}^{i},\bm{\theta}_{t}^{i}\right\}_{i=1}^{N}$. Then, the posterior predictive distribution can be approximated by $\left\{W_{t}^{i},(\hat{\sigma}^{2}_{t+1})^{i}\right\}_{i=1}^{N}$ where $(\hat{\sigma}^{2}_{t+1})^{i}\sim p(\sigma^{2}_{t+1}|y_{1:t},\sigma_{1:t}^{2},\bm{\theta}_{t}^{i})$ for $i=1,\ldots,N$. From the weighted sample, a point estimate (e.g. mean or median) and predictive interval can be easily estimated. The computation associated with constructing the approximation to the predictive posterior can be easily vectorised and parallelised.

4.2 Forecast through maximum likelihood estimation (MLE) approach

We compare the Bayesian SMC approach for forecasting with MLE. In order to form one-step ahead forecasts with MLE, at time $t$ we first find the maximum likelihood estimate $\bm{\theta}_{\mathrm{MLE}}^{t}$ based on the observations up to time $t$. MLE parameter estimates are obtained by maximising each model’s log-likelihood provided in Section 2 via numerical optimisation. After obtaining the MLE parameter estimates at time $t$, a point estimate of the prediction of the conditional variance from the three GARCH-type models at time $t+1$ can be obtained by utilizing one of Equations (4), (6) or (10) in Section 2.

4.3 Forecast evaluation

In order to assess which model provides more accurate forecasts, we adopt the leave-future-out cross-validation (LFO-CV) method (Bürkner et al.,, 2019), which is a modification of the popular leave-one-out cross-validation (LOO-CV) method. LOO-CV is a well-known approach to assess out-of-sample predictive performance using in-sample fits. For an out-of-sample data set $y_{t+1},\dots,y_{T}$ and a in-sample data set $y_{1},\dots,y_{t}$, LOO-CV requires to refit the model $T-t$ times based on all the data except the forecast data point on which the predictive performance is assessed. Bürkner et al., (2019) state that LOO-CV is not suitable for assessing predictive performance when the data is sequentially ordered in time since future observations may provide some information that affects in-sample model fit, and thus the predictions may be affected if only leaving one observation out at a time. Therefore, when making one-step ahead forecasts at time $t+1$ we use LFO-CV that leaves out all the future values after time $t$ to assess predictive performance for time series models.

To measure the predictive accuracy, the expected log predictive density (elpd; Vehtari et al.,, 2017) is approximated by:

[TABLE]

(Bürkner et al.,, 2019) where $\tau$ is time point that we decide to start assessing the predictions, $T$ is the number of observations in the dataset, and the density $p(y_{t+1}|y_{1:t})$ can be approximated by evaluating an estimated density, which is constructed using kernel density estimation based on the weighted sample $\left\{W_{t}^{i},\hat{y}_{t+1}^{i}\right\}_{i=1}^{N}$, at the true observation $y_{t+1}$. The standard approach to LFO-CV would be to re-compute the posterior distribution (referred to here as a refit) as each new data point is introduced, however, this would be computationally intensive. Bürkner et al., (2019) propose an algorithm to approximate LFO-CV and reduce the number of refits reducing computation time. However, our data annealing SMC naturally re-estimates the posterior distribution in a computationally efficient manner as more data are introduced, and hence, we demonstrate it is valuable for performing LFO-CV. The weighted sample from the posterior $\left\{W_{t}^{i},\bm{\theta}_{t}^{i}\right\}_{i=1}^{N}$ based on different subsets of observations $y_{1:t}$ is a by-product of posterior simulation via data annealing SMC, which means there is no need to refit the time series model again to perform and test predictions.

A simulation study and an empirical application are performed in Section 6.3 to evaluate the predictive performance for conditional variance and real returns, respectively, and explore the relative merits of the fully Bayesian approach relative to the MLE approach. To assess predictive accuracy we provide an estimated $\mathrm{elpd}_{\mathrm{LFO}}$ for each candidate model with different datasets. We also illustrate that the predictive distribution of the conditional volatility can provide us with different levels of confidence about the magnitude of the volatility at the next time point.

5 Model choice

Practitioners are not only interested in single model analysis, but also the relative performance of competing models. Under a fully Bayesian framework, the relative performance of one model against another candidate model is assessed using the Bayes factor (Jeffreys,, 1935; Kass and Raftery,, 1995). The Bayes factor compares two different models, $m_{1}$ and $m_{2}$, with parameter vectors $\bm{\theta}_{m_{1}}$ and $\bm{\theta}_{m_{2}}$ is given by:

[TABLE]

which is just the ratio of normalising constants under each model. The normalising constant is also referred to as the marginal likelihood or evidence. The larger the evidence in favour of model $m_{1}$ relative to $m_{2}$ is, the greater the Bayes factor $\mathrm{BF}_{m_{1},m_{2}}$. One challenge of implementing Bayes factor is the intensive computation of the evidence. A review of some commonly used methods of estimating the model evidence is given in Friel and Wyse, (2012).

Fortunately, an estimate of the evidence can be obtained as a by-product from SMC methods. The estimate of evidence $Z$ can be computed as $Z=Z_{T}/Z_{0}=\prod_{t=1}^{T}Z_{t}/Z_{t-1}$ with $Z_{0}=1$ (Del Moral et al.,, 2006). Under likelihood annealing SMC, the ratio of normalizing constants $Z_{t}/Z_{t-1}$ is given by:

[TABLE]

Then the estimate of log evidence $\log{Z}$ is given by:

[TABLE]

where $w_{t}^{i}$ is the unnormalised weight from Algorithm 3.1. A similar estimator under data annealing can be easily derived. It should be noted that due to the adaptive techniques employed in the SMC methods, Beskos et al., (2016) show that the estimator of log evidence is consistent but not unbiased. The estimates of log evidence of the three GARCH models obtained through both likelihood annealing and data annealing SMC are provided in Section 6.4.

Other widely used model selection tools such as the Bayesian information criterion (BIC) (Schwarz,, 1978), which provides a crude approximation to the logarithm of Bayes factor (Kass and Raftery,, 1995), and deviance information criterion (DIC) (Spiegelhalter et al.,, 2002) have restrictive assumptions and limitations. For instance, the BIC will behave poorly when the implicit prior assumed by BIC is not a multivariate normal distribution (Kuha,, 2004), and DIC is only valid when the posterior can be approximated well with a multivariate normal distribution (Gelman et al.,, 2014). However, the fully Bayesian approach based on the Bayes factor does not rely on the restrictive approximations. A simulation study provided in Section 6.4 demonstrates the Bayesian model evidence’s ability to select the true model for the motivating GARCH-type models of interest in this paper.

6 Results

The plan of this section is as follows. Section 6.1 compares the results for the BEGE model with the biased and unbiased likelihood. In Section 6.2, the performance of the MCMC and SMC methods are compared, and the results for the posterior distributions for the GARCH-type models are presented. The one-step ahead forecasts of conditional variances by using both classical MLE and Bayesian approaches are illustrated in Section 6.3. The model selection results obtained from the fully Bayesian approach are presented in Section 6.4. Computer code for implementing our methods is available at https://github.com/DanLAct/GARCH-SMC.

In this study we use $N=10,000$ particles in SMC to generate posterior distributions. The real time series data are monthly log stock returns for the S&P 500 Composite Index from July 1926 to January 2018 from the Center of Research in Securities Prices (CRSP). The parameters in the GARCH-type models we aim to estimate respectively are $\bm{\theta}_{\mathrm{Garch}(1,1)}=(\mu^{\mathrm{G}},\alpha_{0}^{\mathrm{G}},\alpha_{1},\beta_{1},\nu^{\mathrm{G}})$, $\bm{\theta}_{\mathrm{GJR}}=(\mu^{\mathrm{GJR}},\alpha_{0}^{\mathrm{GJR}},\beta_{\sigma},\phi,\phi^{-},\nu^{\mathrm{GJR}})$, and $\bm{\theta}_{\mathrm{BEGE}}=(\mu^{\mathrm{B}},p_{0},n_{0},\rho_{p},\rho_{n},\phi_{p}^{+},\phi_{n}^{+},\phi_{p}^{-},\phi_{n}^{-},\sigma_{p},\sigma_{n})$. Some constraints on the parameters are imposed initially, so that the priors are set as below. The priors on the parameters of the $\rm{GARCH}(1,1)$ model are given by:

[TABLE]

where $\mathcal{U}(\mathord{\cdot})$ denotes a uniform probability density function, $\Gamma(a,b)$ represents a gamma distribution with shape $a$ and scale $b$, and the degrees of freedom $\nu^{\mathrm{G}}>2$. To satisfy stationarity it is necessary to enforce $\alpha_{1}+\beta_{1}\leq 0.9999$. The priors on the parameters of the GJR-GARCH$(1,1)$ model are given by:

[TABLE]

To satisfy stationarity we need $\phi+\frac{1}{2}\phi^{-}+\beta_{\sigma}\leq 0.9999$. The priors on the parameters of the BEGE model are given by:

[TABLE]

Additionally, to satisfy stationarity we require that $\rho_{p}+\frac{1}{2}\phi_{p}^{+}+\frac{1}{2}\phi_{p}^{-}\leq 0.995$ and $\rho_{n}+\frac{1}{2}\phi_{n}^{+}+\frac{1}{2}\phi_{n}^{-}\leq 0.995$. All parameters are assumed to be independent in the prior.

6.1 Comparison of Biased and Unbiased BEGE

The simulated data generated from the BEGE model are used to study the comparative performance of the biased and unbiased likelihood of the BEGE model. Figure 1 demonstrates that the approximate likelihood of $f_{\mathrm{BEGE}}(u_{t})$ can lead to biased posteriors, especially for $\phi_{p}^{-}$ and $\mu^{B}$. Even though the unbiased estimators of the likelihood produced by Monte Carlo integration and importance sampling produce similar posterior distributions, the importance sampling estimating method is more efficient.

To test the efficiency of these two estimation methods, we choose three sets of parameters and 1000 independent likelihood estimates for the two estimators with different Monte Carlo sample sizes $M$. As shown in Table 1, the Monte Carlo integration requires to draw at least 2,000 times from a density to guarantee the standard deviation is lower than 1, while the importance sampling only needs 1,000 samples to achieve a similar variance. Therefore, we adopt the importance sampling estimator when using SMC to sample from the posterior for the BEGE model below.

Both the MH MCMC and SMC methods are used to sample the posterior distributions for the BEGE model with simulated data, and the results are shown in Figure 2. The red dashed lines are the parameters used to generate the simulated data and are the initial values for the MCMC. It is evident that the estimates of the posterior distribution produced by MH MCMC are less smooth than SMC and the values used to generate simulated data are also contained by the posterior distributions. This simulation study can be viewed as a validation of the effectiveness of the proposed SMC methods.

6.2 Comparison of MCMC and SMC

In this section, the real data are used to compare the performance of the MH MCMC and SMC methods. The likelihood function for the BEGE model is computationally expensive. We significantly accelerate the likelihood calculations for the proposed parameters of all particles in one iteration of the MCMC step of the SMC algorithm by taking advantage of parallel computing. More specifically, we allocate a set of these proposed parameters to each core on a multi-core machine. Using the importance sampling estimator of the likelihood and 12 cores, the computation time for the SMC algorithm with data annealing applied to the BEGE model is 22.4 hours which runs 6.7 times faster than a completely serial SMC algorithm. Furthermore, the likelihood annealing SMC applied to the BEGE model runs 6.4 times faster than a completely serial SMC algorithm that takes 103.5 hours. Although this empirical application has focused on monthly log returns, the proposed approach is also feasible for larger datasets such as daily returns. However, in the context of daily returns and unbiased BEGE likelihood, the estimation with SMC can be computationally intensive, and the computation time may vary with different datasets.

Figure 3 shows the posterior distributions for the BEGE model explored by using MH MCMC and SMC methods. The MH MCMC requires tuning of the proposal distribution. To make implementing the MH MCMC easier, the final population of SMC particles are used to form the proposal distribution. More specifically, we use a multivariate normal random walk with a covariance estimated from the SMC particles. However, even though we give MH MCMC an advantage by avoiding the tuning step, the posteriors obtained by MCMC still perform poorly compared with the ones from SMC methods in Figure 3. For instance, the posterior distribution for parameter $\sigma_{p}$ obtained by MCMC has an obvious unexpected bump, and the trace plot in Figure 4 shows that MCMC can get stuck and the Markov chain explores the parameter space slowly. The runtime of the MH MCMC algorithm with 1 million iterations is more than 2.2 times slower than the data annealing SMC that exploits parallel computing.

The univariate posterior distributions obtained by our SMC algorithms for the parameters of the GARCH$(1,1)$ and GJR-GARCH$(1,1)$ are shown in Figures 5 and 6, respectively. As can be seen in the figures, each parameters’ posterior distribution derived by data annealing SMC is very close to that obtained under likelihood annealing SMC. Except for the parameters of the fixed mean value $\bf{\mu}=(\mu^{\mathrm{G}},\mu^{\mathrm{GJR}})$ under GARCH$(1,1)$ and GJR-GARCH$(1,1)$ model, all other posterior distributions demonstrate that they are not symmetric. Particularly, for the BEGE case in Figure 3, the posterior distributions are far from normally distributed, which contradicts the assumption of normally distributed parameter estimates under standard asymptotic inference.

6.3 Posterior Predictive Distribution

6.3.1 Simulation Study

In this section we perform a simulation study which shows the relative merits of the Bayesian approach over the classical approach for forecasting the conditional variance. Firstly, we simulate the conditional variance of stock returns $\sigma_{1:1099}^{2}$ from the GARCH$(1,1)$ model, the GJR-GARCH$(1,1)$ model and the BEGE model, respectively. The parameters used to generate the three simulated data sets are randomly selected from the final population of SMC particles when modelling real time series data with the three GARCH-type models. The initial values of conditional variance at time $t=1$ of the three GARCH-type models are $(\sigma_{1}^{2})_{\mathrm{GARCH}}=0.0023$, $(\sigma_{1}^{2})_{\mathrm{GJR-GARCH}}=0.0022$ and $(\sigma_{1}^{2})_{\mathrm{BEGE}}=\sigma_{p}^{2}p_{1}+\sigma_{n}^{2}n_{1}=0.0013$. The total number of simulated data in our study is 1099. We consider the performance of one-step ahead prediction based on simulated conditional variances $t=201$ to $t=1099$. Besides using the true model to fit the simulated data, we also perform a study of model misspecification by using misspecified models to fit the simulated data and make predictions.

Here, we present the results for the one-step ahead forecasts of conditional variances of the three GARCH-type models obtained via data annealing SMC and MLE approaches. In order to provide a clearer visual interpretation, in Figure 7 we illustrate the results for forecasted conditional volatilities, which are the square roots of the forecasted conditional variances. Figure 7 shows the 95% prediction interval bounds for the last 100 predicted volatilities $\hat{\sigma}_{1000:1099}$ simulated from the three GARCH-type models by using our Bayesian approach. The 95% prediction interval is built simply by computing the 0.025 and 0.975 quantiles of the square roots of the posterior predictive sample of the conditional variances. However, the classical approach cannot produce such prediction intervals for the conditional volatility. It is obvious that these prediction intervals capture the simulated data successfully. We take posterior medians as point estimators for the Bayesian approach. After illustrating these posterior medians together with the classical point estimators of three GARCH-type models in Figure 7, we can see that the differences between the Bayesian and classical point estimators of predictions are small.

Each models’ estimated $\mathrm{elpd}_{\mathrm{LFO}}$ for the one-step ahead forecasted conditional variance are given in Table 2. Besides checking the accuracy of Bayesian predictive distributions for conditional volatility by using the true model, we also use misspecified models to fit the simulated data and compare their performance with the true model. The model with the highest value of $\mathrm{elpd}_{\mathrm{LFO}}$ is the one that has the best predictive performance. For the data simulated from the BEGE model, the log predictive density is much higher than the other two misspecified models that provide poor predictive performance. For the data simulated from the GARCH$(1,1)$ model the true model also provides the highest $\mathrm{elpd}_{\mathrm{LFO}}$. It is worth noting that the BEGE model performs well for forecasting conditional variances from data simulated from the GJR-GARCH$(1,1)$ model. As the BEGE model allows for the type of asymmetries embedded within the GJR-GARCH$(1,1)$ model, it is to be expected that the BEGE model will perform well in this context. Therefore, even though the real volatility of the real stock returns is unobservable, this simulation study provides some confidence that the distribution of the real volatility can be adequately captured by adopting the Bayesian approach.

6.3.2 Empirical Results

With the advantages of the fully Bayesian approach demonstrated in the simulation study, we implement the methods to make forecasts for real data. Figure 8 shows the time series for real stock returns and one-step ahead forecasts for conditional variances obtained via the Bayesian approach. The Bayesian posterior median value of the conditional variance is set as the point estimate for the one-step ahead prediction. It is immediately apparent that the Bayesian approach can provide efficient predictions from three GARCH-type models, which successfully capture the features of the data such as the evident growth of volatility following the negative shock of 1987 crash.

Table 3 provides each models’ estimated $\mathrm{elpd}_{\mathrm{LFO}}$ for the one-step ahead predicted real return, where it is noted that the BEGE model has the highest predictive density indicating superior predictive performance. To examine this relative performance further, Figure 9 shows the log predictive density for each time step. For much of the sample, during periods of lower volatility, all three models exhibit, high, and similar values for the log predictive density. However, in most cases where volatility spikes and the log predictive density is lower, the GARCH$(1,1)$ and GJR-GARCH$(1,1)$ models consistently produce lower log predictive density than the BEGE model. This indicates that the superior performance of the BEGE model stems from periods of relatively high volatility.

Another promising aspect of the Bayesian approach is reflected in the provision of the predictive distribution of conditional volatility, which quantifies the uncertainty in the one-step ahead predicted volatility. However, this predictive distribution cannot be analytically derived. In Figure 10, the first column shows some distributions of the predicted volatilities from the BEGE model, and the predictive distributions of returns are shown on the right. From the predictive distribution, more information relating to the one-step forward unobserved volatility can be extracted. For instance, the probability of the conditional volatility at the next time point being within any particular range can be computed. The first plot shows a distribution of the predicted conditional volatility of the BEGE model at time $t=743$ conditional on real data $t=1$ to $t=742$. From this distribution, we can see that the conditional volatilities with values around 0.065 have a high probability of being observed. In contrast, the last two plots on the left show that there is a lot of uncertainty about the relatively high and low volatility at time $t=580$ and $t=460$, respectively. For the corresponding distribution of predicted returns on the right, it is evident to see that the return at time $t=460$ is a relatively skewed distribution.

6.4 Model Choice

6.4.1 Simulation Study

In this section we perform a simulation study to assess the performance of the SMC estimator of the Bayesian model evidence for model selection for the GARCH-type models. We first simulate five sets of data from each competing model of interest, producing 15 datasets in total. The associated parameter values used to generate simulated data are randomly selected from the SMC posterior samples of the three GARCH-type models in Section 6.1. After obtaining simulated data, we use SMC to estimate the log evidence for each candidate model given different sets of simulated data.

The results of each models’ estimation of log evidence given the simulated data are presented in Table 4. As the data generating process (DGP) is known, the performance of the evidence can be assessed by checking if the true model is correctly selected according to the estimated value of the log evidence. As can be seen from Table 4, for each simulated dataset the model with the highest estimated log evidence corresponds to the DGP. This demonstrates the effectiveness of the fully Bayesian approach for model choice and the SMC evidence estimator for the GARCH-type models and length of time series investigated here. It is worth noting that for simulated data from the GARCH$(1,1)$ model, the GJR-GARCH$(1,1)$ model has only a slightly smaller log evidence value compared to the DGP. However, when data is generated from the GJR-GARCH$(1,1)$ model, the log evidence of the GARCH$(1,1)$ model is much smaller than that of the DGP. When the GARCH$(1,1)$ model is true, there is a possibility that the GJR-GARCH$(1,1)$ model can mimic the behaviour of the GARCH$(1,1)$ model since the GARCH$(1,1)$ model is nested within the GJR-GARCH$(1,1)$ model.

6.4.2 Empirical Results

We also applied our methods to the real data. Each models’ estimated log evidence for the real time series data is provided in Table 5. There is strong evidence that the BEGE model is preferred over the other two GARCH-type models for this dataset. This is consistent with the result from Bekaert et al., (2015) and reflects the fact that the BEGE model endowed with non-Gaussian innovations outperforms the traditional GARCH model and the asymmetric GJR-GARCH model, which are equipped with Student-t innovations, when applied to equity market returns. Relative to the traditional GARCH models, the BEGE framework offers a great deal more flexibility. Not only does the BEGE model allow for volatility components driving positive and negative returns, but the conditional shapes of the two distributions can also vary in much more complex ways that are impossible under more standard GARCH models. Different volatility processes for good and bad volatility and less restrictive distributional assumptions lead to a better explanation of observed returns.

Moreover, the BEGE model can also produce the one-step ahead time-varying skewness and kurtosis. However, such dynamics for the higher conditional moments cannot be generated from the other candidate models with zero conditional skewness and constant conditional kurtosis. The time series of the median of the predicted conditional skewness of the BEGE model, together with the 95% credible interval, is illustrated in Figure 11. The time-varying skewness shows very different dynamics from the fixed zero skewness assumed in GARCH$(1,1)$ and GJR-GARCH$(1,1)$ for the prediction period from $t=201$ to $t=1099$. It demonstrates that the dynamic skewness is negative for this specific period, whereas the constant zero skewness is even above the upper bound of the credible interval most of the time. The time series of the predicted excess kurtosis for the same prediction period is shown in Figure 12. It is worth noting that the one-step ahead predicted excess kurtosis of GARCH$(1,1)$ and GJR-GARCH$(1,1)$, which is a function of the degrees of freedom of the assumed Student-t innovation, is not constant due to the changing of the parameters $\nu_{t}^{\mathrm{G}}$ and $\nu_{t}^{\mathrm{GJR}}$. Therefore, more flexibility due to the two separate ‘good environment’ component with positive skewness and the ‘bad environment’ component with negative skewness, the richer time-varying higher-order moments and conditional non-Gaussianity captured by the BEGE model make this model outperforms the alternatives evidently when fitting the real data.

7 Discussion

In this article, a novel Bayesian framework for the analysis of GARCH-type models has been developed, that allows for parameter estimation, model selection and prediction. To this end, an efficient fully Bayesian approach has been developed using SMC. To obtain exact Bayesian inference for the BEGE model, an unbiased likelihood estimator was proposed to be used in place of the biased approximation adopted in Bekaert et al., (2015). This approach, in comparison to traditional estimation inference and forecasting techniques, offers a number of important advantages. Overall, in contrast to traditional methods, the Bayesian framework offers a direct method for quantifying the uncertainty surrounding both model parameter estimates and volatility forecasts. It is shown that SMC offers a number of advantages over standard Monte-Carlo methods. For such applications on long time series data, the SMC approach proposed here offers an efficient method for generating the full predictive distribution of volatility, computation of the evidence for model comparison and undertaking leave-one-out cross-validation. Although this study has focused on GARCH-type models, the benefits of the proposed approach would translate to other time series models for alternative economic or financial problems. We have demonstrated here that SMC data annealing can be an efficient approach for comparing the predictive performance of time series models. An interesting avenue for future research would be to explore this idea more comprehensively across different time series models and to draw comparisons with the approximate approach of Bürkner et al., (2019).

Acknowledgements

The authors would like to thank the reviewers and associate editor for helpful comments to improve the paper. We would also like to thank Leah South for the involvement in vectorising the code for the BEGE model and the corresponding implementation of SMC, which provided a good basis for developing this work. The first author was supported by a scholarship from the Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS). CD was supported by an Australian Research Council Discovery Project (DP200102101). Computational resources and services used in this work were provided by the HPC and Research Support Group, Queensland University of Technology, Brisbane, Australia.

• Alexander and Lazar, (2006)

Alexander, C. and Lazar, E. (2006).

Normal mixture GARCH(1, 1): Applications to exchange rate modelling.

Journal of Applied Econometrics, 21(3):307–336.

• Andrieu and Roberts, (2009)

Andrieu, C. and Roberts, G. O. (2009).

The pseudo-marginal approach for efficient Monte Carlo computations.

The Annals of Statistics, 37(2):697–725.

• Bekaert et al., (2015)

Bekaert, G., Engstrom, E., and Ermolov, A. (2015).

Bad environments, good environments: A non-Gaussian asymmetric volatility model.

Journal of Econometrics, 186(1):258–275.

• Beskos et al., (2016)

Beskos, A., Jasra, A., Kantas, N., and Thiery, A. (2016).

On the convergence of adaptive sequential Monte Carlo methods.

The Annals of Applied Probability, 26(2):1111–1146.

• Black, (1976)

Black, F. (1976).

Studies of stock market volatility changes.

Proceedings of the American Statistical Association, the Business and Economics Section.

• Bollerslev, (1986)

Bollerslev, T. (1986).

Generalized autoregressive conditional heteroskedasticity.

Journal of Econometrics, 31(3):307–327.

• Bollerslev, (1987)

Bollerslev, T. (1987).

A conditionally heteroskedastic time series model for speculative prices and rates of return.

Review of economics and statistics, 69(3):542–547.

• Bürkner et al., (2019)

Bürkner, P.-C., Gabry, J., and Vehtari, A. (2019).

Approximate leave-future-out cross-validation for time series models.

arXiv preprint arXiv:1902.06281.

• Chib and Greenberg, (1995)

Chib, S. and Greenberg, E. (1995).

Understanding the Metropolis-Hastings algorithm.

The American Statistician, 49(4):327–335.

• Chopin, (2002)

Chopin, N. (2002).

A sequential particle filter method for static models.

Biometrika, 89(3):539–552.

• Chopin et al., (2013)

Chopin, N., Jacob, P. E., and Papaspiliopoulos, O. (2013).

SMC2: an efficient algorithm for sequential analysis of state space models.

Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(3):397–426.

• Corradi et al., (2009)

Corradi, V., Distaso, W., and Swanson, N. R. (2009).

Predictive density estimators for daily volatility based on the use of realized measures.

Journal of Econometrics, 150(2):119–138.

• Corradi et al., (2011)

Corradi, V., Distaso, W., and Swanson, N. R. (2011).

Predictive inference for integrated volatility.

Journal of the American Statistical Association, 106(496):1496–1512.

• Del Moral et al., (2006)

Del Moral, P., Doucet, A., and Jasra, A. (2006).

Sequential Monte Carlo samplers.

Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3):411–436.

• Duan and Fulop, (2015)

Duan, J.-C. and Fulop, A. (2015).

Density-tempered marginalized sequential Monte Carlo samplers.

Journal of Business & Economic Statistics, 33(2):192–202.

• Engle, (2001)

Engle, R. (2001).

GARCH 101: The use of ARCH/GARCH models in applied econometrics.

Journal of Economic Perspectives, 15(4):157–168.

• Engle, (1982)

Engle, R. F. (1982).

Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation.

Econometrica: Journal of the Econometric Society, pages 987–1007.

• Engle and Ng, (1993)

Engle, R. F. and Ng, V. K. (1993).

Measuring and Testing the Impact of News on Volatility.

The Journal of Finance, 48(5):1749–1778.

• Francq and Zakoïan, (2011)

Francq, C. and Zakoïan, J. M. (2011).

GARCH Models: Structure, Statistical Inference and Financial Applications.

John Wiley & Sons.

• Friel and Wyse, (2012)

Friel, N. and Wyse, J. (2012).

Estimating the evidence–a review.

Statistica Neerlandica, 66(3):288–308.

• Gelman et al., (2014)

Gelman, A., Hwang, J., and Vehtari, A. (2014).

Understanding predictive information criteria for Bayesian models.

Statistics and Computing, 24(6):997–1016.

• Glosten et al., (1993)

Glosten, L. R., Jagannathan, R., and Runkle, D. E. (1993).

On the relation between the expected value and the volatility of the nominal excess return on stocks.

The Journal of Finance, 48(5):1779–1801.

• Griffin et al., (2016)

Griffin, J., Liu, J., and Maheu, J. (2016).

Bayesian nonparametric estimation of ex-post variance.

Technical report, University Library of Munich, Germany.

• Haario et al., (2006)

Haario, H., Laine, M., Mira, A., and Saksman, E. (2006).

DRAM: Efficient adaptive MCMC.

Statistics and Computing, 16(4):339–354.

• Haario et al., (2001)

Haario, H., Saksman, E., and Tamminen, J. (2001).

An adaptive Metropolis algorithm.

Bernoulli, 7(2):223–242.

• Hamilton, (1994)

Hamilton, J. D. (1994).

Time Series Analysis, volume 2.

Princeton university press Princeton, NJ.

• Hastings, (1970)

Hastings, W. K. (1970).

Monte Carlo sampling methods using Markov chains and their applications.

Biometrika, 57(1):97–109.

• Hentschel, (1995)

Hentschel, L. (1995).

All in the family: Nesting symmetric and asymmetric GARCH models.

Journal of Financial Economics, 39(1):71–104.

• Jeffreys, (1935)

Jeffreys, H. (1935).

Some tests of significance, treated by the theory of probability.

In Mathematical Proceedings of the Cambridge Philosophical Society, volume 31, pages 203–222. Cambridge University Press.

• Kass and Raftery, (1995)

Kass, R. E. and Raftery, A. E. (1995).

Bayes factors.

Journal of the American Statistical Association, 90(430):773–795.

• Kim et al., (1998)

Kim, S., Shephard, N., and Chib, S. (1998).

Stochastic volatility: likelihood inference and comparison with ARCH models.

The Review of Economic Studies, 65(3):361–393.

• Kuha, (2004)

Kuha, J. (2004).

AIC and BIC: Comparisons of assumptions and performance.

Sociological Methods and Research, 33(2):188–229.

• Liu, (2008)

Liu, J. S. (2008).

Monte Carlo strategies in scientific computing.

Springer Science & Business Media.

• Metropolis et al., (1953)

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953).

Equation of state calculations by fast computing machines.

The Journal of Chemical Physics, 21(6):1087–1092.

• Neal, (2001)

Neal, R. M. (2001).

Annealed importance sampling.

Statistics and Computing, 11(2):125–139.

• Nelson, (1991)

Nelson, D. B. (1991).

Conditional Heteroskedasticity in Asset Returns: A New Approach.

Econometrica: Journal of the Econometric Society, pages 347–370.

• O’Hagan, (2004)

O’Hagan, A. (2004).

Bayesian statistics: principles and benefits.

Frontis, pages 31–45.

• Robert and Casella, (2004)

Robert, C. and Casella, G. (2004).

Monte Carlo statistical Methods.

Springer New York.

• Roberts and Rosenthal, (2009)

Roberts, G. O. and Rosenthal, J. S. (2009).

Examples of adaptive MCMC.

Journal of Computational and Graphical Statistics, 18(2):349–367.

• Salomone et al., (2018)

Salomone, R., South, L. F., Drovandi, C. C., and Kroese, D. P. (2018).

Unbiased and consistent nested sampling via sequential Monte Carlo.

arXiv preprint arXiv:1805.03924.

• Schwarz, (1978)

Schwarz, G. (1978).

Estimating the Dimension of a Model.

The Annals of Statistics, 6(2):461–464.

• Sherlock and Roberts, (2009)

Sherlock, C. and Roberts, G. (2009).

Optimal scaling of the random walk Metropolis on elliptically symmetric unimodal targets.

Bernoulli, pages 774–798.

• South et al., (2019)

South, L., Pettitt, A., and Drovandi, C. C. (2019).

Sequential Monte Carlo samplers with independent Markov chain Monte Carlo proposals.

Bayesian Analysis, 14(3):773–796.

• Spiegelhalter et al., (2002)

Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and Van Der Linde, A. (2002).

Bayesian measures of model complexity and fit.

Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4):583–639.

• Takaishi, (2009)

Takaishi, T. (2009).

An adaptive Markov chain Monte Carlo method for GARCH model.

In International Conference on Complex Sciences, pages 1424–1434. Springer.

• Vehtari et al., (2017)

Vehtari, A., Gelman, A., and Gabry, J. (2017).

Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.

Statistics and Computing, 27(5):1413–1432.

• Vrontos et al., (2000)

Vrontos, I. D., Dellaportas, P., and Politis, D. N. (2000).

Full Bayesian inference for GARCH and EGARCH models.

Journal of Business & Economic Statistics, 18(2):187–198.

Appendix: Importance sampling estimator of likelihood of the BEGE model

As described in Section 2.4.1, we need an importance distribution that will produce a low variance estimator of the likelihood with a small Monte Carlo sample size. This appendix shows that the likelihood requires to be estimated only when both $p_{t}>1$ and $n_{t}>1$, and in other cases an exact likelihood can be computed. For simplicity, we only consider here the likelihood for a single observation. Four cases are considered to compute the likelihood $f_{\mathrm{BEGE}}(u)$ with the following expression

[TABLE]

where $p\geq 1$, $n\geq 1$, $\overline{\omega}_{p}=-p\sigma_{p}$, and $\overline{\omega}_{n}=-n\sigma_{n}$. The lower limit of integration in the expression for $f_{\mathrm{BEGE}}(u)$ is $\omega_{p}$, which must be greater than $\overline{\omega}_{p}$ and $\overline{\omega}_{n}+u$. Therefore, we define the lower limit of integration as $\delta=\max\left(\overline{\omega}_{p},\overline{\omega}_{n}+u\right)$.

Firstly, define $\sigma=\frac{1}{\sigma_{p}}+\frac{1}{\sigma_{n}}$. When $p=n=1$ an exact expression for $f_{\mathrm{BEGE}}(u)$ can be found by

[TABLE]

Then, the log-likelihood of the BEGE model in this case is $\log{\left[f_{\mathrm{BEGE}}(u)\right]}=\log{k_{1}}-\log{\sigma}-\sigma\delta$. When $p=1$ and $n>1$, an exact expression for the likelihood can be found as

[TABLE]

define $k_{2}=\frac{1}{\sigma_{p}}\frac{1}{\Gamma(n)\sigma_{n}^{n}}\exp\left(\frac{\overline{\omega}_{p}}{\sigma_{p}}\right)$, and $c_{1}=k_{2}\exp\left(-\frac{1}{\sigma_{p}}\left(u+\overline{\omega}_{n}\right)\right)$, then the likelihood can be written as

[TABLE]

We set $x_{p}=\omega_{p}-u-\overline{\omega}_{n}=\omega_{p}-(u+\overline{\omega}_{n})$, which gives $\omega_{p}=x_{p}+u+\overline{\omega}_{n}$ and $\mathrm{d}\omega_{p}=\mathrm{d}x_{p}$. The likelihood then can be written as

[TABLE]

where $F$ is the CDF of a gamma distribution with parameters $n$ and $\sigma$ over the interval of $\left(0,\overline{\omega}_{p}-(u+\overline{\omega}_{n})\right)$. Similarly, in the case of $n=1$ and $p>1$ an exact expression for the likelihood can be found as

[TABLE]

define $k_{3}=\frac{1}{\sigma_{n}}\frac{1}{\Gamma(p)\sigma_{p}^{p}}\exp\left(\frac{u+\overline{\omega}_{n}}{\sigma_{n}}\right)$ and $c_{2}=k_{3}\exp\left(-\frac{\overline{\omega}_{p}}{\sigma_{n}}\right)$, then the likelihood can be written as

[TABLE]

We set $x_{p}=\omega_{p}-\overline{\omega}_{p}$, which gives $\omega_{p}=x_{p}+\overline{\omega}_{p}$ and $\mathrm{d}\omega_{p}=\mathrm{d}x_{p}$. The likelihood then can be written as

[TABLE]

where $F$ is the CDF of a gamma distribution with parameters $p$ and $\sigma$ over the interval of $\left(0,u+\overline{\omega}_{n}-\overline{\omega}_{p}\right)$.

In the case of $p>1$ and $n>1$, we use importance sampling method to estimate likelihood of the BEGE model. Here we attempt to find a parametric importance distribution $g(\omega_{p})$ that matched closely the target distribution $h(\omega_{p})$. Firstly, we determine the mode and curvature of $h(\omega_{p})$. The expression for $h(\omega_{p})$ is given by

[TABLE]

We set $m(\omega_{p})=\log\left[h(\omega_{p})\right]$, which is given by

[TABLE]

The first and second derivatives of $m(\omega_{p})$ are given by

[TABLE]

The mode $\hat{\omega}_{p}$ of $m(\omega_{p})$ can be found by setting $m^{\prime}(\omega_{p})=0$, then we get

[TABLE]

By setting $a=\sigma$, $b=-\sigma(u+\overline{\omega}_{n}+\overline{\omega}_{p})+(2-n-p)$, and $c=\sigma\overline{\omega}_{p}(u+\overline{\omega}_{n})+(p-1)(u+\overline{\omega}_{n})+(n-1)\overline{\omega}_{p}$, then we can write $a\omega_{p}^{2}+b\omega_{p}+c=0$ and the mode $\hat{\omega}_{p}=\frac{-b+\sqrt{b^{2}–4ac}}{2a}$ can be easily derived. The variance of $m(\omega_{p})$ is estimated by using the second-derivative of $m(\omega_{p})$ around the mode $\hat{\omega}_{p}$, which is given by

[TABLE]

We propose a gamma distribution, $X\sim\Gamma(i,j)$ for $x>0$, with mode $\hat{\omega}_{p}$ and variance $\hat{V}$ to be the importance distribution. The mode of this gamma distribution is $(i-1)j$ and the variance is $ij^{2}$. As $\omega_{p}>\delta$, we need to make the importance distribution greater than $\delta$ and a shift gamma distribution, $Y=X+\delta$ for $y>\delta$, is selected. As $x=y-\delta$, we can show that $\left|\frac{\mathrm{d}x}{\mathrm{d}y}\right|=1$. Then an expression for the density of this importance distribution $g(y)$ is given by

[TABLE]

where the mode of this distribution is $(i-1)j+\delta$ and the variance is $ij^{2}$. The value of parameter $i$ and $j$ can be derived by solving $(i-1)j+\delta=\hat{\omega}_{p}$ and $ij^{2}=\hat{V}$, and the results are given by

[TABLE]

After obtaining the importance distribution $g(\cdot)$, we can unbiasedly estimate $f_{\mathrm{BEGE}}(u)$ by following Equation 15. It is worth noting that this importance distribution can be used only when the mode $\hat{\omega}_{p}$ is greater than $\delta$, which is the lower limit of integration in the expression for $f_{\mathrm{BEGE}}(u)$. For this case, it can be proved that $\hat{\omega}_{p}$ is always greater than $\delta$, and the details of proof are shown below.

[TABLE]