1. Introduction

Measuring volatility is of great importance in empirical finance to manage risk, such as portfolio allocation and to derivative asset pricing. Since volatility is not directly observable in financial markets, we must rely on some measurement techniques that estimate volatility. One of the measurement techniques, presumably the most popular, is the autoregressive conditional heteroscedasticity (ARCH) model (Engle, 1982) and its generalized version, the generalized autoregressive conditional heteroscedasticity (GARCH) model (Bollerslev, 1986). These "ARCH-type" models can successfully capture some of stylized facts of asset returns such as volatility clustering and the fat-tailed nature of return distributions. Researchers proposed numerous extensions to the ARCH or GARCH models designed to capture more features of asset returns and with the expectation of estimating more accurate volatilities. The asymmetric GARCH models such as the exponential GARCH (EGARCH) (Nelson, 1991) and GJR-GARCH (Glosten, 1993) models aim to capture the asymmetric property of the volatility process. Certain assets, such as stock prices, display asymmetric volatility. It is crucial to introduce the asymmetric volatility process in the model to obtain more accurate estimates for these assets. The rational GARCH (RGARCH) model (Takaishi, 2017) is another asymmetric model that uses a rational form in the volatility process. Most volatility processes of ARCH-type models are constructed as a polynomial of past volatilities and returns. The polynomial can be viewed as a Taylor expansion of the yet unknown true volatility process (Sentana, 1995) and the GARCH (1, 1) model corresponds to the lowest order Taylor expansion for the symmetric volatility process. The RGARCH model is based on a rational polynomial or "Padé approximants, " which could approximate the true volatility process better than Taylor polynomials do. Various fields use rational polynomials to obtain a better approximate functions. For instance, in lattice quantum chromodynamics simulations, the power of a fermion matrix is approximated by a rational function, which is used in the hybrid Monte Carlo method (Clark, 2006, 2007). This then shows that the rational form is superior to the polynomial version of the hybrid Monte Carlo method (de Forcrand, 1997; Takaishi, 2001b, 2001c, 2001a, 2002). In finance, some studies applied Padé approximants to describe return probability distributions (Nuyts, 2001; Alderweireld, 2004) and error terms of the GARCH model (Takaishi, 2012; Chen, 2013). The RGARCH model has been examined for stock returns by the goodness of fit test with the deviance information criterion (DIC) (Spiegelhalter, 2002), and some claim that the RGARCH model can be equally effective with other asymmetric GARCH models. In this study, we further examine the performance of volatility estimation with the RGARCH model. We compare the accuracy of volatility estimations by a loss function using realized volatility (RV) (Andersen, 1998; Barndorff-Nielsen, 2001) as a proxy of true volatility.

The rest of this paper is organized as follows. Section 2 introduces the RGARCH model and the other GARCH-type models we use in this study. In Section 3, we describe the Bayesian inference we use for the parameter estimation of the models and in Section 4, we present the comparative accuracy results of the volatility estimations. Finally, we conclude in Section 5.

---

2. Rational GARCH model

The RGARCH model captures volatility asymmetry in a rational form as

$ \sigma_t^2 = \frac{\omega + \alpha r_{t-1}^2 + \beta \sigma_{t-1}^2}{1+\delta r_{t-1}}, $

(1)

where $\sigma_t^2$ and $r_t$ is volatility and return at time $t$, respectively, and $\alpha, \beta, \omega$ and $\delta$ are the model parameters. $\delta$ introduces asymmetry against $r_{t-1}$ in the volatility process. If $\delta$ is positive, then the volatility will be higher for negative returns. For $\delta = 0$, the RGARCH model reduces to the well-known symmetric model, that is, the GARCH (1, 1) model. The return $r_t$ at time $t$ is defined by $r_t = \sigma_t \epsilon_t$, where $\epsilon_t \sim N(0, 1)$.

