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.
59,238 characters · 19 sections · 45 citation commands
Stochastic Volatility in Mean: Efficient Analysis by a Generalized Mixture Sampler
{\bf JEL classification}: C11, C15, C32, C58 \\ {\bf Keywords}: Excess Holding Yield; Markov chain Monte Carlo; Mixture Sampler; Risk Premium; Stochastic Volatility in Mean
In financial time series, volatility clustering, the phenomenon of persistent volatility, is well known to exist. One way to model time-varying volatility, or volatility clustering, is by using the stochastic volatility (SV) model of Taylor(08). In the simplest version of this model, the standard deviation of the outcome is given by an exponential transformation of an unobserved log-variance variable $h_t$, where $h_t$ in turn is modeled by a stationary first-order autoregressive process. The model in this basic form can be viewed as a state-space model in which the measurement equation is nonlinear in the latent variance $h_t$. A significant variant of this standard SV model is one in which the standard deviation of the outcome $\exp(h_t/2)$ appears as a predictor variable in the mean of the measurement equation. This is called the SV in mean model (SVM) and is similar in spirit to the ARCH in mean (ARCH-M) model introduced by EngelLilienRobins(87). Like the standard SV model, the SVM model has also been used in various fields, including macroeconomics and finance (see KooopmanHol(02), BerumentYalcinYildirim(09), MumtazZanetti(13), CrossHouKoopPoon(23)).
In the Markov chain Monte Carlo (MCMC) estimation of model parameters for the SV-type models, it is often observed that sampling one latent variable conditional on all the other latent variables and parameters, which is referred to as the single-move sampler, is inefficient in the sense that MCMC draws are highly autocorrelated. To address this issue, KimShephardChib(98) introduced the mixture sampler as a highly efficient Bayesian estimation method for the standard SV model. This approach was extended to SV models with jumps and fat-tailed errors in ChibNardariShephard(02) and to SV models with leverage in OmoriChibShephardNakajima(07).
In the existing literature, the SVM model without leverage has been estimated by the multi-move (block) sampler (ShephardPitt(97) and OmoriWatanabe(08)), for example, AbantoMigonLachos(11), AbantoMigonLachos(12), and LeaoAbantoChen(17), and by other similar approaches, for example, Chan(17), AbantoRodriguezGarrafa(21), and AbantoRodriguezHernan(23). There is no known mixture sampler approach for SVM models with leverage. In this paper, we develop efficient MCMC based algorithms for SVM models, with and without leverage, that are based on accurate representations of these models in terms of mixtures of conditionally Gaussian linear state-space models, just as in the approach of KimShephardChib(98). However, due to the $\beta\exp(h_t/2)$ term in the mean equation, instead of characterizing the distribution of a central chi-squared distribution with one degrees of freedom in terms of mixtures of normal distributions, a mixture representation to the distribution of a non-central chi-squared distribution with one degrees of freedom is needed, in which the non-centrality parameter is $\beta^2$. We show that the latter distribution has an infinite series expansion. Our estimation approach uses a truncated version of this series expansion to develop a highly efficient fitting algorithm. It can be viewed as a generalized mixture sampler. Furthermore, the small approximation error due to the truncation can be corrected by a data augmentation method by incorporating a pseudo-target probability density whose marginal probability density is the exact conditional posterior density.
We apply our proposed method to excess holding yields and S&P500 returns data, and show that SVM models outperform other SV models based on marginal likelihoods. The rest of this paper is organized as follows. Section 2 introduces the SVM model and describes the novel mixture sampler as an efficient sampling method for such models. The MCMC simulation and the particle filter are described in Section 3. Section 4 illustrates the performance of this sampling method using the simulated data for several cases. Section 5 further extends the SVM model to incorporate the leverage effect. Finally, in Section 6, we apply our proposed SVM model to financial data and perform a model comparison. Conclusion and remarks are given in Section 7.
We define the stochastic volatility in mean (SVM) model as follows:
where $N(m,S)$ denotes normal distribution with mean vector $m$ and covariance matrix $S$, $\theta = (\mu, \phi, \sigma^2, \beta)$ is a model parameter vector of interest, and $h = (h_1, \dots, h_n)'$ is the logarithm of the latent volatility vector. We use the standard deviation $\exp(h_t/2)$ in the mean equation, rather than the variance $\exp(h_t)$, to match the units of the outcome variable. For $\beta \neq 0$, we denote it as the stochastic volatility in mean (SVM) model. The standard stochastic volatility (SV) model is obtained as a special case with $\beta =0$. For $\theta$, we assume the prior distribution
where $Beta(a,b)$ denotes beta distribution with parameters $(a,b)$ and $IG(a,b)$ denotes inverse gamma distribution with parameters $(a,b)$ whose probability density function is
We let $f(y, h|\theta)$ and $\pi(\theta)$ denote the probability density function of $(y, h)$ given $\theta$ and the prior probability density function of $\theta$ where $y \equiv (y_1, \dots, y_n)'$. The posterior density function of $(h, \theta)$ is given in Appendix (ref). \\
{\it Remark}. It is straightforward to include the constant term in the measurement equation. However, noting that $\beta \times \exp(h_t/2) \approx \beta \times (1 + h_t/2)$, it is often confounded with $\beta$ and therefore omitted in this paper.
To sample $h$ from its conditional distribution, we transform Equation ((ref)) as below:
Since $(\beta + \epsilon_t) \sim N(\beta, 1)$, its square $(\beta + \epsilon_t)^2 \sim \chi_1^2(\beta^2)$ where $\chi_1^2(\beta^2)$ denotes the non-central chi-square distribution with the non-centrality parameter $\beta^2$ and one degrees-of-freedom. The special case with $\beta=0$ is considered in KimShephardChib(98) and OmoriChibShephardNakajima(07), who introduced the idea of accurately approximating the probability of the logarithm of the central chi-square distribution with one degrees of freedom, $\log \chi_1^2(0)$,
by a mixture of normal distributions. Below, we elaborate a highly accurate approximation of the distribution of $\epsilon_t^*$ given $\beta \neq 0$ by the mixture of normal distributions. Let $p(x;\nu, \lambda)$ be the probability density function of $\chi_\nu^2(\lambda)$. It can be expressed as an infinite mixture of central $\chi^2$ probability density functions (see, e.g. JohnsonKotzBalakrishnan(95)):
Setting $\nu=1$ and noting that
we obtain the expression
Let $f(u;\lambda)$ denote the probability density function of $U \sim \log \chi_1^2(\lambda)$. Using ((ref)), it follows that
As in OmoriChibShephardNakajima(07), we consider the mixture of ten normal distributions to approximate $f(u;0)$, the probability density function of $\log \chi_1^2(0)$,
where $\phi(\cdot)$ denotes the probability density function of the standard normal distribution. The values of $(p_i, m_i, v_i^2)$ are taken from OmoriChibShephardNakajima(07) and are reproduced in Table (ref). The columns labeled $a_i$ and $b_i$ will be used when we consider the model with leverage in Section (ref).
By substituting Equation ((ref)) to Equation ((ref)), we obtain
where
and $\Tilde{m}_{i,j} = m_i+j v_i^2$. In the last equality, we truncate the summation at $j=J$ and normalize $\Tilde{p}_{i,j}$ to ensure that $\sum_{i=1}^{K} \sum_{j=0}^{J} \Tilde{p}_{i,j} = 1$.
This expression implies that the probability density function of $\log \chi_1^2(\lambda)$ is approximated by the mixture of $K(J+1)$ normal distributions. Especially, when $\lambda = 0$ and $J=0$, the approximation ((ref)) reduces to ((ref)). The extent of the impact on the approximation is based on the value $w_{i,j}$ which includes $\lambda=\beta^2$. The coefficient $\beta$ of volatility $\exp(h_t/2)$ is estimated to be less than one in past empirical studies. When $\lambda = \beta^2$ is less than 1.0, $J = 1$ or $2$ makes $\sum_{i=1}^{K} \sum_{j=0}^{J} w_{i,j}$ more than 0.9. For $J=2$, Figures (ref) and (ref) show the true and approximate densities of $\log \chi_1^2(\beta^2)$ for $\beta = 0.3, 0.5$ and $0.7$ (equivalently, $\beta = -0.3, -0.5$ and $-0.7$) and the difference between the two densities, respectively. Since these differences are quite small and the approximation almost overlaps the true probability density of $\log \chi_1^2(\beta^2)$ even for $\beta=0.7$, we employ $J = 2$ in this paper. That is, we approximate the probability density of $\epsilon_t^*|\beta \sim \log \chi_1^2(\beta^2)$ by
where
Let $s_t = (s_{1t}, s_{2t}) \in \{ (i,j) | i=1,...,K, j=0,...,J \}$ denote the component of the mixture of the normal densities in ((ref)) at time $t$. Given $s_t = (i,j)$, we have $\epsilon_t^* | s_t=(i,j) \sim N(\Tilde{m}_{i,j}, v_i^2)$ and we see that the SVM model can be approximated by the linear Gaussian state space form
where $y_t^* = \log(y_t^2)$. In the following sections, $\Tilde{m}_{s_{1t}, s_{2t}}$ and $\Tilde{p}_{s_{1t}, s_{2t}}$ are abbreviated as $\Tilde{m}_{s_t}$ and $\Tilde{p}_{s_t}$, and we write $v_{s_t}^2$ instead of $v_{s_{1t}}^2$.
Algorithm 1 (Generalized mixture sampler, GMS). Let us denote $\theta = (\alpha,\beta)$ where $\alpha = (\mu,\phi,\sigma^2)$. The Markov chain Monte Carlo simulation is implemented in four blocks:
The conditional posterior distribution of $\beta$ is normal with mean $b_1$ and variance $B_1$ where
and
Thus we generate $\beta \sim N(b_1,B_1).$
As discussed in the previous section, we sample $h$ using the mixture sampler using the mixture of normal distributions. Since our approximation is highly accurate, we can use this mixture approximation directly in Step 3, as in Step 2 of Algorithm 1 in ChibNardariShephard(02). However, one can remove the small approximation error with an additional MH step, as detailed in Algorithm 2, GMS with MH algorithm (GMH), given in the Appendix (ref), but due to the fact that the tailored mixture very closely fits the non-central chi-squared distribution, this additional step would be rarely necessary.
Let $f_N(\cdot|m, s^2)$ denote the probability density of $N(m, s^2)$, and let $\pi(\alpha)$ denote the prior density of $\alpha$. Define our target density in Step 3 as
where
and $\Tilde{m}_{s_t} = \Tilde{m}_{s_{1t}, s_{2t}}$ and $\Tilde{p}_{s_t} = \Tilde{p}_{s_{1t}, s_{2t}}$ are defined in ((ref)) and $(p_{s_{1t}}, m_{s_{1t}}, v_{s_{1t}}^2)$ are given in Table (ref). Note that $m(y^*|\alpha,s,\beta)$ is a normalizing constant for $\pi^*(h|\alpha, s, \beta, y^*)$ and is evaluated using the Kalman filter algorithm. Only $\Tilde{p}_{s_t}$, which depends on $\beta$, needs to be updated according to the formula in ((ref)) before sampling. Note that the target density $\pi^*(\alpha,h|\beta,y)$ approximates the true conditional density $\pi(\alpha,h|\beta,y)$ accurately. We generate the sample $(\alpha, h, s)$ in two steps.
We describe how to compute the likelihood $f(y|\theta)$
numerically as it is necessary to obtain the marginal likelihood, $f(y) = \int f(y|\theta) \pi(\theta) d\theta$ and Bayes factor for the model comparison. The filtering and the associated particle computations are carried out by the auxiliary particle filter (see e.g. PittShephard(99), OmoriChibShephardNakajima(07)). Let us denote $Y_t = (y_1, \dots, y_n)$, and
and consider the importance function for the auxiliary particle filter
where
This leads to the following particle filtering.
It can be shown that as $I \rightarrow \infty$, $\overline{w}_{t+1} \xrightarrow{p} f(y_{t+1}| Y_t, \theta)$, and $\overline{W}_{t+1} \xrightarrow{p} F(y_{t+1}| Y_t, \theta)$. Therefore, it follows that
is a consistent estimate of the conditional log-likelihood and can be used as an input in the calculation of the marginal likelihood by the method of Chib(95), as extended to M-H sampler output by ChibJeliazkov(01).
This section illustrates our proposed estimation method using the simulated data. We generate $y_t$ $(t = 1, \dots, 1000)$ by setting
To avoid the case $y_t = 0$ which leads to $\log (y_t^2) = -\infty$, we introduce very small value $c$ and use $y_t^* = \log (y_t^2 + c)$. We set $c$ equal to $1.0 \times 10^{-7}$. For $\beta$, we consider three cases $\beta = 0.3, 0.5$ and $0.7$ to investigate the effect of the approximation error. The common random numbers are used to generate $y_t$'s for three cases. In these simulation studies, we specify the prior as
The prior on $(\phi+1)/2$ is set to ensure the stationarity of the latent volatility process. We iterated MCMC simulation 50,000 times after discarding initial 10,000 MCMC draws as burn-in period.
These simulation results show that our proposed sampling method works well for those $\beta$'s found in the past empirical studies.
Next, we compare the sampling efficiency of our proposed method with the method in Chan(17), which was shown to be more efficient than the more complex methods of Mccausland(12) and Andrieu(10). Chan's method needs to implement the accept-reject Metropolis-Hastings (ARMH) algorithm ChibGreenberg95 (the usual MH algorithm could be used, but then the rejection rate becomes much higher), using a proposal distribution based on a multivariate normal approximation of the full conditional distribution, obtained by a second-order Taylor expansion.
Table (ref) shows the IFs for selected values ($h_t$, $t=100,200,\ldots,1000$), the mean ($\overline{h}$) and the median ($h_{med}$) of $h_t$ for $\beta=0.3, 0.5$ and $0.7$. Among the three algorithms, the generalized mixture sampler (Algorithm 1, denoted by GMS) is the most efficient with IFs less than 10. These are substantially smaller than those of the CHN method. The IFs for Algorithm 2, the generalized mixture sampler with an additional Metropolis-Hastings step (denoted as GMH), are somewhat larger, but still smaller than those from Chan's method (denoted by CHN). We note that the IFs of GMH become larger as the absolute value of $\beta$ increases. However, as shown in the results of the simulation experiment above and in Appendix (ref), the resulting posterior estimates by GMS and GMH are almost the same. The takeaway from this experiment is the striking efficiency of the GMS method.
In addition, we compare the computational times required in each experiment for the three algorithms. As shown in Table (ref), our proposed methods (GMS and GMH) are much faster than Chan's method (CHN). This is mainly because (i) GMS and GMH sample $h$ using the highly efficient and fast simulation smoother, whereas the CHN method is based on an expensive tailoring step that entails the inversion of a high dimensional covariance matrix in each MCMC iteration; (ii) the AR step of the ARMH algorithm degrades with many rejections when $h$ is high-dimensional and $\beta$ is large; and (iii) the search for the mode in the tailoring step by an iterative optimization methods consumes considerable time for longer time series (i.e., as the dimension of $h$ increases). Our proposed methods encounters none of these difficulties.
In this section, we consider the SVM model with leverage which we call SVML model. The leverage effect implies the decrease in the return at time $t$ followed by the increase in the volatility at time $t+1$. Thus we incorporate the correlation $\rho$ between $y_t$ and $h_{t+1}$ and replace ((ref)) by
The negative correlation, $\rho<0$, indicates the existence of the leverage effect. Next, we construct the linear and Gaussian state space model that approximates the SVM model with leverage using the mixture of the normal densities given in ((ref)). We first let $d_t = I(y_t \geq 0) - I(y_t < 0)$ where $I(A) = 1$ if $A$ is true and $I(A) = 0$ otherwise. Noting that
we rewrite the conditional distribution as
Let $s_t = (s_{1t}, s_{2t}) \in \{ (i,j) | i=1,...,K, j=0,...,J \}$ denote the component of the mixture of normal densities in ((ref)) at time $t$. Given $s_t = (i,j)$, we have $\epsilon_t^* | s_t=(i,j) \sim N(\Tilde{m}_{i,j}, v_i^2)$. Furthermore, we approximate $\exp(\epsilon_t^*/2)$ by $\exp(\Tilde{m}_{i,j}/2) \{ a_i + b_i(\epsilon_t^* - \Tilde{m}_{i,j}) \}$ with $a_i = \exp(v_i^2/8)$, $b_i = \frac{1}{2} \exp(v_i^2/8)$, as in Table (ref), which minimize the mean square norm
Thus, the approximate conditional distribution is
Given $s = (s_1, ..., s_n)$, we find that the SVM model can be approximated by the linear Gaussian state space form
where $y_t^* = \log(y_t^2)$ and $d_t = I(y_t \geq 0) - I(y_t < 0)$. In the following subsections, $\Tilde{m}_{s_{1t}, s_{2t}}$ and $\Tilde{p}_{s_{1t}, s_{2t}}$ are abbreviated as $\Tilde{m}_{s_t}$ and $\Tilde{p}_{s_t}$, and we write $v_{s_t}^2, a_{s_t}$, and $b_{s_t}$ instead of $v_{s_{1t}}^2, a_{s_{1t}}$, and $b_{s_{1t}}$, respectively. The MCMC algorithm and the particle filter are detailed in Appendix (ref).
Furthermore, we will give a couple of examples to illustrate how to extend our proposed SVM model to the multivariate models. As a first example, consider the factor multivariate stochastic volatility (MSV) model proposed by ChibNardariShephard(06) (see IshiharaOmori(17) for the model with leverage). Let $\bm{y}_t=(y_{1t},\ldots,y_{pt})'$ and $\bm{f}_t=(f_{1t},\ldots,f_{qt})'$ denote the dependent and factor variables ($q<p$) respectively. The basic factor MSV model is given by
where $\mathbf{B}$ is a $p\times q$ factor loading matrix subject to constraints ($b_{ij}=0$ for $j>i$ and $b_{ii}=1$ for $i \leq q$). To incorporate the SVM in the MSV model, we replace the factor equation by
where $\bm{\beta}=(\beta_1,\ldots,\beta_q)'$ is a $q\times 1$ vector of weights.\\ As a second example, we may consider the MSV model with the SVM that is common within a group. Let $\bm{y}_t=(y_{1t},\ldots,y_{pt})'$ and consider
where $\bm{\beta}=(\beta_1,\ldots,\beta_p)'$. To linearize the measurement equation with respect to $\tilde{h}_t$, we define $\bm{y}^*=(y_{1t}^*,\ldots,y_{pt}^*)'$ and $\bm{\epsilon}_t^*=(\epsilon_{1t}^*,\ldots, \epsilon_{pt}^*)'$ where $y_{it}^*=\log y_{it}^2$ and $\epsilon_{it}^*=\log (\beta_i/\sigma_i+\epsilon_{it}/\sigma_i)^2$. Thus the transformed measurement equation is
Noting that $\epsilon_{it}^*\sim \log \chi_1^2(\beta_i^2/\sigma_i^2)$, we can approximate the noncentral chisquare distribution by the mixture of normal distributions as we have described.
{\it Remark}. Some previous studies considered alternative mean specifications of $y_t$ using $\exp(h_t)$ or $h_t$, as in EngelLilienRobins(87) for ARCH-M model. We note that our proposed sampler could be utilized to generate $h_t's$ efficiently, but details are left for our future work.
This section applies the SVM model with leverage and several alternative models to three excess holding yields data. The descriptions of the data (labeled as TB, DGS and S&P500) are given below\footnote{The data are obtained from the website of Federal Reserve Bank of St. Louis.}.
The time series plots of three datasets are shown in Figures (ref), (ref) and (ref). Volatility clustering phenomena is observed for all three series, suggesting that the stochastic volatility models are appropriate to describe these excess holding yield data.
Using the same prior distributions as in illustrative examples in Section (ref), the proposed SVM models with leverage are fitted to TB, DGS and S&P500 data. We iterated MCMC simulation 50,000 times after discarding initial 10,000 MCMC draws as burn-in period using Algorithm 3 in Appendix (ref). The acceptance rates of the MH algorithms for $\alpha$ and $(\alpha, h)$ are 72.0% and 20.8% with TB data, 69.5% and 12.8% with DGS data, and 77.3% and 58.9% with S&P500 data, respectively.
Table (ref) shows posterior means, standard deviations, 95% credible intervals, inefficiency factors for parameters (IF), and the posterior probability that the parameter is positive for three datasets. Further, the IF's for log volatilities, $h_t$'s are found to be less than 80, which implies that our mixture sampler is highly efficient as shown in Section (ref). The coefficient $\beta$ is estimated to be greater than 0.6 for TB and DGS data, albeit close to zero for S&P500 data. In all cases, we find strong evidence of the positive risk premium since the posterior probability that $\beta$ is positive is almost one for TB and DGS data and 0.97 for S&P500 data. The autoregressive parameter $\phi$ for the log volatility process is estimated to be more than 0.9, suggesting the high persistence in the volatility as found in the various empirical studies in the previous literature. The correlation parameter $\rho$ is estimated to be negative for TB data and S&P500 data, implying a strong evidence of the leverage effect since the posterior probability that $\rho$ is negative is almost one, $Pr(\rho < 0|y) \approx 1.000$. On the other hand, there is no evidence that $\rho$ is negative for DGS data. Finally, Figures (ref), (ref) and (ref) show 95% credible intervals and posterior medians for $h$ in the case of TB, DGS and S&P500 data. Noting that
we further plotted moving average $\sum_{s=-10}^{10} z_{t+s}/21$ for reference where $z_t = \log (y_t^2) - E(\log \chi_1^2(\beta^2))$ is evaluated at the posterior means of $\beta$. The expected values of $\log \chi_1^2(\beta^2)$ are computed numerically as $-0.88, -0.78$, and $-1.27$ for TB, DGS, and S&P500 using Monte Carlo integration. The traceplot of the estimated log volatilities is similar to that of the moving average series taking account of 95% credible intervals. The large volatilities around years 1980, 2008 and 2020/3-2020/4 are well captured by the proposed model as shown in Figures (ref), (ref) and (ref) respectively.
In this section we use Bayesian marginal likelihoods to conduct a comparison of different stochastic volatility models. We calculate the marginal likelihood using the method of Chib(95). Suppressing the model index, this method is based on an identity introduced in that paper:
where the first term on the right side is the log of the likelihood. the second term is the prior, and the third is the posterior density. We evaluate each of these terms with the posterior mean of $\theta$. For each model, we calculate the first term using the particle filter method given in Section (ref), setting $I = 80,000$. To compute the posterior density ordinate, we apply ChibJeliazkov(01) to the MCMC draws from Algorithm 4 in Appendix (ref).
The SVML model had the highest log marginal likelihood for TB and S&P500 data, while the SVM model had the highest for DGS data. This implies our proposed model best describes the risk premium and the time varying volatility among competing models including the standard SV and SVL models. These results are also consistent with high posterior probabilities of $Pr(\beta > 0|y) $ for TB, DGS and S&P500 data, and of $Pr(\rho< 0|y)$ for TB and S&P500 data as given in Table (ref) of Section (ref).
In this paper, we have successfully extended the mixture sampler for the SV model to the SVM model of which the mean equation is described by the standard deviation of the error term as an independent variable. Our main point is the approximation of the distribution of $\log \chi_1^2(\beta^2)$ by mixture of normal distributions which is dependent on the parameter $\beta$. This approximation facilitates efficient sampling, leveraging well-established methods for the linear Gaussian state-space model. It is shown in simulation studies that our proposed method is implemented easily and works fast and efficiently. In the empirical studies of the excess holding yield data and S&P500 data, we conducted the model comparison among the SV and SVM models with and without leverage, and found that our proposed SVM models outperform other models in terms of the marginal likelihood. It shows that there exists the positive risk premiums and time-varying volatilities for all data, while the leverage effects are found to exist for TB and S&P500 data.
The authors thank the editor and anonymous referees for their helpful comments. This work is partially supported by JSPS KAKENHI [Grant number: 24H00142]. The computational results are obtained using Rcpp and Ox (see Doornik(07)).