Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
64,882 characters · 17 sections · 60 citation commands
Stochastic volatility model with range-based correction and leverage
\allowdisplaybreaks \baselineskip=20pt
\thispagestyle{empty}
This study presents contemporaneous modeling of asset return and price range within the framework of stochastic volatility with leverage. A new representation of the probability density function for the price range is provided, and its accurate sampling algorithm is developed. A Bayesian estimation using Markov chain Monte Carlo (MCMC) method is provided for the model parameters and unobserved variables. MCMC samples can be generated rigorously, despite the estimation procedure requiring sampling from a density function with the sum of an infinite series. The empirical results obtained using data from the U.S. market indices are consistent with the stylized facts in the financial market, such as the existence of the leverage effect. In addition, to explore the model's predictive ability, a model comparison based on the volatility forecast performance is conducted.
{\it Key words}: alternating series, leverage, Markov chain Monte Carlo, stochastic volatility, range.
\setcounter{page}{1}
Estimating the variance or volatility of an asset return observed in the financial market has been an attractive research topic for many authors. Stochastic volatility (SV) model and generalized autoregressive conditional heteroskedasticity (GARCH) model are two major statistical models for asset return through describing the dynamics of the variance using state space representation (Bauwens_etal(2012)).
Many scholars have studied a model-free approach for estimating the volatility of asset return using price range since the 1980s. Price range is defined as the logarithmic difference between the highest and lowest asset prices of daily transactions, and Feller(1951) derived the distribution as the sum of an infinite series. Parkinson(1980) derived the $p$-th moment of the range and provided an estimator of the variance, although GarmanKlass(1980) pointed out that the ranges observed in the financial market may have (downward) bias due to price discretization noise and non-trading hours. GarmanKlass(1980), RogersSatchell(1991), Kunitomo(1992), and YangZhang(2000) explored such bias and proposed volatility estimators based on a range (also see Molnar(2012)). On the other hand, some authors took a model-based approach for volatility estimation using the price range. For example, Gallant_etal(1999), Alizadeh_etal(2002), and BrandtJones(2005) investigated SV models using a daily range of the log asset process. Horst_etal(2011) presented a joint model for daily asset return, the highest and lowest log prices within the context of SV models. Chou(2005), BrandtJones(2006), and Chen_etal(2008) presented GARCH-related models that make use of the range. I notice that these articles approximate the likelihood function as an easy density or sum of a finite series since the distribution function of the range is complicated, as described before, and exact computation of the likelihood is difficult. To the best of my knowledge, the rigorous evaluation of the likelihood is not found in the past papers.
Due to the availability of high-frequency data in the financial market, the realized measure has been provided by many authors as the proxy for a latent variable of interest using such data in recent decades. AndersenBollerslev(1998) proposed a daily volatility estimator known as “realized volatility,” the sum of the finitely sampled squared returns over a day. Many pieces of literature showed that it is a very informative measure for volatility. As pointed out by many authors (Bauwens_etal(2012)), however, it is known that high-frequency data is contaminated by noise, and proxies for unobserved components using such data are biased. For example, MartensDijk(2007) and ChristensenPodolskij(2007) discussed the bias of high-low range for the short interval of a day. As a result, many studies addressed this topic and presented various realized measures for accurately estimating unobserved variables of interest. MartensDijk(2007), ChristensenPodolskij(2007), and Christensen_etal(2009) proposed “realized range” volatility estimators based on short interval ranges. Other than this model-free approach, some authors took a model-based approach for estimating the variance of daily asset return using the return and the corresponding realized measure simultaneously. Takahashi_etal(2009) and DobrevSzerszen(2010) extended the SV model framework by adding an observation equation that describes the linkage between realized volatility and unobserved volatility. It is referred as realized stochastic volatility (RSV) model, and related models were studied by KoopmanScharth(2013), Shirota_etal(2014), ZhengSong(2014), and Takahashi_etal(2016). Realized GARCH model was introduced as an extension of the GARCH model to relate latent volatility and the corresponding realized measure, similarly as RSV models (e.g., Hansen_etal(2012), Watanabe(2012), Wang_etal(2019), and Chen_etal(2021)).
Although approaches based on realized measures successfully and accurately estimate the latent volatility of an asset return, DegiannakisLivada(2013) noted that realized volatility is a more accurate estimator than range-based one only when the intraday dataset is available at high sampling frequency. That is, the dependability of these approaches is determined by the transaction frequency of the asset. On the other hand, even if the asset is not liquid and has low transaction frequency, the financial market always has the open, highest, lowest, and closing prices of the asset. Thus, the range-based approach for the estimation of latent volatility is appealing in light of data availability. Additionally, collecting and storing the few transaction prices is a very light burden compared with acquiring and storing the detailed intraday transaction information and constructing the corresponding realized measure from it.
As stated above, there are two major problems about the usage of the range for volatility estimation in the previous studies: (1) leverage effect and (2) calculation of the likelihood function. In this study, I tackle these problems and propose a new framework that models the return and the price range jointly in the context of the SV model with leverage.
The leverage effect is a phenomenon in the stock market: a negative relationship between the return of an asset at date $t$ and the variance of the return at day $t+1$. The Bayesian estimation method of the SV model with leverage was discussed by Jacquier_etal(2004) and Omori_etal(2007). As far as the author knows, past research did not deal with SV model using range data that incorporates the leverage effect.
In this paper, I express the range's density function as the sum of an infinite series other than the one derived by Feller(1951). The two representations enable us to generate samples from the density rigorously without computing the sum of the infinite series. Using a Bayesian approach, I apply this sampling method to the proposed model and implement a generation algorithm for latent volatilities via the Markov chain Monte Carlo (MCMC) method.
The rest of this article is organized as follows. In Section 2, I investigate the density of the range and sampling from the density. Bias correction of the range is discussed, and the SV model with range-based correction and leverage is proposed. Section 3 introduces a Bayesian estimation method for the model that make use of an efficient MCMC algorithm. Section 4 contains an empirical analysis that uses financial return and price range data from the U.S. stock indices. The findings of the paper are presented with a conclusion in Section 5.
I consider a continuous time process,
where $p(s)$ is the logarithm of the price of an asset observed in the financial market at time $s$, $\sigma^2(s)$ denotes the instantaneous volatility, and $W(s)$ is a standard Wiener process. I also assume that $\sigma^2(s)$ is statistically independent of $W(s)$.
The main objective of this study is to estimate the true volatility or integrated volatility for a day $t$ defined as
Let the highest price of the day $t$ be $H_t$, and the lowest price of the day $t$ be $L_t$. I define the price range in day $t$, $r_t$, as the log difference between the highest and the lowest value, i.e., $\log H_t -\log L_t$.
Feller(1951) showed that the density function of range $r_t$ given $\sigma^2_t$ is
and that
where $L(x) =(2\pi)^{\frac{1}{2}} x^{-1} \sum_{n_d=1}^\infty \exp\big\{ -\frac{(2n_d-1)^2\pi^2}{8x^2} \big\}$ is the Kolmogorov-Smirnov distribution function.
Feller(1948) stated that $L(x)$ can be written in two equivalent forms as
It is straightforward that the density function ((ref)) corresponds to the distribution function represented in equation ((ref)). I also obtain another representation of the density as follows:
Summarizing the above discussion, I obtain the next proposition.
{\bf Proposition 1.} \\ The probability density of $r_t$ given $\sigma^2_t$ is expressed as follows: \\ (a) (Feller(1951)) \[ f_{\mathrm{range}}(r_t|\sigma^2_t) =8\sum_{n_d=1}^\infty (-1)^{n_d-1} \frac{n_d^2}{\sqrt{2\pi\sigma_t^2}}\exp\bigg( -\frac{n_d^2r_t^2}{2\sigma_t^2} \bigg). \] (b) \[ f_{\mathrm{range}}(r_t|\sigma^2_t) =8\sum_{n_d=1}^\infty \bigg\{ \frac{(2n_d-1)^2\pi^2\sigma_t^4}{r_t^5} -\frac{\sigma_t^2}{r_t^3} \bigg\} \exp\bigg\{ -\frac{(2n_d-1)^2\pi^2\sigma_t^2}{2r_t^2} \bigg\}. \]
In this subsection, I consider the generation of a random variable with density stated in the last subsection. I develop a random sampling algorithm for $x=r_t^2/\sigma_t^2$ using an alternating series method based on Devroye(1986).
Assume that the density function for a random variable $z$ is represented as $f(z)=ch(z)\sum_{i=0}^\infty a_i(z)$, where $c$ is a constant, a sequence of functions $\{a_i(z)\}_{i=0}^\infty$ ($a_0(z)\equiv 1$) is decreasing, and $h$ is a density function. I also assume that it is easy to generate a random sample from $h$. Then, the following inequality for $f$ and the partial sums $A_j(z) =ch(z)\sum_{i=0}^j a_i(z), \; j=0, 1, \ldots,$ holds: \[ A_0(z) \ge A_2(z) \ge \cdots \ge f(z) \ge \cdots \ge A_3(z) \ge A_1(z). \] A random sample can be generated from density function $f$ using an acceptance-rejection algorithm as shown below.
Notice that the density function for $x=r_t^2/\sigma_t^2$ is represented in two forms:
where $a_{n_d} = (n_d+1)^2 \exp\{-\frac{1}{2}((n_d+1)^2-1)x\}$, and
where
I have the following lemma. (The proof is very similar to that of Lemma 5.1 in Devroye(1986).)
{\bf Lemma 1.}
(a) $a_{n_d}$'s in equation ((ref)) are monotone decreasing for $x>\frac{4}{3}$.
(b) $a_{n_d}$'s in equation ((ref)) are monotone decreasing for $x<\pi^2$.
{\bf Proof:}
For the first series expansion ((ref)), I have
if $x>\frac{4}{3}$. In the second series expansion ((ref)), for $n_d$ even and for $x<\pi^2$, I have
and
$\Box$
This lemma implies that for at least one of the representations the condition stated below is satisfied: \[ A_0(x) \ge A_2(x) \ge \cdots \ge f(x) \ge \cdots \ge A_3(x) \ge A_1(x), \] where $x$ is real and positive. The lemma also shows that the inequality holds true for neither of the representations on the entire support of $f$. It suggests that I can use the alternating series method with
and
where the threshold $c_{\mathrm{th}}$ is in between the range $(\frac{4}{3}, \pi^2)$. Thus, I find that as a proposal $X$ may be sampled from a mixture distribution of chi squared and inverse gamma:
where $p =\int_{c_{\mathrm{th}}}^{\infty} (2\pi)^{-\frac{1}{2}}x^{-\frac{1}{2}}\exp(-\frac{x}{2})dx$, $q =\int_{0}^{c_{\mathrm{th}}} \pi^2 x^{-3} \exp ( -\frac{\pi^2}{2x} ) dx $, and $\mathrm{I}(\cdot)$ is an indicator function.
Summarizing above, I obtain the algorithm that generates a random sample $x$ as follows:
Based on Feller(1951), Parkinson(1980) showed that the $p$-th moment of daily return is represented as
where $\Gamma(\cdot)$ denotes gamma function and $\zeta(\cdot)$ is Riemann zeta function. Note that $\mathrm{E}(r_t) =\sqrt{8\sigma^2_t/\pi}$ and $\mathrm{E}(r_t^2) =(4\log 2)\sigma^2_t$ and that $r_t^2/(4\log 2)$ is often used as a volatility estimator. Contrary to the assumption that the assets are driven by a Brownian motion and have a continuous price, the asset prices observed in the market take discrete values. In addition, the market is not open 24/7 for trading. Thus, several authors suggested that the observed high-low range $r_t$ is regarded as biased.
In this paper, the high-low range $r_t$ observed at the market is assumed to be proportional to the true range $\widetilde{r}_t$ as
where a gamma distribution is denoted as $\mathrm{G}(\alpha, \beta)$. If the scaling parameter $\lambda_t$ is smaller than $1$, it implies that the observed range $r_t$ has a downward bias.
The asymmetric SV model is represented using state space form as
where $y_t$ is a daily asset return, $\sigma^2_t$ is the unobserved volatility, $\mathrm{N}(\bm{m}, S)$ denotes a normal distribution with mean vector $\bm{m}$ and covariance matrix $S$, $\bm{0}_p$ is a $p$-dimensional vector of zeros and $|\phi|<1$. If $\omega_{\eta\epsilon}$ takes a negative value, it means the existence of asymmetry or leverage effect as stated in Section 1.
I add the following observation equation to the state space model described above:
where $f_{\mathrm{range}}$ is the density function defined in Proposition 1. The above Equations ((ref))--((ref)) define the stochastic volatility model with range-based correction and leverage (SVRG).
Let $\widetilde{\sigma}^2_t =\lambda_t\sigma^2_t$ and notice that $r_t \sim f_{\mathrm{range}}(r_t |\widetilde{\sigma}^2_t)$. Thus, I may rewrite the proposed SVRG model as follows:
I assume the prior of $\phi$ that
where a beta distribution with parameters $(a,b)$ is denoted as $\mathrm{Be}(a,b)$.
For $\Omega$, let
and notice that $\omega^{(\epsilon\epsilon)} =1 +(\omega^{(\epsilon\eta)})^2 (\omega^{(\eta\eta)})^{-1}$. Thus, I assume the prior distribution for $\Omega^{-1}$ such that
Let $\bm{\nu}=(\nu_1, \nu_2)'$ and an assumption for the prior of $\bm{\nu}$ is made as follows:
Define $\bm{\vartheta}=(\phi, \Omega, \bm{\nu})$ and let us denote the prior probability density of $\bm{\vartheta}$ as $\pi(\bm{\vartheta})$. I obtain the joint posterior probability density that is expressed as
where
The MCMC algorithm is implemented in five blocks:
The conditional posterior density function of $\sigma^2_t$ is given by
Let $h_{\sigma^2_t} =\log\sigma^2_t$. Using Taylor expansion around, say, $c_{\sigma^2_t}$, I obtain an approximation for $\sigma_t^{-1}$ as
It is recommended to set $c_{\sigma^2_t} =\log\{ r_t^2/(4\lambda_t\log 2) \}$. For $t=1, \ldots, n-1,$ I may consider an approximation of $f(\sigma^2_{t+1} |y_t, \sigma^2_{t}, \bm{\vartheta})$ as
where
Consider the following function,
and note that $g_{\sigma^2_t}$ is approximated by a function proportional to the probability density of log-normal distribution as shown below:
Let $\alpha_{\sigma^2_t} =(\exp(s_{\sigma^2_t})-1)^{-1}+2$, and $\beta_{\sigma^2_t} =\exp(m_{\sigma^2_t}+\frac{1}{2}s_{\sigma^2_t}) \{(\exp(s_{\sigma^2_t})-1)^{-1}+1\}$. Notice that a random variable with the inverse gamma density function with parameters $(\alpha_{\sigma^2_t}, \beta_{\sigma^2_t})$ has mean $m_{\sigma^2_t}$ and variance $s_{\sigma^2_t}$. It implies that the log-normal probability density function with (log) mean $m_{\sigma^2_t}$ and (log) variance $s_{\sigma^2_t}$ is approximated by inverse gamma distribution with parameters $(\alpha_{\sigma^2_t}, \beta_{\sigma^2_t})$.
Then, an approximation of $f(\sigma^2_{t} |\{y_s\}_{s=1}^n, \{r_s\}_{s=1}^n, \{\sigma^2_{s}\}_{s\neq t}, \{\lambda_s\}_{s=1}^n, \bm{\vartheta})$ is given by the probability density function $f_a(\sigma^2_{t} |\{y_s\}_{s=1}^n, \{\widetilde{r}_s\}_{s=1}^n, \{\sigma^2_{s}\}_{s\neq t}, \bm{\vartheta})$, where
and where
The right hand of equation ((ref)) is rewritten as
where
I notice that $f_a(\sigma^2_t|\cdot)$ is represented using a constant $c$, an easy density function $h$ and a series $b_{n_d}$ as $f_a(\sigma^2_t|\cdot)=ch(\sigma^2_t|\cdot)\sum_{n_d=0}^\infty (-1)^{n_d} b_{n_d}$, where
or
I have the following lemma.
{\bf Lemma 2.}
(a) $b_{n_d}$'s in equation ((ref)) are monotone decreasing for $\sigma^2_t<\frac{3}{4}\widetilde{r}_t^2$.
(b) $b_{n_d}$'s in equation ((ref)) are monotone decreasing for $\sigma^2_t>\frac{\widetilde{r}_t^2}{\pi^2}$.
I may prove this lemma exactly the same as Lemma 1.
This lemma implies that for at least one of the representations the inequality stated below holds true: \[ B_0(\sigma^2_t) \ge B_2(\sigma^2_t) \ge \cdots \ge f_a(\sigma^2_t) \ge \cdots \ge B_3(\sigma^2_t) \ge B_1(\sigma^2_t), \] where $B_j(z) =ch(z)\sum_{i=0}^j b_i(z)$. The lemma also shows that for neither of the representations the inequality holds true on the entire support of $f_a$. This suggests that the alternating series method can be used with
and
where the threshold $c_{\mathrm{th}}$ is in between the range $(\pi^{-2}\widetilde{r}_t^2, \frac{3}{4}\widetilde{r}_t^2)$. Thus, I find that as a proposal $\sigma^2_t$ may be sampled from a mixture distribution of inverse gamma and generalized inverse Gaussian:
Note that
where $\Gamma(\alpha, x)$ is an incomplete gamma function defined as
and
where $K_\nu(x,y)$ denotes an incomplete Bessel function defined as
A generalized inverse Gaussian distribution with parameters $(\nu, \delta, \gamma)$ is denoted as $\mathrm{GIG}(\nu, \delta, \gamma)$. These incomplete integrals appeared in $p_t$ and $q_t$ are calculated numerically (see Harris(2008)). An efficient sampling algorithm from a truncated gamma distribution is available from Philippe(1997). I generate a sample from a (truncated) generalized inverse Gaussian density using the algorithm proposed by Dagpunar(1989).
Now an algorithm for generating $\sigma^2_t$ is obtained as follows:
The conditional posterior density function of $\lambda_t$ is derived as
Let us denote $\log\lambda_t$ by $h_{\lambda_t}$. I derive the Taylor expansion of $\lambda_t^{\frac{1}{2}}$ with respect to $h_{\lambda_t}$ around, say, $c_{\lambda_t}$, and obtain the approximation as
It is recommended to use the prior mode $\log\{ (\frac{\nu_1}{2}-1)/(\frac{\nu_2}{2}) \}$ as $c_{\lambda_t}$. Thus, for $t=1, \ldots, n-1,$ an approximation of $f(\widetilde{\sigma}^2_{t+1} |y_t, \lambda_{t+1}, \widetilde{\sigma}^2_{t}, \lambda_t, \bm{\vartheta})$ may be considered as
Define the following function,
and note that I can approximate $g_{\lambda_t}$ to a function proportional to the probability density of log-normal distribution:
I can show that a random variable with the gamma density function with parameters $(\alpha_{\lambda_t}, \beta_{\lambda_t})$, has mean $m_{\lambda_t}$ and variance $s_{\lambda_t}$, where $\alpha_{\lambda_t} =(\exp(s_{\lambda_t})-1)^{-1}, \;\; \beta_{\lambda_t} =\exp(-m_{\lambda_t}-\frac{1}{2}s_{\lambda_t}) (\exp(s_{\lambda_t})-1)^{-1}$. I observe that the log-normal density function is approximated by the gamma distribution with parameters $(\alpha_{\lambda_t}, \beta_{\lambda_t})$.
Then, $f(\lambda_{t} |\{y_s\}_{s=1}^n, \{r_s\}_{s=1}^n, \{\widetilde{\sigma}^2_{s}\}_{s=1}^n, \{\lambda_{s}\}_{s\neq t}, \bm{\vartheta})$ is approximated to the probability density function $f_a(\lambda_{t} |\{y_s\}_{s=1}^n, \{r_s\}_{s=1}^n, \{\widetilde{\sigma}^2_{s}\}_{s=1}^n, \{\lambda_{s}\}_{s\neq t}, \bm{\vartheta})$, where
I generate a candidate $\lambda_t^\dagger \sim \mathrm{G}(\frac{\alpha_{t, 1}}{2}, \frac{\beta_{t, 1}}{2})$ and it is accepted with probability
where $\lambda_t$ is a current MCMC sample.
{\it Generation of $\phi$.} The conditional posterior density function of $\phi$ is given by
where
Sample $\phi^\dagger\sim\mathrm{TN}_{(-1,1)}(m_\phi, s_\phi)$ as a candidate, where a normal distribution with mean $m_\phi$ and variance $s_\phi$ truncated to the region $(-1, 1)$ is denoted as $\mathrm{TN}_{(-1,1)}(m_\phi, s_\phi)$. It is accepted with probability $\min\{1, g(\phi^\dagger)/g(\phi)\}$, where $\phi$ is a current MCMC sample.
{\it Generation of $\Omega$.} I can obtain the conditional posterior density of $(\omega^{(\epsilon\eta)}, \omega^{(\eta\eta)})$ as follows:
where
(For the multivariate case, see Ishihara_etal(2016).) $\Omega^\dagger$ is generated as a candidate as the following:
{\it Generation of $\bm{\nu}$.} For $\bm{\nu}$, the conditional posterior density is given by
I transform $\bm{\nu}$ such that $\widetilde{\bm{\nu}}=(\log\nu_1, \log\nu_2)'$. Taylor expansion around the mode $\widehat{\bm{\nu}}$ gives the following approximation for the conditional posterior density of $\widetilde{\bm{\nu}}$ as
where
I generate a candidate $\widetilde{\bm{\nu}}^\dagger$ from the normal density $f_{\mathrm{N}}$ with mean vector $\bm{m}_{\bm{\nu}}$ and covariance matrix $S_{\bm{\nu}}$ and it is accepted with probability
This section provides empirical results using asset returns and price ranges for two stock indices in the U.S. financial market: Dow Jones Industrial Average (DJIA) index and the S&P500 index. The asset return $y_t$ is defined as $(\log p_{t} -\log p_{t-1})\times 100$, where $p_t$ is the asset's closing price at day $t$. The price range $r_t$ at day $t$ is calculated as $(\log H_t -\log L_t)\times 100$, where $H_t$ is the highest and $L_t$ is the lowest asset prices at date $t$. For model comparison, I use corresponding five-minute sub-sampled realized volatility data obtained from Oxford-Man Institute's “Realized Library” (Heber_etal(2009)). The sample period is from January 3, 2012 to December 31, 2020 (2,254 observations for DJIA and 2,256 observations for S&P500, respectively).
I apply the return and price range dataset described in the last subsection to the proposed model with the following priors for the parameters:
Using the algorithm discussed in Section 3, 10,000 MCMC samples are generated after discarding the initial 1,000 samples as the burn-in.
Table (ref) reports the posterior means, 95% credible intervals and estimated inefficiency factors (IF). The inefficiency factor is defined as $1+2\sum_{k=1}^\infty \rho(k)$, where $\rho(k)$ is the autocorrelation function for the parameter at lag $k$ (see Chib(2001)).
For DJIA (S&P500) data, two top panels of Figure (ref) ((ref)) depicts the time series plot of $y_t$ and $r_t$, respectively. The transition for the estimated posterior means of $\sigma_t$ is illustrated in the third top panel of Figure (ref) ((ref)). The bottom panel of Figure (ref) ((ref)) reports the posterior means and the bounds for 95% credible intervals of $\lambda_t$.
The posterior means of the first order autoregressive parameter $\phi$ are over 0.9 and very high, which indicates that the processes of (log) volatilities are strongly persistent. It is consistent with the plot of the posterior means of $\sigma_t$'s recorded in Figures (ref) and (ref). I observe the sharp increase of the trajectories of $\sigma_t$ in March 2020, resulted from COVID-19 pandemic. The trajectories increased in August 2015, corresponding to the global stock market sell-off while those increases in February 2018 and late 2018 are due to the spike in VIX and cryptocurrency crash, respectively.
I notice that the 95% credible intervals of $\omega_{\epsilon\eta}$ do not include $0$. For DJIA (S&P500), the posterior mean of correlation coefficient between $y_t$ and $\log\sigma^2_{t+1}$ is equal to $-0.449$ ($-0.468$). It strongly suggests the existence of asymmetry or leverage effect and is consistent with reports in the vast amount of papers.
When I use DJIA (S&P500) data, the sample means of 95% credible interval upper and lower bounds of $\lambda_t$'s are 1.037 (1.080) and 0.568 (0.606), respectively. I see that the 95% credible intervals of the bias correction coefficients for the price range are large and that many of them include the value of 1. It implies that the measured price range contains a large amount of noise and that the downward bias for the range revealed in the previous studies is not necessarily significant.
For DJIA (S&P500), the acceptance rates of $\{\sigma^2_t\}_{t=1}^n$, $\{\lambda_t\}_{t=1}^{n}$, $\bm{\nu}$, $\phi$ and $\Psi$ in the independent chain Metropolis-Hastings algorithm are 0.952 (0.942), 0.957 (0.950), 0.985 (0.983), 0.984 (0.994) and 0.993 (0.993), respectively. I notice that the proposal density functions approximate the posterior accurately and that the MCMC algorithm can generate samples efficiently.
{\it Generation of volatility from posterior predictive distribution.} The posterior predictive density of $\log\sigma^2_{n+1}$ is given by
I add a step of generation of $\sigma^2_{n+1}$ to each MCMC iteration described in the last section and obtain the estimated posterior predictive density.
{\it Model comparison based on robust loss function.} Suppose that $L(\sigma^2_t, h_{t})$ denotes a loss function, i.e., a measure for error between the true volatility of asset return at $t$-th day $\sigma^2_t$ and a forecast $h_{t}$. In addition, suppose that I have two volatility forecasts $h_{j,t}$ and $h_{k,t}$ based on two statistical models $j$ and $k$. $h_{j,t}$ outperforms $h_{k,t}$ with respect to $L$ if the expectation of the loss function is smaller than the other one, i.e., $\mathrm{E}\{L(\sigma^2_t, h_{j,t})\} < \mathrm{E}\{L(\sigma^2_t, h_{k,t})\}$. Needless to say, the true volatility is not observed in the financial market, and I have to use a volatility proxy. Since a volatility proxy contains noise, it is not always true that the ranking between two volatility forecasts $h_{j,t}$ and $h_{k,t}$ is preserved if one substitute a conditionally unbiased proxy $\widehat{\sigma}^2_t$ to the true volatility, i.e.,
By Patton(2011), it was shown that a class of loss functions holds equivalence of such ordering between two volatility forecasts. The following two loss functions mean squared error (MSE) and quasi-likelihood (QLIKE) are included in Patton(2011)'s class:
Note that these take value of zero if and only if $h_t =\widehat{\sigma}^2_t$. Patton(2011) used realized measure and Parkinson's range-based estimator $r_t^2/(4\log 2)$ as conditionally unbiased volatility proxies.
It is well known that realized measure and range-based estimator for variance of asset return often take smaller value than the true due to the market microstructure noise or overnight effect as stated in Section 1. HansenLunde(2006) recommended that realized volatility $RV_t$ should be scaled using daily returns as
It implies that the mean of scaled realized volatility $RV_t^{\mathrm{scale}}$ takes the same value of the sample variance of corresponding daily returns.
In this subsection, I use realized volatility (RV) and Parkinson's range-based estimator (RG) as volatility proxies which are scaled as HansenLunde(2006). Posterior predictive mean of $\sigma^2_{t+1}$ is adopted as the forecasts. Thus, I may obtain the ranking between two volatility forecasts, which is consistent with that using true volatility with respect to the loss functions described above.
Although the ranking provided as above is robust to the noise in the volatility proxy and is reliable, one may notice that the difference between the two losses is not measured. I use a predictive ability test proposed by GiacominiWhite(2006) to determine whether the difference is statistically significant. I utilize constant and lagged loss difference for the test function and conduct the test in line with the theorem in that publication.
{\it Empirical results.} The rolling forecast is conducted as follows: (i) Apply the estimation algorithm to the data from 2012 to 2018 to (I set the data as $\{y_t\}_{t=1}^n$ and $\{r_t\}_{t=1}^n$) and generate the parameters and unobserved variables as well as one-step-ahead forecast of variance, $\sigma^2_{n+1}$. (6,000 MCMC samples are generated after discarding 1,000 samples from the iterations.) (ii) The new data observed at the followed day is added to the period and the data observed at the oldest day is removed. (iii) Generate the model parameters and unobserved variables for the data (relabeled as $\{y_t\}_{t=1}^n$ and $\{r_t\}_{t=1}^n$). Further, generate $\sigma^2_{n+1}$. (iv) Go to Step (ii). In the end, I obtain posterior predictive means of $\sigma^2_{n+1}$ (499 in total for DJIA data and 497 in total for S&P500 data, respectively).
As a benchmark model, I employ RSV model (Takahashi_etal(2009)). Takahashi_etal(2021) conducted model comparison between RSV model and the other competing volatility models and show that it has an excellent predictive ability. Details for this benchmark model are stated in Appendix. I use daily asset return data and associated realized measure of variance as mentioned in Section 5.1 for this competing model.
Tables (ref) and (ref) report the average values of the MSE and QLIKE loss functions of the forecasts. I use RV and RG as the volatility proxies. Tables (ref) and (ref) also show the $p$-values of Giacomini and White's test with the null hypothesis that the predictive ability is equal to that of SVRG model.
For DJIA data, both MSE and QLIKE are lower for SVRG model than for RSV model. It indicates that SVRG model ranks better than RSV model, whereas the results of Giacomini and White's test show that significant difference is not found between the two models at significant level of 0.01. For S&P500 data, Table (ref) reports similar empirical results as DJIA data. In summary, I conclude that the predictive ability of the model is equivalent to or better than that of RSV model. Without high frequency data in the financial market, I may obtain an accurate forecast of volatility of an asset return by using the proposed model.
This paper introduces a novel SV model in the state space form by adding new observation equation which relates the variance and the price range. I investigate the probability distribution of range for daily stock exchange and the rigorous sampling algorithm from the density function. An MCMC algorithm is developed based on the Bayesian framework to obtain estimation results for the proposed SVRG model.
An empirical study is presented using daily returns and price ranges of DJIA and S&P500 indices in the U.S. stock market. I find the persistence in the (log) volatility processes, the existence of the leverage effect, and a non-negligible amount of noise in the observed price range as suggested by the previous empirical papers. The model comparison is based on out-of-sample predictive performance and the results reveal that the forecasts generated by the model are as accurate (if not more so) as those of RSV model.
Although the proposed univariate model can be extended to a multivariate model, building such models presents various challenges. Many financial market researchers have found dynamic correlations between asset returns, which should be considered in the multivariate volatility model. I must guarantee the covariance matrix of multivariate asset returns is positive definite even if dynamic correlation structure is incorporated. Furthermore, there are usually too many model parameters and latent variables to obtain stable estimation results for such multivariate or high-dimensional time-varying covariance structure models. In the context of SV models, e.g., KuroseOmori(2016) proposed parsimonious modeling of multivariate asset returns and it would be worthy of considering such extension. This is beyond the scope of this paper and the author leaves the extension to the future study.
The author wishes to thank Kaoru Irie, Tsunehiro Ishihara, Yasuhiro Omori, and Makoto Takahashi for the comments. This research was supported by JSPS KAKENHI Grant Numbers 26245028, 19H00588, 20K19751. The computational results in this article were generated using Ox (see Doornik(2009)).