Equation (1) is not well-defined when the denominator of Equation (1) is negative. The denominator will be negative for $1 < -\delta r_{t-1}$. Although we did not experience a negative denominator in this study, it may occur for huge $r_{t-1}$ or outliers. To avoid such cases, we also propose an alternative form:

$ \sigma_t^2 = \frac{\omega + \alpha r_{t-1}^2 + \beta \sigma_{t-1}^2}{\exp(\delta r_{t-1})}. $

(2)

We call this model the "RGARCH-Exp" model. For a small $\delta r_{t-1}$, equation (2) reduces to equation (1).

To compare the accuracy of the volatility estimation, we also consider the following GARCH-type models.

● GARCH model (Bollerslev, 1986)

$ \sigma_t^2 = \omega + \alpha r_{t-1}^2 + \beta \sigma_{t-1}^2. $

(3)

● Exponential GARCH (EGARCH) model (Nelson, 1991)

$ \ln \sigma_t^2 = \omega + \beta \ln \sigma_{t-1}^2 + \theta z_{t-1} +\gamma (|z_{t-1}| -E(|z_{t-1}|)), $

(4)

where $z_{t-1} = r_{t-1}/\sigma_{t-1}$.

● GJR model (Glosten, 1993)

$ \sigma^2_t = \left\{ \begin{array}{ll} \omega + \alpha r_{t-1}^2 + \beta \sigma_{t-1}^2, & r_{t-1} \ge 0\\ \omega + (\alpha+\rho) r_{t-1}^2 + \beta \sigma_{t-1}^2, & r_{t-1} < 0. \end{array} \right. $

(5)

The GARCH model is the symmetric volatility model and the EGARCH and GJR models correspond to the asymmetric volatility models.

---

3. Bayesian inference

In this study, we employ Bayesian inference for the parameter estimation of the models. From the Bayes theorem, we obtain the posterior density $P(\theta|r)$ with $T$ observation $r = (r_1, r_2, \dots, r_T)$ given by

$ P(\theta|r) \propto L(r|\theta)\pi(\theta), $

(6)

where $L(y|\theta)$ represents the likelihood function and $\pi(\theta)$ is the prior density for $\theta$. Here, $\theta$ represents model parameters such as $\theta = (\alpha, \beta, \omega)$ for the GARCH model. We assume that the prior density is constant. The likelihood function is given by

$ L(r|\theta) = \prod\limits_{i = 1}^T (2\pi\sigma_t^2)^{1/2} \exp\left(-\frac{r_i^2}{\sigma_t^2}\right). $

(7)

Using the posterior density $P(\theta|r)$, we infer the model parameters as the expectation values by

$ \left\langle \theta \right\rangle = \frac1Z \int \theta P(\theta|r)d\theta, $

(8)

where $Z = \int P(\theta|r)d\theta$ is the normalization constant.

Generally, we cannot solve equation (8) analytically. Hence, we estimate equation (8) by the Markov Chain Monte Carlo (MCMC) method using the Metropolis-Hastings (MH) algorithm (Hastings, 1970). We generate samples from the probability distribution $P(\theta|r)$. In the MH algorithm, starting from $\theta$, we first propose a candidate $\theta^\prime$ from a certain proposal density $g(\theta^\prime|\theta)$. Then, we accept the candidate $\theta^\prime$ with a probability of

$ P_{MH}(\theta^\prime|\theta) = \min \left[ 1, \frac{P(\theta^\prime|r)g(\theta|\theta^\prime)}{P(\theta|r)g(\theta^\prime|\theta)}\right]. $

(9)

If the candidate $\theta^\prime$ is rejected, we keep $\theta$.

Following Refs.(Mitsui, 2003; Asai, 2006; Takaishi, 2009a), we use a multivariate Student t-distribution for the proposal density. The $p$-dimensional multivariate Student t-distribution is given by

$ \begin{eqnarray*} g(\theta|\theta^\prime) & = & \frac{\Gamma((\nu+p)/2)/\Gamma(\nu/2)}{\det \Sigma^{1/2} (\nu\pi)^{p/2}} \nonumber \\ & & \times \left[1+\frac{(\theta-M)^t \Sigma^{-1}(\theta-M)}{\nu}\right]^{-(\nu+p)/2}, \end{eqnarray*} $

(10)

where $\theta$ and $M$ are column vectors,

$ \theta = \left[ \begin{array}{c} \theta_1 \\ \theta_2 \\ \vdots \\ \theta_p \end{array} \right], M = \left[ \begin{array}{c} M_1 \\ M_2 \\ \vdots \\ M_p \end{array} \right], $

(11)

$ M_i = E(\theta_i), $

(12)

and $p$ corresponds to the number of model parameters. $\Sigma$ is the covariance matrix defined as

$ \frac{\nu\Sigma}{\nu-2} = E[(\theta-M)(\theta-M)^t]. $

(13)

$\nu$ is a parameter to tune the shape of the Student t-distribution. In this study, we choose $\nu = 10$ (Takaishi, 2009a). Note that since equation (10) does not depend on the previous value of $\theta^\prime$, we obtain $g(\theta|\theta^\prime) = g(\theta)$.

To utilize equation (10) in the MH algorithm, we need to specify the unknown parameters $M$ and $\Sigma$. We determine $M$ and $\Sigma$ through MCMC simulations(Takaishi, 2009a, 2009b, 2009c, 2010). First, we make a short pilot run using the Metropolis algorithm (Metropolis, 1953) and accumulate some data. We then estimate $M$ and $\Sigma$ by equations (12) and (13). Second, substituting the estimated $M$ and $\Sigma$ into eq.(10), we perform the MH simulation. After accumulating more sample data, we recalculate $M$ and $\Sigma$, and update $M$ and $\Sigma$ in equation (10).

---

4. Empirical results

We use daily closing stock price time series data from the Tokyo Stock Exchange from January 5,2004 to December 30,2015 for the following six stocks: 1 Japan Tobacco (JT) Inc., 2 Canon Inc., 3 Toray Industries Inc., 4 Panasonic Co., 5 Sumitomo Metal Industries Ltd., and 6 Nippon Steel Co. The daily return $R_t$ is defined by $100\times (ln P(t+1)-ln P(t))$, where $P(t)$ is the daily closing price at day $t$. For the RV calculations, we use 1-min high frequency return data for the same period. We evaluate the performance of the volatility estimate compared to the RV as a proxy of true volatility. To be specific, we calculate the following QLIKE loss function (Patton, 2011):

$ QLIKE(\Delta) = \frac1T\sum\limits_{t = 1}^T (\ln \sigma_t^2 -\frac{RV_{t, \Delta}}{\sigma_t^2}), $

(14)

where $RV_{t, \Delta}$ is the RV at time $t$ calculated at the sampling frequency $\Delta$. The $RV_{t, \Delta}$ is defined by the sum of the squared intraday returns,

$ RV_{t, \Delta} = \sum\limits_{i = 1}^{n(\Delta)} r_{t+i\Delta, \Delta}^2. $

(15)

where $r_{t, \Delta}$ is an intraday return sampled at the frequency of $\Delta$ at $t$ and $n(\Delta)$ is the number of intraday returns in a day. We use $RV_{t, \Delta}$ at $\Delta = 1, 2, ..., 40$. Thus, there are 40 RVs in total at various sampling frequencies. In ideal circumstances, equation (15) goes to the true volatility in the limit of $\Delta \rightarrow 0$. However, in practice, due to the microstructure noise and non-trading hours, the actual RV is biased. We modify the RV using the bias correction factor introduced by Hansen and Lunde (Hansen, 2005) given by

$ C_\Delta = \frac{\sum\limits_{i = 1}^T (R_t^D-\bar{R^D})^2}{\sum\limits_{i = 1}^T RV_{t, \Delta}}, $

(16)

where $R_t^D$ is the daily return, $\bar{R^D}$ is the average of the daily returns, and $T$ is the number of daily returns. Using $C_\Delta$, the modified RV, we obtain $RV_{t, \Delta}^\prime$ using $RV_{t, \Delta}^\prime = C_\Delta RV_{t, \Delta}$.

We estimate the model parameters with the Bayesian inference described in Section 3. We perform the Bayesian inference using the MH algorithm with a multi-dimensional Student t proposal density. First, we estimate the parameters of the proposal density with the Metropolis algorithm as a pilot run, and then switch to the MH algorithm. We recalculate the parameters of the proposal density every 1,000 Monte Carlo (MC) updates using all of the MC data. After the 5000 repetition burn-in process, we collect 30,000 data points for analysis. Figure 1 shows the acceptance of the MH algorithm every 1,000 MC updates for the RGARCH-Exp model with JT stock data. The initial acceptance is low since the parameters of the proposal density are still not estimated accurately. Then, as the MC run proceeds, the acceptance increases quickly to about 70%. We obtain similar results for other cases.

Figure 2 represents daily return time series of JT, the volatility $\sigma_t^2$ estimated by the MCMC simulation with the RGARCH-Exp model and the modified RV at $\Delta = 1min$ as representative. The estimated volatility is obtained by an average over the MC samples and its accuracy compared to the RV is examined by the loss function.

Table 1 summarizes the results for the model parameters and the QLIKE loss function. We provide the QLIKE values as an average over 40 values of $QLIKE(\Delta)$ at $\Delta = 1, 2, ..., 40$. The values in square parentheses indicate the ranking based on QLIKE. Since the ranking varies depending on the stock price, we also give an average over six stocks in square parentheses in the first column at QLIKE. The average QLIKEs for all asymmetric GARCH-type models are less than that of the GARCH model, which indicates that for the stock price time series, introducing asymmetry in the volatility process is crucial to obtaining accurate volatilities. According to the average QLIKE, the RGARCH-Exp model should rank first. However, for individual stocks, the other symmetric GARCH-type models can take the rank first. For instance, the RGARCH, EGARCH, and GJR models rank first for Toray, Canon, and, Panasonic respectively. Thus, we conclude that all asymmetric GARCH-type models examined in this study perform equally well. In Table 1, we also provide the DIC (Spiegelhalter, 2002) which measures the goodness-of-fit of the models. The values in square parentheses indicate the ranking results based on DIC. We find that the DIC results show the similar ranking to the QLIKE loss function. Namley all asymmetric GARCH-type models are superior to the GARCH model and the RGARCH-Exp model ranks first.

---

5. Conclusions

We introduce the RGARCH model as an alternative GARCH model that can capture the asymmetric property of volatility. In addition to the RGARCH model, we propose a model called the RGARCH-Exp model that is more stable for outliers. We examine the performance of the volatility estimation of the RGARCH and RGARCH-Exp models with the QLIKE loss function using realized volatility to proxy for true volatility. We compare the performance of these models with that of the EGARCH, GJR, and GARCH models. We conduct empirical analyses for stock price time series for six stocks listed on the Tokyo Stock Exchange. Based on the QLIKE loss function, we find that the RGARCH, RGARCH-Exp, and other asymmetric GARCH-type models have superior performance in terms of volatility estimation compared to the GARCH model. This indicates that it is crucial to include the asymmetric volatility process in a model for accurate volatility estimations. Since our ranking varies depending on individual stock prices, and all of the asymmetric volatility models we use in this study rank first for different stocks, we conclude that performance of the RGARCH and RGARCH-Exp models are comparable to that of the asymmetric GARCH-type models.

---

Acknowledgment

Numerical calculations in this work were conducted at the Yukawa Institute Computer Facility and the facilities of the Institute of Statistical Mathematics.

---

Conflict of interest

All authors declare no conflict of interest in this paper.

---