# A Fast and Efficient Markov Chain Monte Carlo Method for Market Microstructure Model

Authors: Sun Yapeng; Peng Hui; Xie Wenbiao
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

The non-linear market microstructure (MM) model for financial time series modeling is a flexible stochastic volatility model with demand surplus and market liquidity. The estimation of the model is difficult, since the unobservable surplus demand is a time-varying stochastic variable in the return equation, and the market liquidity arises both in the mean term and in the variance term of the return equation in the MM model. A fast and efficient Markov Chain Monte Carlo (MCMC) approach based on an efficient simulation smoother algorithm and an acceptance-rejection Metropolis–Hastings algorithm is designed to estimate the non-linear MM model. Since the simulation smoother algorithm makes use of the band diagonal structure and positive definition of Hessian matrix of the logarithmic density, it can quickly draw the market liquidity. In addition, we discuss the MM model with Student-t heavy tail distribution that can be utilized to address the presence of outliers in typical financial time series. Using the presented modeling method to make analysis of daily income of the S&P 500 index through the point forecast and the density forecast, we find clear support for time-varying volatility, volatility feedback effect, market microstructure theory, and Student-t heavy tails in the financial time series. Through this method, one can use the estimated market liquidity and surplus demand which is much smoother than the strong stochastic return process to assist the transaction decision making in the financial market.


---

## Body

## 1. Introduction

The prevention and avoidance of financial risk has always been the core issue of investment theory and practice. Risk can be measured by the time-varying variance, i.e., volatility of return on investment. The main models to analyze time-varying volatility of financial time series are SV (stochastic volatility) type models (Taylor) [1] and GARCH (generalized autoregressive conditional heteroscedasticity) type models (Bollerslev) [2]. GARCH type models can describe some behavioral characteristics of volatility, but they regard volatility as a deterministic function of historical information, which makes them unable to flexibly fit financial time series. In theory, due to the introduction of stochastic volatility, SV type models are more in line with the characteristics of financial market and are then more suitable for the analysis and prediction of financial time series. In the fields of option pricing, risk management, and portfolio, SV type models have very important theoretical and practical significance in capturing the relation of asset return and volatility and thus can provide market information and decision-making reference for investors and market managers.The general SV models focus on modeling the dynamics of asset price itself and especially its volatility. Unlike those models, Bouchaud and Cont [3] proposed a phenomenological model called market microstructure (MM) model that describes the relationship between price and market liquidity as well as surplus demand. Based on the ideal of the MM model, Iino and Ozaki [4] proposed a continuous-time MM model that is a SV model. The estimation of the MM model is not straightforward, since the market liquidity appears both in the mean term and in the variance term of the model, and the surplus demand is also time-varying stochastic variable. In previous works, the non-linear Kalman filter and maximum likelihood method (see, e.g., [5–10]) were used to estimate the parameters of the MM model, but this kind of method has some shortcomings in dealing with strong non-Gaussian non-linear financial data. Thus, we consider using the MCMC method to improve estimation of the model.The MCMC method was promoted by Gelfand and Smith [11], and it is equitable to say that the MCMC method has revitalized (perhaps even revolutionized) Bayesian statistics. Bardsley and Cui [12] developed a scalable optimization-based MCMC method for solving hierarchical Bayesian inverse problems with non-linear parameter-to-observable maps and a broader class of hyperparameters. Zhou and Tartakovsky [13] proposed an adaptive MCMC method, in which a transmission model is replaced by a fast and accurate substitutive model that adopts a deep convolutional neural network to reduce calculation burden. For the shortcoming of MCMC, i.e., the huge amount of computation when the dataset is too large, Nemeth and Fearnhead [14] studied the latest development of scalable Monte Carlo algorithm for reducing computational cost of MCMC. In this paper, we focus on the stochastic gradient MCMC (SGMCMC), which uses the data subsampling technique to reduce the cost of each iteration of MCMC and to achieve good results.In this paper, unlike the non-linear Kalman filter and maximum likelihood method-based estimation way for the MM model (see, e.g., Peng et al. [5]; Peng et al. [6]; Peng et al. [7]; and Qin et al. [9]), we develop an efficient MCMC method to estimate the market microstructure model. Note that the Kalman filter-based simulation smoother algorithm (see, e.g., De Jong [15]; Koopman [16]; De Jong and Shephard [17]; and Durbin and Koopman [18]) and the auxiliary mixture sampler (see, e.g., Kim et al. [19]; Omori et al. [20]; Nakajima and Omori [21]; Delatola and Griffin [22]; and Jensen and Maheu [23]) cannot be directly applied to the MM model, since the market liquidity enters the mean term of return equation. To make an efficient posterior inference easier using WinBUGS (Windows version of Bayesian Analysis Using Gibbs Sampler, a statistical software for Bayesian analysis using MCMC methods), Xi et al. [8] and Xi et al. [10] developed a MCMC method to estimate a MM model. But the status of the chain after a huge number of steps is then used as a sample of the desired distribution, which may lead to complex desired distribution and an inaccurate estimate result. Thus, the modification of the simulation smoother algorithm is required for the estimation of the MM model.We extend the precision-based simulation smoother algorithm in Chan and Jeliazkov [24] and McCausland et al. [25] to the MM model and other SV models, which was originally designed for linear Gaussian state-space models. The precision-based simulation smoother includes drawing state variables or disturbances in discrete-time state-space model from their conditional distribution given parameters and observations, which can be approximate to the conditional distribution of the state in a non-Gaussian or non-linear model by its counterpart in a linear Gaussian model. The proposed MCMC method in this paper is based on the precision-based simulation smoother algorithms instead of the conventional Kalman filter. The literatures using the former methods include Rue [26] for linear Gaussian Markov random fields, Chan and Jeliazkov [24] and McCausland et al. [25] for linear Gaussian state-space models, Rue et al. [27] for non-linear Markov random fields, and McCausland [28] and Djegnéné and Mccausland [29] for non-linear state space models. McCausland et al. [25] provided a detailed comparison between the Kalman filter-based algorithm and precision-based algorithm and showed that the latter is mainly more efficient and faster, since it draws state variables based on band diagonal structure and the Hessian matrix of the logarithmic density is positive definite. Our presented method inherits and sublimes the characteristics of the original precision-based method.In contrast to the earlier works [3, 30–34], we show that the proposed model in this paper can capture the volatility feedback effect and the effect of the dynamic offset between supply and demand. The proposed method directly simulates smoother log-volatility from its posterior distribution using a proposal distribution, which approximates the target distribution and is different from the auxiliary mixture sampler [19] in which the return equation is first transformed into a linear equation. Our method takes all the market liquidity in the MM model as a single MCMC block, which is different from the single-move sampler [20–24], so it results in efficiency improvement when market liquidity is posterior serial dependence. Unlike the Kalman filter-based simulation smoother [15–18], our method achieves efficiency gains by exploiting the special structure of the problem, particularly in which the Hessian of the log-conditional density of the market liquidity given the return and the surplus demand is a diagonal band and positive matrix. Unlike the Hessian method approximation [28, 29] for the log-likelihood with respect to the log-volatility with the first five derivatives, our method obtains a Gaussian approximation with only the first two derivatives and correct the approximation error with an acceptance-rejection Metropolis–Hastings (ARMH) step [35–38].The remainder of this paper is structured as follows. In Section2, we introduce six SV type models including two newly proposed MM models. In Section 3, we present the estimation approach to the market liquidity in the MM model. To this end, we propose a precision-based simulation smoother to obtain a Gaussian approximation of the market liquidity density, and the acceptance-rejection Metropolis–Hastings (ARMH) step is used to improve efficiency of the posterior simulator. In Section 4, using the method proposed in Section 3, we present the efficient posterior simulators for other states and parameters of the SV type models. In Section 5, we show the simulation study result for SV model and MM model using the proposed estimation method. In Section 6, we give the modeling method application results for real data using the six SV type models that are estimated by the proposed estimation methods. Some conclusions are given in Section 7.


## 2. Preliminary: Problem Description

This paper studies the financial time series modeling issues using stochastic process models. The proposed estimation method for the stochastic models may be also used to the stochastic signal processing problems in other fields (see, e.g., Valipour [39]; Valipour et al. [40]; and Valipour [41]). Financial time series analysis is mainly to analyze statistical properties of return. In finance, one-period simple return from date t−1 to date t is Rt=pt−pt−1/pt−1, where pt is the asset price. The one-period simple compound return is 1+Rt=pt/pt−1. For convenience, the log return rt=logpt/pt−1=logpt−logpt−1 is usually used in modeling; especially in case study, one often uses rt=100logpt/pt−1 to represent the return for clearly showing the number.Since the avant-garde works of Box and Jenkins, the autoregressive (AR) model and autoregressive moving average (ARMA) models have become typical tools for modeling and detecting time series (see, e.g., Valipour et al. [42]). Although the theoretical proof and empirical results of such models have been verified, a large number of literatures emphasized the importance of introducing time-varying volatility into financial time series for estimation and forecasting. In finance, volatility refers to the degree of variation of transaction price series over time, which is measured by the standard deviation of return. Standard AR and ARMA models with constant variance (or constant volatility) are seemingly not flexible enough for financial time series modeling. Because the stochastic volatility (SV) model (Taylor [1]) supports time-varying stochastic volatility (exponential variance), the SV model is widely used in financial data modeling. The general form of the SV model is to treat the logarithmic volatility as a first-order autoregressive process, as shown in the following equation:(1)rt=μ+eσt/2εtr,(2)σt=γ+δσt−1−γ+ηεtσ,where εtr and εtσ are mutually and serially uncorrelated Gaussian noise with zero mean and unit variance, μ is the mean of return, γ is the mean of long-term volatility in (2), δ is the rate where the volatility reverts toward its long-term mean, and η is the volatility during the volatility process in (2). In the model, δ<1 and the logarithmic volatility is initialized with NIDγ,η2/1−δ2, where NID⋅,⋅ denotes the normal and independent Gaussian distribution.The relationship between the expected stock index returns and conditional volatility has attracted much attention in financial economic literature, in which the volatility feedback effect is an important property. Financial volatility feedback is the effect of volatility to return in financial market. Schwert et al. [30] and Campbell and Hentschel [31] gave the explanation on the volatility feedback effect as follows. If volatility is priced, an expected increase in volatility would rise at the requisite rate of return; conversely, an immediate stock-price decline is needed in order to get higher future returns. Therefore, the potential causality of the volatility feedback effect runs from volatility to prices. A positive relation between expected returns and expected volatility conforms to the capital asset pricing model (CAPM), since rational risk-averse investors need higher expected returns during more variable periods (see, e.g., Koopman and Uspensky [32]). Unfortunately, the SV model cannot capture the financial volatility feedback effect and measure the relationship between expected return and expected volatility. To solve the problem, Koopman and Uspensky [32] proposed the stochastic volatility-in-mean (SVM) model, including the unobservable volatility square as an explanatory variable in the mean term of the return equation as follows:(3)rt=ϕeσt+eσt/2εtr,(4)σt=γ+δσt−1−γ+ηεtσ,where the regression coefficient ϕ measures the volatility-in-mean effect or volatility feedback effect; when ϕ=0, the SVM model becomes SV models (1) and (2). The SVM model builds a strong ex ante relationship between expected return and expected volatility, and the volatility feedback effect could be interpreted as a stronger ex ante relationship between the expected returns and expected volatility.To extract more information from financial data, in this paper, we study the market microstructure model proposed by Bouchaud and Cont [3], which is described by the following equations:(5)dPt=λϕtdt,(6)ϕt=ϕt+−ϕt−,where Pt=logpt is the log price of an asset, ϕt is the market surplus demand, and λ is the market liquidity. At any given instant of time, there is an instantaneous demand ϕt+ and an instantaneous supply ϕt− for the asset. The dynamical equations (5) and (6) describe the impact of supply-demand balance, i.e., ϕt tends to push the price up (if ϕt+>ϕt−) or down (if ϕt+<ϕt−) in the other case, so it is named as surplus demand. λ is a measurement of market depth, i.e., the surplus demand is required to remove the price by a unit. When inverse of λ is high, the market is able to “absorb” supply/demand offsets by the little price changes, so it is named as market liquidity, which is a market’s capability to buy or sell an asset without making big change in the asset’s price (see, e.g., Bouchaud and Cont [3]). The dynamic equation (5) assumes that the price Pt is driven by the surplus demand ϕt, and the range of price changes depends on the market liquidity 1/λ.However, in the original setting, models (5) and (6) only provide an abstract description to the dynamics of stock market because the surplus demand ϕt and the market liquidity 1/λ cannot be easily measured directly. Iino and Ozaki [4] treated the market liquidity 1/λt and the surplus demand ϕt as stochastic processes and proposed a market microstructure model, which consists of a set of continuous time stochastic differential equations. Peng et al. (see e.g., Peng et al. [5]; Peng et al. [6]; Peng et al. [7]; Xi et al. [8]; Qin et al. [9]; Xi et al. [10]), also extended the market microstructure model in continuous-time domain and discrete-time domain. More generally, in this paper, we build the continuous-time market microstructure model as follows:(7)dPt=λtϕtdt+λtdwP,t,dϕt=α+βϕt−αdt+κdwϕ,t,dlogλt=γ+δlogλt−γdt+ηdwλ,t,where wP,t, wϕ,t, and wλ,t are the independent Wiener processes with unit variance and α, γ, β, δ, κ, and η are the constant parameters. In the related literatures, a continuous model is often discretized to make estimation easy. In mathematics, the Euler–Maruyama method is a simple and common method for solving approximate numerical solutions of stochastic differential equations. Using the Euler–Maruyama discrete-time approximation method, we can build a discrete-time market microstructure (MM) model as follows:(8)rt=ϕteσt+eσt/2εtr,(9)ϕt=α+βϕt−1−α+κεtϕ,(10)σt=γ+δσt−1−γ+ηεtσ,where σt≡logλt, rt≡Pt−Pt−1, and εtr, εtϕ, and εtσ are mutually and serially uncorrelated Gaussian noises with zero mean and unit variance. Equation (8) characterizes the financial asset price process. It is obvious that the conditional expected mean value and the conditional variance of the financial asset price process in (8) are given by ϕteσt and eσt, respectively. Equations (9) and (10) describe the unobservable surplus demand ϕt and market liquidity σt as an autoregressive process, respectively. The MM models (10)–(12) are similar to the SVM model, in which the unobservable logarithmic volatility or market liquidity is described by an autoregressive process.In contrast to the SV model, MM models (10)–(12) place a strong ex ante relation between expected return and expected volatility and can capture the volatility feedback effect. Different from the general SVM model, MM models (10)–(12) place a time-varying surplus demand ϕt instead of a constant coefficient ϕ used in SVM model, and thus MM models (10)–(12) can capture the effect of the dynamic offset between supply and demand. So, MM models (10)–(12) are more flexible stochastic volatility models, and one may utilizes the surplus demand process that is a more stable and predictable process, rather than the highly random return process, to assist decision making in financial investment process. It is therefore suggested that one may allocate assets by monitoring and making use of the surplus demand information. This methodology was implemented in, e.g., Peng et al. [5], Peng et al. [6], Peng et al. [7], and Qin et al. [9].Besides stochastic volatility, another prominent characteristic in typical financial time series is the existence of outliers. The stochastic volatility model with Student-t thick tail distribution (SVt) is a common method to study this characteristic. Evidence in favor of fat tails in financial data has been uncovered by Geweke [33], Gallant et al. [34], and Jacquier et al. [43]. Since the Student-t distribution can be written as a scale mixture of Gaussian distributions (see, e.g., Geweke [44]), the potential variable representation of the SVt model is as follows:(11)rt=μ+eσt/2λt1/2εtr,(12)σt=γ+δσt−1−γ+ηεtσ,(13)λt|υ∼IGυ2,υ2,where λ=λ1,λ2,…,λTT are conditionally independent, given the model parameters (including υ) and data, and IG⋅,⋅ denotes the inverse-gamma distribution. Because the Student-t distribution has a heavier tail than the Gaussian distribution, outliers are allowed to occur more frequently in the SVt model than in the standard SV model. In particular, by utilizing the Student-t heavy tails, we can generalize the original specification in the SVM model and MM model and obtain their counterparts, i.e., the SVM plus Student-t heavy tail (SVMt) model:(14)rt=ϕeσt+eσt/2λt1/2εtr,(15)σt=γ+δσt−1−γ+ηεtσ,(16)λt|υ∼IGυ2,υ2,and the MM plus Student-t heavy tail (MMt) model:(17)rt=ϕteσt+eσt/2λt1/2εtr,(18)ϕt=α+βϕt−1−α+κεtϕ,(19)σt=γ+δσt−1−γ+ηεtσ,(20)λt|υ∼IGυ2,υ2.All six models studied in this paper are summarized in Table1. It is worth noting that the structures of MM models (10)–(12) and MMt models (19)–(22) are different from those of MM and MMt models previously studied [5, 8–10]. In the return equation of the proposed MM type models in this paper, the mean term is a stochastic variable. However, in the previous works [5, 8–10], the mean term is built by a deterministic function of historical information because the MM model is estimated by a Kalman filter and maximum likelihood method-based approach, which cannot cope with the stochastic mean term. Therefore, the proposed MM models in this paper can better describe the dynamic behavior between return and market liquidity as well as surplus demand of financial market.Table 1
List of different stochastic volatility (SV) models.
ModelDescriptionSV (1) and (2)Standard SV models in whichσt follows a stationary AR(1)SVt (13)–(15)SV models whose return innovation follows a Student-t distributionSVM (3) and (4)SV models butσt enters the mean term of return equation as a covariateSVMt (16)–(18)Same as SVM but the return innovation follows a Student-t distributionMM (10)–(12)Similar to SVM butϕt follows a stationary AR(1)MMt (19)–(22)Same as MM but the return innovation follows a Student-t distributionFor the estimation of the SVM model, Koopman and Hol Uspensky [32] developed a simulated maximum likelihood estimator built on the Kalman filter-based simulation smoother algorithm and approximated the Gaussian model to fit the SVM model. Since there are only one state and a constant ϕ in the SVM model, it is able to maximize the likelihood numerically to get the maximum likelihood estimate of σt. However, this approach cannot be easily generalized to the time-varying surplus demand ϕt in MM models (10)–(12) and MMt models (19)–(22) due to the likelihood evaluation in the case that involves two kinds of status, ϕt and σt. Meanwhile, applying the Kalman filter-based smoother to obtain a smoothed estimate for σt leads to an iterative process, which needs a large computation burden. More details were given by Durbin and Koopman [45].In this paper, we present a fast and efficient MCMC approach based on an efficient simulation smoother algorithm and an acceptance-rejection Metropolis–Hastings algorithm to estimate the models listed in Table1. First, we directly simulate smoother log-volatility from its posterior distribution by using the proposal distribution, which approximates the target distribution. It is unlike Kalman filter-based simulation smoother and the auxiliary mixture sampler, in which the return equation is firstly transformed into a linear equation, and any non-Gaussian distributions in this transformed model may be close to Gaussian mixtures. However, through modeling the distribution of rt∗=logrt2, rather than rt directly, information about the sign of the return rt is lost. Therefore, the distribution of rt or εtr can only be recovered by assuming the distribution of the sign of the return (see, e.g., Kim et al. [19]; Omori et al. [20]; Nakajima and Omori [21]; Delatola and Griffin [22]; and Jensen and Maheu [23]). Unlike the above approach, the proposed precision-based simulation smoother has a Gaussian approximating model for the conditional return distribution, and it preserves the parametric specification of the log-volatility process.Second, we take all the market liquidityσ=σ1,σ2,…,σTT as a single MCMC block. It brings efficiency improvement when market liquidity is posterior serial dependence. This is achieved by constructing a Gaussian approximation of the target distribution and correcting the approximation error using the acceptance-rejection Metropolis–Hastings step. This Gaussian approximation is closer to the target distribution than that obtained by any multivariate Gaussian approximation approach that updates the state vector but usually only for about 10–50 observations at a time (see, e.g., [35–37]).Third, the new method can improve efficiency by exploiting the special structure of the problem, particularly in which the Hessian of the log-conditional density of the market liquidityσt given the return rt and the excess demand ϕt is a diagonal band and positive matrix. This feature turns out to be crucial in developing a simple, efficient, and fast precision-based simulation smoother for the market liquidity σt. Using forward- and back-substitution to simulate smoother market liquidity σt, the Cholesky decomposition diagonal band matrix Hessian can be obtained in OT operations instead of OT3 operations for full matrices. Meanwhile, applying the precision-based smoother to obtain a smoothed estimate for market liquidity σt leads to a very simple and fast Newton–Raphson iterative procedure, since the Hessian of the log-conditional density of the market liquidity σt is a positive matrix. These two features are essential for reducing the computational costs.The last and most important point is the Gaussian approximation of the conditional density of logarithmic volatility. The HESSIAN method of McCausland [28] and Djegnéné and Mccausland [29] provided an excellent approximation for the log-likelihood with respect to the log-volatility with the first five derivatives. However, it also needs a substantial burden on the end user, which refers to the first five derivatives, so it is not easily extended to complicated volatility-in-mean structure and the multivariable MM model. We present a simple, efficient, and fast precision-based simulation smoother algorithm to obtain a Gaussian approximation for the conditional density of the log-volatility with only the first two derivatives and to correct the approximation error with an acceptance-rejection Metropolis–Hastings (ARMH) step (see, e.g., Tierney [38]; Chib and Jeliazkov [46]; Dimitrakopoulos and Kolossiatis [47]; and Blasques et al. [48]).


## 3. Efficient Posterior Simulator for the Market Liquidity in MMt Model

In MMt models (19)–(22), since the market liquidity σt appears in both mean term and variance term of return equation (17) as volatility and the model is non-linear in the market liquidity σt, one cannot directly sample the market liquidity σt from the complicate joint conditional density. In this part, we first propose a precision-based simulation smoother to obtain a Gaussian approximation of the market liquidity density for the non-linear non-Gaussian MMt model. Then, the acceptance-rejection Metropolis–Hastings (ARMH) step is discussed, which may significantly improve the posterior simulator efficiency, since it corrects the approximation error. The efficient posterior simulators for other parameters and other states are presented in Section 4.


### 3.1. Gaussian Approximation of Joint Conditional Density

There are usually two approaches that are commonly used for estimating non-linear non-Gaussian state-space models. The first is the so-called sequential Monte Carlo method or more popularly known as particle filter [49] that involves sequential importance sampling and bootstrap resampling. Despite recent advances, particle filter is still quite computationally intensive and is not designed for efficient simulation smoother for the states. In fact, in posterior simulation with a particle filter, it is a common practice to generate candidate draws on parameters via a random walk sampler and then use a particle filter to compute the acceptance probability of the candidate draw. More recently, there has been work on using a particle filter to obtain candidate draws for the states (see e.g., Andrieu et al. [50]). This direction looks promising, but most applications using these methods to date are limited to univariate state-space models due to computational limitations.The alternative approach to estimating general state space models is built on fast approximations of the conditional state density, where the approximating density is used for posterior simulation via the independence-chain Metropolis–Hastings (MH) algorithm. Of course, the main challenges are just crucial to have a fast routine to obtain a good approximation for the density, and it should be easy to generate candidate draws from the approximating density. Durbin and Koopman [45] and Shephard and Pitt [51] considered approximating the log target density around the mode by a second-order Taylor expansion. This approximation provides a Gaussian density, where its mean is the mode of the target density and its precision or inverse variance equals the negative Hessian evaluated at the mode. Candidate draws for the states are then generated via the Kalman filter-based simulation smoother. However, one problem with sampling states in a single block with the Gaussian proposal density is that the acceptance rate in the Metropolis–Hastings (MH) step can be quite low, at least in the context of stochastic volatility models. It is therefore suggested that the states are divided into blocks, and each block is sampled sequentially via the MH step. This methodology was implemented in, e.g., Watanabe and Omori [35], Omori and Watanabe [36], and Abanto-Valle et al. [37]. In addition, the research of Chan [52], Chan and Strachan [53], Zhang [54], and Zhang [55] also brings inspiration to this paper.Departing from the obvious Gaussian approximation, McCausland [28] and Djegnéné and Mccausland [29] introduced the HESSIAN method that provides an excellent approximation, and it results in a highly efficient simulation smoothing sampling algorithm. However, the HESSIAN method requires the first five derivatives of the log-likelihood with respect to the states, which gives a substantial burden on the end user. A more severe restriction is that currently it can only be applied to univariate state-space models.In contrast to the HESSIAN method, we present a more simple precision-based simulation smoother algorithm to obtain a Gaussian approximation for the conditional density of the log-volatility with only the first two derivatives. By exploiting the diagonal band and positive matrix structure of the precision matrix, the precision-based algorithm is more efficient than Kalman filter-based methods in general. This feature is crucial since one needs to obtain the approximating density tens of thousands times in a full Bayesian analysis via MCMC. More importantly, the marginal cost of obtaining additional draws under the precision-based algorithm is much smaller compared to Kalman filter-based methods. We exploit this important feature for two purposes: efficient simulation of the states and evaluation of the integrated likelihood. In addition, different from being divided into blocks, we sample the market liquidity in a single block and develop an acceptance-rejection Metropolis–Hastings (ARMH) algorithm to correct the approximation error.


### 3.2. Precision-Based Simulation Smoother of the Market Liquidity

In this section, a precision-based simulation smoother is built, and a simple, efficient, and fast precision-based smoother and Newton–Raphson iterative procedure are presented to get a Gaussian approximation density for the joint conditional density of the states. For notational convenience, letr=r1,r2,…,rTT, ϕ=ϕ1,ϕ2,…,ϕTT, σ=σ1,σ2,…,σTT, and λ=λ1,λ2,…,λTT.Now, we discuss a quick way to get the Gaussian approximation density for the joint conditional densitypσ|r,ϕ,λ,θ where θ=γ,δ,η2. By Bayes’ theorem, one has pσ|r,ϕ,λ,θ∝pr|σ,ϕ,λpσ|θ, and the prior density pσ|θ that characterizes the state equation (19) is Gaussian density. By approximating the observation density pr|σ,ϕ,λ, which characterizes return equation (17) via a Gaussian density in σ, one obtains a Gaussian approximation of pσ|r,ϕ,λ,θ at once. To this end, following equation (17), the log-density logprt|σt,ϕt,λt related to σt is given as follows:(21)logprt|σt,ϕt,λt=−σt2−rt−ϕteσt22λteσt+c1,where c1 is independent of σt. It is easy to examine that the first and second partial derivatives of the log-density with respect to σt are given by(22)gt=∂logprt|σt,ϕt,λt∂σt=−12+rt2e−σt−ϕt2eσt2λt,ht=−∂2logprt|σt,ϕt,λt∂σt2=−∂gt∂σt=rt2e−σt+ϕt2eσt2λt.As the first and second derivatives can be calculated quickly, we should note that when a pointσ^t is given, one can estimate the log conditional likelihood logprt|σt,ϕt,λt by making use of a second-order Taylor expansion around σ^t to obtain(23)logqrt|σt,ϕt,λt=logprt|σ^t,ϕt,λt+σt−σ^tgt−σt−σ^t2ht2.Letlogqr|σ,ϕ,λ=∑t=1Tlogqrt|σt,ϕt,λt, which can be written in matrix form as follows:(24)logqr|σ,ϕ,λ=logpr|σ^,ϕ,λ+σ−σ^TG−σ−σ^THσ−σ^2=−σTHσ2+σTG+Hσ^+c2,where c2 is independent of σ, G=g1,g2,…,gTT, H=diagh1,h2,…,hT, and σ^=σ^1,σ^2,…,σ^TT. Next, rewrite state equation (19) in matrix form below:(25)Aσσ=μσ+εσ,εσ∼NID0,Λσ,where μσ=γ,1−δγ,…,1−δγT, Λσ=diagη2/1−δ2,η2,…,η2, and(26)Aσ=100…0−δ10…00−δ1…0⋮⋱⋱⋱⋮00…−δ1.By a change of variable, one has(27)pσ|θ∼NIDAσ−1μσ,AσTΛσ−1Aσ−1,with log-density(28)logpσ|θ=−σ−Aσ−1μσTAσTΛσ−1Aσσ−Aσ−1μσ2+c3=−σTAσTΛσ−1Aσσ2+σTAσTΛσ−1μσ+c4,where c3 and c4 are independent of σ. Then, considering the expressions in (24) and (28), one has(29)logqσ|r,ϕ,λ,θ=logqr|σ,ϕ,λ+logpσ|θ=−σTH+AσTΛσ−1Aσσ2+σTG+Hσ^+AσTΛσ−1μσ+c5,where c5 is independent of σ. Since the expression in (29) is the log-kernel of NIDμ^σ,Λ^σ in σ, where Λ^σ−1=H+AσTΛσ−1Aσ is the precision matrix or inverse of the variance and μ^σ=Λ^σG+Hσ^+AσTΛσ−1μσ is the mean vector, pσ|r,ϕ,λ,θ can be calculated by the Gaussian density(30)qσ|r,ϕ,λ,θ∼NIDμ^σ,Λ^σ.This Gaussian approximation density is then put into use as the proposal density in the acceptance-rejection Metropolis–Hastings (ARMH) step.Note that using forward- and back-substitution to simulate smoother market liquidityσ, one can calculate the (banded) Cholesky decomposition Cσ of Λ^σ such as CσCσT=Λ^σ and obtain a random draw from NIDμ^σ,Λ^σ efficiently such that σ=μ^σ+Cσu and u∼NID0,IT. A crucial result is that the Cholesky decomposition Cσ can be obtained in OT operations rather than OT3 operations for full matrices, since Λ^σ is a band diagonal matrix (see, e.g., Chan and Jeliazkov [24]).To this end, it remains to select the pointσ^ for the Taylor expansion in (21), which is chosen to be the mode of pσ|r,ϕ,λ,θ and can be quickly obtained by the numerical optimization. Using the Newton–Raphson method, one can quickly locate the maximum of qσ|r,ϕ,λ,θ, or equivalently that of logqσ|r,ϕ,λ,θ, or to select the point σ^ as the mode of pσ|r,ϕ,λ,θ or smooth σ. It follows from (29) that negative Hessian of logqσ|r,ϕ,λ,θ evaluated at σ=σ^ is Λ^σ−1 and the gradient at σ=σ^ is −Λ^σ−1σ^+Λ^σ−1μ^σ. Hence, one can carry out the Newton–Raphson method as follows: initializing with σ=σ1 for some lasting vector σ1. For i=1,2,⋯, assume σ^=σi in the evaluation of Λ^σi−1 and −Λ^σi−1σi+Λ^σi−1μ^σi and compute(31)σi+1=σi+Λ^σi−Λ^σi−1σi+Λ^σi−1μ^σi=μ^σi.Note that sinceqσ|r,ϕ,λ,θ is Gaussian density, its mode σ^ equals its mean vector μ^σ. Repeating this procedure until some convergence criterion is reached, e.g., we get the mode σ^=σi+1 or μ^σ=μ^σi when σi+1−σi<ε for some prefixed tolerance level ε>0. Also, note that the negative Hessian Λ^σ−1 is positive definite for all σ. This guarantees a fast convergence of the Newton–Raphson method.


### 3.3. Acceptance-Rejection Metropolis–Hastings (ARMH) Algorithm

As mentioned before, one problem with sampling states in a single block with the Gaussian proposal density is that the accepting rate in the MH step can be very low, and it is probably because that the Gaussian approximation density is not sufficiently precise. In this section, different from being divided into blocks, implemented in, e.g., Watanabe and Omori [35], Omori and Watanabe [36], and Abanto-Valle et al. [37], we sample the market liquidity σ in a single block and exploit an acceptance-rejection Metropolis–Hastings (ARMH) algorithm to correct the approximation error.As its name implies, the ARMH algorithm (Tierney [38]) is a MCMC sampling procedure that links classic accept-reject sampling with the Metropolis–Hastings algorithm. In the setting, the target density is the joint conditional density pσ|r,ϕ,λ,θ∝pr|σ,ϕ,λpσ|θ of σ and we have the Gaussian approximation proposal density qσ|r,ϕ,λ,θ from which we generate candidate draws σ. In the classic accepting-rejecting sampling, an important requirement is that there exists a constant c, which is selected by equation (34), such that(32)pσ|r,ϕ,λ,θ<cqσ|r,ϕ,λ,θ,for σ in support of pσ|r,ϕ,λ,θ. The ARMH releases the domination condition (32) such that when it is not satisfied for σ, one may resort to the MH algorithm. To show Algorithm 1, it is useful to define the set(33)ℂ=σ:pσ|r,ϕ,λ,θ<cqσ|r,ϕ,λ,θ,and let ℂc denote the complement code. Then, the ARMH procedure is as follows.Algorithm 1: Acceptance-rejection Metropolis–Hastings with Gaussian proposal density [8, 24, 38, 56].
(1)
Accept-reject step: generate a drawσ∗∼qσ|r,ϕ,λ,θ∼NIDμ^σ,Λ^σ and accept σ∗ with probability min1,pσ∗|r,ϕ,λ,θ/cqσ∗|r,ϕ,λ,θContinue the above process until a drawσ∗ is accepted.(2)
Metropolis–Hastings step: given the currentσ and the proposed draw σ∗(a) ifσ∈ℂ, accept σ∗ with probability 1;(b) ifσ∈ℂc and σ∗∈ℂ, accept σ∗ with probabilitymin1,cqσ|r,ϕ,λ,θ/pσ|r,ϕ,λ,θ(c) ifσ∈ℂc and σ∗∈ℂc, accept σ∗ with probabilitymin1,qσ|r,ϕ,λ,θpσ∗|r,ϕ,λ,θ/qσ∗|r,ϕ,λ,θpσ|r,ϕ,λ,θotherwise returnσ.Adjusting the original Gaussian proposal densityqσ|r,ϕ,λ,θ by the accept-reject step, a better approximation of the target density pσ|r,ϕ,λ,θ is achieved, i.e., the new proposed density is identical with the target density on the set ℂ (albeit different normalizing constants c), whereas on ℂc, the new proposed density is decreased to the original one.Chib and Jeliazkov [46] presented a practical way to select the constant c and the trade-off in such a choice, which is outlined here. Notice that if a bigger c is chosen, the set ℂ is larger and one is more likely to accept the candidate σ∗. On the other hand, selecting a smaller c will reduce the sampling efficiency, so it will generate more rejections during the sampling steps of accepting-rejecting, resulting in increased invalid and large calculations. A practical way to strike a balance between these two conflicting considerations is to set(34)c=bpσ^|r,ϕ,λ,θqσ^|r,ϕ,λ,θ,where σ^ is the mode of the conditional density qσ|r,ϕ,λ,θ and b is, say, between 1 and 5. Such kind of choice would make sure that c is sufficiently small to decrease the required number of draws from qσ|r,ϕ,λ,θ, while it is also big enough so that the set ℂ contains the mode σ^ and its neighboring points.


## 3.1. Gaussian Approximation of Joint Conditional Density

There are usually two approaches that are commonly used for estimating non-linear non-Gaussian state-space models. The first is the so-called sequential Monte Carlo method or more popularly known as particle filter [49] that involves sequential importance sampling and bootstrap resampling. Despite recent advances, particle filter is still quite computationally intensive and is not designed for efficient simulation smoother for the states. In fact, in posterior simulation with a particle filter, it is a common practice to generate candidate draws on parameters via a random walk sampler and then use a particle filter to compute the acceptance probability of the candidate draw. More recently, there has been work on using a particle filter to obtain candidate draws for the states (see e.g., Andrieu et al. [50]). This direction looks promising, but most applications using these methods to date are limited to univariate state-space models due to computational limitations.The alternative approach to estimating general state space models is built on fast approximations of the conditional state density, where the approximating density is used for posterior simulation via the independence-chain Metropolis–Hastings (MH) algorithm. Of course, the main challenges are just crucial to have a fast routine to obtain a good approximation for the density, and it should be easy to generate candidate draws from the approximating density. Durbin and Koopman [45] and Shephard and Pitt [51] considered approximating the log target density around the mode by a second-order Taylor expansion. This approximation provides a Gaussian density, where its mean is the mode of the target density and its precision or inverse variance equals the negative Hessian evaluated at the mode. Candidate draws for the states are then generated via the Kalman filter-based simulation smoother. However, one problem with sampling states in a single block with the Gaussian proposal density is that the acceptance rate in the Metropolis–Hastings (MH) step can be quite low, at least in the context of stochastic volatility models. It is therefore suggested that the states are divided into blocks, and each block is sampled sequentially via the MH step. This methodology was implemented in, e.g., Watanabe and Omori [35], Omori and Watanabe [36], and Abanto-Valle et al. [37]. In addition, the research of Chan [52], Chan and Strachan [53], Zhang [54], and Zhang [55] also brings inspiration to this paper.Departing from the obvious Gaussian approximation, McCausland [28] and Djegnéné and Mccausland [29] introduced the HESSIAN method that provides an excellent approximation, and it results in a highly efficient simulation smoothing sampling algorithm. However, the HESSIAN method requires the first five derivatives of the log-likelihood with respect to the states, which gives a substantial burden on the end user. A more severe restriction is that currently it can only be applied to univariate state-space models.In contrast to the HESSIAN method, we present a more simple precision-based simulation smoother algorithm to obtain a Gaussian approximation for the conditional density of the log-volatility with only the first two derivatives. By exploiting the diagonal band and positive matrix structure of the precision matrix, the precision-based algorithm is more efficient than Kalman filter-based methods in general. This feature is crucial since one needs to obtain the approximating density tens of thousands times in a full Bayesian analysis via MCMC. More importantly, the marginal cost of obtaining additional draws under the precision-based algorithm is much smaller compared to Kalman filter-based methods. We exploit this important feature for two purposes: efficient simulation of the states and evaluation of the integrated likelihood. In addition, different from being divided into blocks, we sample the market liquidity in a single block and develop an acceptance-rejection Metropolis–Hastings (ARMH) algorithm to correct the approximation error.


## 3.2. Precision-Based Simulation Smoother of the Market Liquidity

In this section, a precision-based simulation smoother is built, and a simple, efficient, and fast precision-based smoother and Newton–Raphson iterative procedure are presented to get a Gaussian approximation density for the joint conditional density of the states. For notational convenience, letr=r1,r2,…,rTT, ϕ=ϕ1,ϕ2,…,ϕTT, σ=σ1,σ2,…,σTT, and λ=λ1,λ2,…,λTT.Now, we discuss a quick way to get the Gaussian approximation density for the joint conditional densitypσ|r,ϕ,λ,θ where θ=γ,δ,η2. By Bayes’ theorem, one has pσ|r,ϕ,λ,θ∝pr|σ,ϕ,λpσ|θ, and the prior density pσ|θ that characterizes the state equation (19) is Gaussian density. By approximating the observation density pr|σ,ϕ,λ, which characterizes return equation (17) via a Gaussian density in σ, one obtains a Gaussian approximation of pσ|r,ϕ,λ,θ at once. To this end, following equation (17), the log-density logprt|σt,ϕt,λt related to σt is given as follows:(21)logprt|σt,ϕt,λt=−σt2−rt−ϕteσt22λteσt+c1,where c1 is independent of σt. It is easy to examine that the first and second partial derivatives of the log-density with respect to σt are given by(22)gt=∂logprt|σt,ϕt,λt∂σt=−12+rt2e−σt−ϕt2eσt2λt,ht=−∂2logprt|σt,ϕt,λt∂σt2=−∂gt∂σt=rt2e−σt+ϕt2eσt2λt.As the first and second derivatives can be calculated quickly, we should note that when a pointσ^t is given, one can estimate the log conditional likelihood logprt|σt,ϕt,λt by making use of a second-order Taylor expansion around σ^t to obtain(23)logqrt|σt,ϕt,λt=logprt|σ^t,ϕt,λt+σt−σ^tgt−σt−σ^t2ht2.Letlogqr|σ,ϕ,λ=∑t=1Tlogqrt|σt,ϕt,λt, which can be written in matrix form as follows:(24)logqr|σ,ϕ,λ=logpr|σ^,ϕ,λ+σ−σ^TG−σ−σ^THσ−σ^2=−σTHσ2+σTG+Hσ^+c2,where c2 is independent of σ, G=g1,g2,…,gTT, H=diagh1,h2,…,hT, and σ^=σ^1,σ^2,…,σ^TT. Next, rewrite state equation (19) in matrix form below:(25)Aσσ=μσ+εσ,εσ∼NID0,Λσ,where μσ=γ,1−δγ,…,1−δγT, Λσ=diagη2/1−δ2,η2,…,η2, and(26)Aσ=100…0−δ10…00−δ1…0⋮⋱⋱⋱⋮00…−δ1.By a change of variable, one has(27)pσ|θ∼NIDAσ−1μσ,AσTΛσ−1Aσ−1,with log-density(28)logpσ|θ=−σ−Aσ−1μσTAσTΛσ−1Aσσ−Aσ−1μσ2+c3=−σTAσTΛσ−1Aσσ2+σTAσTΛσ−1μσ+c4,where c3 and c4 are independent of σ. Then, considering the expressions in (24) and (28), one has(29)logqσ|r,ϕ,λ,θ=logqr|σ,ϕ,λ+logpσ|θ=−σTH+AσTΛσ−1Aσσ2+σTG+Hσ^+AσTΛσ−1μσ+c5,where c5 is independent of σ. Since the expression in (29) is the log-kernel of NIDμ^σ,Λ^σ in σ, where Λ^σ−1=H+AσTΛσ−1Aσ is the precision matrix or inverse of the variance and μ^σ=Λ^σG+Hσ^+AσTΛσ−1μσ is the mean vector, pσ|r,ϕ,λ,θ can be calculated by the Gaussian density(30)qσ|r,ϕ,λ,θ∼NIDμ^σ,Λ^σ.This Gaussian approximation density is then put into use as the proposal density in the acceptance-rejection Metropolis–Hastings (ARMH) step.Note that using forward- and back-substitution to simulate smoother market liquidityσ, one can calculate the (banded) Cholesky decomposition Cσ of Λ^σ such as CσCσT=Λ^σ and obtain a random draw from NIDμ^σ,Λ^σ efficiently such that σ=μ^σ+Cσu and u∼NID0,IT. A crucial result is that the Cholesky decomposition Cσ can be obtained in OT operations rather than OT3 operations for full matrices, since Λ^σ is a band diagonal matrix (see, e.g., Chan and Jeliazkov [24]).To this end, it remains to select the pointσ^ for the Taylor expansion in (21), which is chosen to be the mode of pσ|r,ϕ,λ,θ and can be quickly obtained by the numerical optimization. Using the Newton–Raphson method, one can quickly locate the maximum of qσ|r,ϕ,λ,θ, or equivalently that of logqσ|r,ϕ,λ,θ, or to select the point σ^ as the mode of pσ|r,ϕ,λ,θ or smooth σ. It follows from (29) that negative Hessian of logqσ|r,ϕ,λ,θ evaluated at σ=σ^ is Λ^σ−1 and the gradient at σ=σ^ is −Λ^σ−1σ^+Λ^σ−1μ^σ. Hence, one can carry out the Newton–Raphson method as follows: initializing with σ=σ1 for some lasting vector σ1. For i=1,2,⋯, assume σ^=σi in the evaluation of Λ^σi−1 and −Λ^σi−1σi+Λ^σi−1μ^σi and compute(31)σi+1=σi+Λ^σi−Λ^σi−1σi+Λ^σi−1μ^σi=μ^σi.Note that sinceqσ|r,ϕ,λ,θ is Gaussian density, its mode σ^ equals its mean vector μ^σ. Repeating this procedure until some convergence criterion is reached, e.g., we get the mode σ^=σi+1 or μ^σ=μ^σi when σi+1−σi<ε for some prefixed tolerance level ε>0. Also, note that the negative Hessian Λ^σ−1 is positive definite for all σ. This guarantees a fast convergence of the Newton–Raphson method.


## 3.3. Acceptance-Rejection Metropolis–Hastings (ARMH) Algorithm

As mentioned before, one problem with sampling states in a single block with the Gaussian proposal density is that the accepting rate in the MH step can be very low, and it is probably because that the Gaussian approximation density is not sufficiently precise. In this section, different from being divided into blocks, implemented in, e.g., Watanabe and Omori [35], Omori and Watanabe [36], and Abanto-Valle et al. [37], we sample the market liquidity σ in a single block and exploit an acceptance-rejection Metropolis–Hastings (ARMH) algorithm to correct the approximation error.As its name implies, the ARMH algorithm (Tierney [38]) is a MCMC sampling procedure that links classic accept-reject sampling with the Metropolis–Hastings algorithm. In the setting, the target density is the joint conditional density pσ|r,ϕ,λ,θ∝pr|σ,ϕ,λpσ|θ of σ and we have the Gaussian approximation proposal density qσ|r,ϕ,λ,θ from which we generate candidate draws σ. In the classic accepting-rejecting sampling, an important requirement is that there exists a constant c, which is selected by equation (34), such that(32)pσ|r,ϕ,λ,θ<cqσ|r,ϕ,λ,θ,for σ in support of pσ|r,ϕ,λ,θ. The ARMH releases the domination condition (32) such that when it is not satisfied for σ, one may resort to the MH algorithm. To show Algorithm 1, it is useful to define the set(33)ℂ=σ:pσ|r,ϕ,λ,θ<cqσ|r,ϕ,λ,θ,and let ℂc denote the complement code. Then, the ARMH procedure is as follows.Algorithm 1: Acceptance-rejection Metropolis–Hastings with Gaussian proposal density [8, 24, 38, 56].
(1)
Accept-reject step: generate a drawσ∗∼qσ|r,ϕ,λ,θ∼NIDμ^σ,Λ^σ and accept σ∗ with probability min1,pσ∗|r,ϕ,λ,θ/cqσ∗|r,ϕ,λ,θContinue the above process until a drawσ∗ is accepted.(2)
Metropolis–Hastings step: given the currentσ and the proposed draw σ∗(a) ifσ∈ℂ, accept σ∗ with probability 1;(b) ifσ∈ℂc and σ∗∈ℂ, accept σ∗ with probabilitymin1,cqσ|r,ϕ,λ,θ/pσ|r,ϕ,λ,θ(c) ifσ∈ℂc and σ∗∈ℂc, accept σ∗ with probabilitymin1,qσ|r,ϕ,λ,θpσ∗|r,ϕ,λ,θ/qσ∗|r,ϕ,λ,θpσ|r,ϕ,λ,θotherwise returnσ.Adjusting the original Gaussian proposal densityqσ|r,ϕ,λ,θ by the accept-reject step, a better approximation of the target density pσ|r,ϕ,λ,θ is achieved, i.e., the new proposed density is identical with the target density on the set ℂ (albeit different normalizing constants c), whereas on ℂc, the new proposed density is decreased to the original one.Chib and Jeliazkov [46] presented a practical way to select the constant c and the trade-off in such a choice, which is outlined here. Notice that if a bigger c is chosen, the set ℂ is larger and one is more likely to accept the candidate σ∗. On the other hand, selecting a smaller c will reduce the sampling efficiency, so it will generate more rejections during the sampling steps of accepting-rejecting, resulting in increased invalid and large calculations. A practical way to strike a balance between these two conflicting considerations is to set(34)c=bpσ^|r,ϕ,λ,θqσ^|r,ϕ,λ,θ,where σ^ is the mode of the conditional density qσ|r,ϕ,λ,θ and b is, say, between 1 and 5. Such kind of choice would make sure that c is sufficiently small to decrease the required number of draws from qσ|r,ϕ,λ,θ, while it is also big enough so that the set ℂ contains the mode σ^ and its neighboring points.


## 4. Bayesian Estimation for Stochastic Volatility Models

Bayesian analysis is a model of inductive inference, which has been applied in many scientific disciplines (see, e.g., Khoshravesh et al. [57]). A typical characteristic of the Bayesian approach is that it allows researchers to use sample (data) and prior (expert judgment) information in a consistent manner when making inferences. In this section, we first discuss efficient posterior simulators for the market microstructure model plus Student-t heavy tails (MMt), which are an extra block to sample the state λ and the parameter ν. Then, we outline the efficient posterior simulators for other parameters and other states, which are built on precision-based simulation smoother algorithm. In order to finish the model specification, we presume independent prior distributions γ∼NIDγ0,Vγ, α∼NIDα0,Vα, δ∼NIDδ0,Vδ, β∼NIDβ0,Vβ, η2∼IGvη2,Sη2, κ2∼IGvκ2,Sκ2, and υ∼UNIυ¯,υ¯, whereUNI⋅,⋅ denotes the uniform distribution.


### 4.1. Posterior Simulators for the Stateλ and the Parameter ν of MMt Models (19)–(22)

Sampling the latent variableλ could be easily done as λ1,λ2,…,λT, which are conditionally independently, to give the parameters and data. Each λt follows an independent inverse-gamma distribution:(35)pλt|rt,σt,ϕt∼IGυ+12,υ+d2,where d=e−σtrt−ϕteσt2.Next, sincepν|λ∝pλ|νpν and the prior for υ∼UNIυ¯,υ¯, one can first derive the log-density logpν|λ following equation (20) to sample pν|λ:(36)logpν|λ=logpλ|ν+logpν=Tν2logν/2−TlogΓν/2−ν/2+1∑i=1Tlogλt−ν/2∑i=1Tλt−1+c6,where c6 is independent of ν and Γ⋅ denotes the gamma function. It is easy to find out that the first and second derivatives of the log-density referring to ν given by(37)g=∂logpν|λ∂ν,=Tlogν2+T2−Tψν2−∑i=1Tlogλt2−∑i=1Tλt−12,h=∂2logpν|λ∂ν2=T2ν−Tψ1ν/24,where ψx=∂logΓx/∂x and ψ1x=∂2logΓx/∂x2 are the digamma function and trigamma function, respectively.Since the first and second derivatives can be calculated in a minute, one can construct a Gaussian proposal distributionqν|λ and maximize logqν|λ using numerical optimization based on the approximation of logpν|λ using second-order Taylor expansion as done in Section 3.2. One can get the mode and the negative Hessian calculated in the mode, denoted as μ^ν and Λ^ν−1, respectively. Then, the independence-chain Metropolis–Hastings step with proposal distribution is implemented as follows:(38)qν|λ∼NIDμ^ν,Λ^ν,and the proposed draw ν∗ is accepted with probability(39)min1,pν∗|λqν|λpν|λqν∗|λ.


### 4.2. Posterior Simulators for the Parametersθ=γ,δ,η2 of MMt Models (19)–(22)

The posterior simulator for the stateσ can be obtained as described in Section 3. So, one can obtain posterior simulators for the parameters θ=γ,δ,η2 following equation (19) as follows.Samplingpη2|σ,γ,δ is straightforward, if we assume the prior for η2∼IGvσ,Sσ. It yields(40)pη2|σ,γ,δ∼IGv^σ,S^σ,where v^σ=vσ+T/2, S^σ=Sσ+σ1−γ21−δ2+∑t=2Tσt−γ−δσt−1−γ2/2.Samplingpγ|σ,δ,η2 can be done easily because the conditional distribution is standard. Especially following pγ|σ,δ,η2∝pγpσ|γ,δ,η2 and the prior for γ∼NIDγ0,Vγ, one has(41)logpγ|σ,δ,η2=logpγ+∑t=1Tlogpσt|γ,δ,η2,=−γ−γ022Vγ−σ1−γ21−δ22η2,−∑t=2Tσt−δσt−1−γ1−δ22η2+c7,=−γ21/Vγ+1−δ2/η2+∑t=2T1−δ2/η22,+γγ0Vγ+σ11−δ2η2+∑t=2Tσt−δσt−11−δη2+c8,where c7 and c8 are independent of γ. As the log-density is quadratic in γ, it is Gaussian and one has(42)pγ|σ,δ,η2∼NIDμ^γ,Λ^γ,where(43)Λ^γ−1=1Vγ+1−δ2η2+∑t=2T1−δ2/η2,μ^γ=Λ^γγ0Vγ+σ11−δ2η2+∑t=2Tσt−δσt−11−δ/η2,where Λ^γ−1 is the precision matrix or inverse of variance and μ^γ is the mean.Sincepσ1|γ,δ,η2 is non-standard in δ, the conditional density pδ|σ,γ,η2 is non-standard, and one can carry out an independence-chain Metropolis–Hastings step to draw δ. Firstly, following pδ|σ,γ,η2∝pδpσ|γ,δ,η2 and the prior for δ∼NIDδ0,Vδ, one has(44)logpδ|σ,γ,η2=logpδ+∑t=1Tlogpσt|γ,δ,η2,≈logpδ+∑t=2Tlogpσt|γ,δ,η2,=−δ−δ022Vδ−∑t=2Tσt−γ−δσt−1−γ22η2+c9,=−δ21/Vδ+∑t=2Tσt−1−γ2/η22,+δδ0Vδ+∑t=2Tσt−γσt−1−γη2+c10,where c9 and c10 are independent of δ. As the approximate log-density is quadratic in δ, it is approximate Gaussian and one has(45)qδ|σ,γ,η2∼NIDμ^δ,Λ^δ,where Λ^δ−1=1/Vδ+∑t=2Tσt−1−γ2/η2 is the precision matrix or inverse of variance and μ^δ=Λ^δδ0/Vδ+∑t=2Tσt−γσt−1−γ/η2 is the mean.Then, using the approximate Gaussian densityqδ|σ,γ,η2∼NIDμ^δ,Λ^δ with δ<1 as proposal density, one can implement an independence-chain Metropolis–Hastings in which the proposal δ∗ is accepted in the feasible region with probability(46)min1,gδ∗|σ1,γ,η2gδ|σ1,γ,η2,where(47)gδ|σ1,γ,η2=−σ1−γ21−δ22η2+log1−δ22−δ−δ022Vδ+c11,where c11 is independent of δ.


### 4.3. Posterior Simulators for the Stateϕ and the Parameters θ=α,β,κ2 of MMt Models (19)–(22)

Sincepϕ|r,σ,λ,θ∝pr|σ,ϕ,λpϕ|θ, pr|σ,ϕ,λ, and pϕ|θ are all Gaussian distributions in ϕ, one immediately obtains a Gaussian distribution pϕ|r,σ,λ,θ. Now, write the return equation (17) and exceed demand equation (18) of ϕt in matrix form(48)r=Arϕ+εr,εr∼NID0,Λr,Aϕϕ=μϕ+εϕ,εϕ∼NID0,Λϕ,where Ar=diageσ1,eσ2,…,eσTT, Λr=diagλ1eσ1,λ2eσ2,…,λTeσTT, μϕ=α,1−βα,…,1−βαT, Λϕ=diagκ2/1−β2,κ2,…,κ2, and(49)Aϕ=100…0−β10…00−β1…0⋮⋱⋱⋱⋮00…−β1.Then,pϕ|θ∼NIDAϕ−1μϕ,AϕTΛϕ−1Aϕ−1 and pr|σ,ϕ,λ∼NIDArϕ,Λr with log-density(50)logpr|σ,ϕ,λ=−ϕTArTΛr−1Arϕ2+ϕTArTΛr−1r+c12,(51)logpϕ|θ=−ϕTAϕTΛϕ−1Aϕϕ2+ϕTAϕTΛϕ−1μϕ+c13,where c12 and c13 are independent of ϕ. Combining (49) and (50), one has(52)logpϕ|r,σ,λ,θ=logpr|σ,ϕ,λ+logpϕ|θ,=−ϕTArTΛr−1Ar+AϕTΛϕ−1Aϕϕ2+ϕTArTΛr−1r+AϕTΛϕ−1μϕ+c14,where c14 is independent of ϕ. As this log-density is quadratic in ϕ, it is Gaussian and one has(53)pϕ|r,σ,λ,θ∼NIDμ^ϕ,Λ^ϕ,where Λ^ϕ−1=ArTΛr−1Ar+AϕTΛϕ−1Aϕ is the precision matrix or inverse of variance and μ^ϕ=Λ^ϕArTΛr−1r+AϕTΛϕ−1μϕ is the mean.Lastly, following equations (18) and (19), α,β, and κ2 are sampled the same as those of γ,δ, and η2, respectively.Forκ2, it has inverse-gamma posterior density(54)pκ2|ϕ,α,β∼IGv^ϕ,S^ϕ,where v^ϕ=vϕ+T/2, S^ϕ=Sϕ+ϕ1−α21−β2+∑t=2Tϕt−α−βϕt−1−α2/2.Forα, one gets the Gaussian posterior density(55)pα|σ,β,κ2∼NIDμ^α,Λ^α,where(56)Λ^α−1=1Vα+1−β2κ2+∑t=2T1−β2κ2,μ^α=Λ^αα0Vα+ϕ11−β2κ2+∑t=2Tϕt−βϕt−11−βκ2,where Λ^α−1 is the precision matrix or inverse of variance and μ^α is the mean.Forβ, it has approximate Gaussian posterior density(57)qβ|σ,α,κ2∼NIDμ^β,Λ^β,where Λ^β−1=1/Vβ+∑t=2Tϕt−1−α2/κ2 is the precision matrix or inverse of variance and μ^β=Λ^ββ0/Vβ+∑t=2Tϕt−αϕt−1−α/κ2 is the mean.Then, one implements an independence-chain Metropolis–Hastings in which the proposalβ∗ is accepted in the feasible region with probability(58)min1,gβ∗|ϕ1,α,κ2gβ|ϕ1,α,κ2,where gβ|ϕ1,α,κ2=−ϕ1−α21−β2/2κ2+log1−β2/2−β−β02/2Vβ.Now an overall description of Algorithm2 is given below.Algorithm 2: The MCMC algorithm for estimation of MMt models (19)–(22) [8, 10, 24, 56].
Step 1. Initialize α,β,κ2,γ,δ,η2,v and σ,ϕ,λ;Step 2. Sample α,β,κ2,ϕ|σ,λ,r by(a)
Samplingϕ|α,β,κ2,σ,λ,r by expression (52),(b)
Samplingκ2|ϕ,α,β by expression (53),(c)
Samplingα|ϕ,β,κ2 by expression (54),(d)
Samplingβ|ϕ,α,κ2 by expression (56);Step 3. Sample γ,δ,η2,σ|ϕ,λ,r by(a)
Samplingσ|γ,δ,η2,ϕ,λ,r by expression (30),(b)
Samplingη2|σ,γ,δ by expression (39),(c)
Samplingγ|σ,δ,η2 by expression (41),(d)
Samplingδ|σ,γ,η2 by expression (44);Step 4. Sample v,λ|σ,ϕ,r by(a)
Samplingλ|v,σ,ϕ,r by expression (34),(b)
Samplingv|λ by expression (37);Step 5. Go to Step 2.


### 4.4. Bayesian Estimation for Other Stochastic Volatility Models

#### 4.4.1. For Standard Stochastic Volatility Models (1) and (2)

We firstly give an overall description of Algorithm3.Algorithm 3: The MCMC algorithm for estimation of SV models (1) and (2) [1, 24, 30, 31, 56].
Step 1. Initialize μ,γ,δ,η2 and σ;Step 2. Sample μ|σ,r;Step 3. Sample γ,δ,η2,σ|μ,r by(a)
Samplingσ|γ,δ,η2,μ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2 above, samplingpμ|r,σ can be done easily, as the conditional distribution is standard. Since the prior for μ∼NIDμ0,Vμ and pμ|r,σ∝pr|σ,μpμ, following (1), one has(59)logprt|σt,μ=−rt−μ2/2eσt+c15,and(60)logpμ|r,σ=logpμ+∑t=1Tlogprt|σt,μ,=−μ−μ022Vμ−∑t=1Trt−μ22eσt+c16,=−μ21/Vμ+∑t=1Te−σt2+μμ0Vμ+∑t=1Trte−σt+c17,where c15, c16, and c17 are independent of μ. As this log-density is quadratic in μ, it is Gaussian and one has(61)pμ|r,σ∼NIDμ^μ,Λ^μ,where Λ^μ−1=1/Vμ+∑t=1Te−σt is the precision matrix or inverse of variance and μ^μ=Λ^μμ0/Vμ+∑t=1Trte−σt is the mean.In Step 3(a), the joint conditional densitypσ|r,μ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,μ follows equation (1) instead of equation (17).(62)logprt|σt,μ=−σt2−rt−μ22eσt+c18,where c18 is independent of σt. Then, starting from equation (61) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


#### 4.4.2. For Stochastic Volatility in Mean

We firstly give an overall description of Algorithm4.Algorithm 4: The MCMC algorithm for estimation of SVM model (3) and (4) [24, 32, 56].
Step 1. Initialize ϕ,γ,δ,η2 and σ;Step 2. Sample ϕ|σ,r;Step 3. Sample γ,δ,η2,σ|ϕ,r by(a)
Samplingσ|γ,δ,η2,ϕ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2, Samplingϕ can be done easily, as the conditional distribution is standard. Since the prior for ϕ∼NIDϕ0,Vϕ and pϕ|r,σ∝pr|σ,ϕpϕ, following (3), one has(63)logprt|σt,ϕ=−rt−ϕeσt22eσt+c19,and(64)logpϕ|r,σ=logpϕ+logpr|σ,ϕ,=logpϕ+∑t=1Tlogprt|σt,ϕ,=−ϕ−ϕ022Vϕ−∑t=1Trt−ϕeσt22eσt+c20,=−ϕ21/Vϕ+∑t=1Teσt2+ϕϕ0Vϕ+∑t=1Trt+c21,where c19, c20, and c21 are independent of ϕ. As the log-density is quadratic in ϕ, it is Gaussian and one has(65)pϕ|r,σ∼NIDμ^ϕ,Λ^ϕ,where Λ^ϕ−1=1/Vϕ+∑t=1Teσt is the precision matrix or inverse of variance and μ^ϕ=Λ^ϕϕ0/Vϕ+∑t=1Trt is the mean.In Step 3(a), the joint conditional densitypσ|r,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,ϕ follows equation (6) instead of equation (17).(66)logprt|σt,ϕ=−σt2−rt−ϕeσt22eσt+c22,where c22 is independent of σt. Then, starting from equation (49) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


#### 4.4.3. For Market Microstructure Model

We firstly give an overall description of Algorithm5.Algorithm 5: The MCMC algorithm for estimation of MM models (10)–(12) [5, 10, 24, 56].
Step 1. Initialize α,β,κ2,γ,δ,η2 and σ,ϕ;Step 2. Sample α,β,κ2,ϕ|σ,r by(a)
Samplingϕ|α,β,κ2,σ,r,(b)
Samplingκ2|ϕ,α,β,(c)
Samplingα|ϕ,β,κ2,(d)
Samplingβ|ϕ,α,κ2;Step 3. Sample γ,δ,η2,σ|ϕ,r by(a)
Samplingσ|γ,δ,η2,ϕ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2(a), the joint conditional densitypϕ|r,σ,γ,δ,η2 is similar to that of the MMt model in which pr|σ,ϕ follows equation (8) instead of equation (17). One can write return equation (8) in matrix form:(67)r=Arϕ+εr,εr∼NID0,Λr,where Ar=Λr=diageσ1,eσ2,…,eσTT. Then, starting from equation (66) instead of equation (47), ϕ can be sampled as that of the MMt model.In Step 3(a), the joint conditional densitypσ|r,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,ϕt follows equation (8) instead of equation (17).(68)logprt|σt,ϕt=−σt/2−rt−ϕteσt2/2eσt+c22,where c22 is independent of σt. Then, starting from equation (67) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


#### 4.4.4. For Stochastic Volatility Model plus Student-t Heavy Tails

We firstly give an overall description of Algorithm6.Algorithm 6: The MCMC algorithm for estimation of SVt models (13)–(15) [8, 24, 43, 44, 56].
Step 1. Initialize μ,γ,δ,η2,v and σ,λ;Step 2. Sample μ|σ,λ,r;Step 3. Sample γ,δ,η2,σ|μ,λ,r by(a)
Samplingσ|γ,δ,η2,μ,λ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Sample v,λ|σ,μ,r by(a)
Samplingλ|v,σ,μ,r,(b)
Samplingv|λ;Step 5. Go to Step 2.In Step 2, samplingpμ|r,σ,λ is similar to that of the MMt model in which the log-density logprt|σt,λt,μ follows equation (11) instead of equation (1).(69)logprt|σt,λt,μ=−rt−μ2/2λteσt+c23,where c23 is independent of μ. Then, starting from equation (68) instead of equation (58), μ can be sampled as that of the SV model.In Step 3(a), the joint conditional densitypσ|r,λ,μ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,λt,μ follows equation (11) instead of equation (17).(70)logprt|σt,λt,μ=−σt2−rt−μ22λteσt+c24,where c24 is independent of σt. Then, starting from equation (69) instead of equation (21), σt can be sampled as that of the MMt model.In Step 4(a), to sampleλt is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:(71)pλt|rt,σt,μ∼IGυ+12,υ+d2,where d=e−σtrt−μ2.Lastly, other steps are the same as those of the MMt model.


#### 4.4.5. For Stochastic Volatility in Mean plus Student-t Heavy Tails

We firstly give an overall description of Algorithm7.Algorithm 7: The MCMC algorithm for estimation of SVMt models (16)–(18) [8, 24, 43, 44, 56].
Step 1. Initialize ϕ,γ,δ,η2,v and σ,λ;Step 2. Sample ϕ|σ,λ,r byStep 3. Sample γ,δ,η2,σ|ϕ,λ,r by(a)
Samplingσ|γ,δ,η2,ϕ,λ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Sample v,λ|σ,ϕ,r by(a)
Samplingλ|v,σ,ϕ,r,(b)
Samplingv|λ;Step 5. Go to Step 2.In Step 2, Samplingϕ is similar to that of the SVM model in which the log-density logprt|σt,λt,ϕ follows equation (14) instead of equation (3).(72)logprt|σt,λt,ϕ=−rt−ϕeσt22λteσt+c25,where c25 is independent of ϕ. Then, starting from equation (71) instead of equation (64), ϕ can be sampled as that of the SVM model.In Step 3(a), the joint conditional densitypσ|r,λ,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,λt,ϕ follows equation (14) instead of equation (17).(73)logprt|σt,λt,μ=−σt2−rt−μ22λteσt+c26,where c26 is independent of σt. Then, starting from equation (73) instead of equation (19), σt can be sampled as that of the MMt model.In Step 4(a), to sampleλt is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:(74)pλt|rt,σt,ϕ∼IGυ+12,υ+d2,where d=e−σtrt−ϕeσt2.Lastly, other steps are the same as those of the MMt model.In summary, this paper studied the estimation problem of six typical SV models listed in Table1, and the difference between previous research and our research in this paper is briefly explained in Table 2.Table 2
Difference between previous research and our research in this paper.
MethodsMain characteristicsKalman filter and maximum likelihood method-based approach [5, 10, 19, 29–32]Needs a large computation burden. Not good for non-Gaussian non-linear data.MCMC method based on Metropolis–Hastings sampling [56]For high-dimensional data, it has two shortcomings: one is to need quite large computation burden and another is difficulty in solving the joint distribution of variables.MCMC method based on Gibbs sampling [8, 10]Gibbs sampling is the application of the Metropolis–Hastings sampling in special cases. It does not need to judge the state reception. It has fast convergence speed and high sampling efficiency.MCMC method based on auxiliary mixture sampler [23, 24, 27]Loses precision due to transforming the return equation to a linear equation.The proposed MCMC method in this paperA fast and efficient MCMC approach based on a precision-based simulation smoother algorithm and an acceptance-rejection Metropolis–Hastings algorithm for estimation of non-linear non-Gaussian SV type models.Uses the band diagonal structure and positive definition of Hessian matrix of the logarithmic density to achieve a small computation burden and high precision.


## 4.1. Posterior Simulators for the Stateλ and the Parameter ν of MMt Models (19)–(22)

Sampling the latent variableλ could be easily done as λ1,λ2,…,λT, which are conditionally independently, to give the parameters and data. Each λt follows an independent inverse-gamma distribution:(35)pλt|rt,σt,ϕt∼IGυ+12,υ+d2,where d=e−σtrt−ϕteσt2.Next, sincepν|λ∝pλ|νpν and the prior for υ∼UNIυ¯,υ¯, one can first derive the log-density logpν|λ following equation (20) to sample pν|λ:(36)logpν|λ=logpλ|ν+logpν=Tν2logν/2−TlogΓν/2−ν/2+1∑i=1Tlogλt−ν/2∑i=1Tλt−1+c6,where c6 is independent of ν and Γ⋅ denotes the gamma function. It is easy to find out that the first and second derivatives of the log-density referring to ν given by(37)g=∂logpν|λ∂ν,=Tlogν2+T2−Tψν2−∑i=1Tlogλt2−∑i=1Tλt−12,h=∂2logpν|λ∂ν2=T2ν−Tψ1ν/24,where ψx=∂logΓx/∂x and ψ1x=∂2logΓx/∂x2 are the digamma function and trigamma function, respectively.Since the first and second derivatives can be calculated in a minute, one can construct a Gaussian proposal distributionqν|λ and maximize logqν|λ using numerical optimization based on the approximation of logpν|λ using second-order Taylor expansion as done in Section 3.2. One can get the mode and the negative Hessian calculated in the mode, denoted as μ^ν and Λ^ν−1, respectively. Then, the independence-chain Metropolis–Hastings step with proposal distribution is implemented as follows:(38)qν|λ∼NIDμ^ν,Λ^ν,and the proposed draw ν∗ is accepted with probability(39)min1,pν∗|λqν|λpν|λqν∗|λ.


## 4.2. Posterior Simulators for the Parametersθ=γ,δ,η2 of MMt Models (19)–(22)

The posterior simulator for the stateσ can be obtained as described in Section 3. So, one can obtain posterior simulators for the parameters θ=γ,δ,η2 following equation (19) as follows.Samplingpη2|σ,γ,δ is straightforward, if we assume the prior for η2∼IGvσ,Sσ. It yields(40)pη2|σ,γ,δ∼IGv^σ,S^σ,where v^σ=vσ+T/2, S^σ=Sσ+σ1−γ21−δ2+∑t=2Tσt−γ−δσt−1−γ2/2.Samplingpγ|σ,δ,η2 can be done easily because the conditional distribution is standard. Especially following pγ|σ,δ,η2∝pγpσ|γ,δ,η2 and the prior for γ∼NIDγ0,Vγ, one has(41)logpγ|σ,δ,η2=logpγ+∑t=1Tlogpσt|γ,δ,η2,=−γ−γ022Vγ−σ1−γ21−δ22η2,−∑t=2Tσt−δσt−1−γ1−δ22η2+c7,=−γ21/Vγ+1−δ2/η2+∑t=2T1−δ2/η22,+γγ0Vγ+σ11−δ2η2+∑t=2Tσt−δσt−11−δη2+c8,where c7 and c8 are independent of γ. As the log-density is quadratic in γ, it is Gaussian and one has(42)pγ|σ,δ,η2∼NIDμ^γ,Λ^γ,where(43)Λ^γ−1=1Vγ+1−δ2η2+∑t=2T1−δ2/η2,μ^γ=Λ^γγ0Vγ+σ11−δ2η2+∑t=2Tσt−δσt−11−δ/η2,where Λ^γ−1 is the precision matrix or inverse of variance and μ^γ is the mean.Sincepσ1|γ,δ,η2 is non-standard in δ, the conditional density pδ|σ,γ,η2 is non-standard, and one can carry out an independence-chain Metropolis–Hastings step to draw δ. Firstly, following pδ|σ,γ,η2∝pδpσ|γ,δ,η2 and the prior for δ∼NIDδ0,Vδ, one has(44)logpδ|σ,γ,η2=logpδ+∑t=1Tlogpσt|γ,δ,η2,≈logpδ+∑t=2Tlogpσt|γ,δ,η2,=−δ−δ022Vδ−∑t=2Tσt−γ−δσt−1−γ22η2+c9,=−δ21/Vδ+∑t=2Tσt−1−γ2/η22,+δδ0Vδ+∑t=2Tσt−γσt−1−γη2+c10,where c9 and c10 are independent of δ. As the approximate log-density is quadratic in δ, it is approximate Gaussian and one has(45)qδ|σ,γ,η2∼NIDμ^δ,Λ^δ,where Λ^δ−1=1/Vδ+∑t=2Tσt−1−γ2/η2 is the precision matrix or inverse of variance and μ^δ=Λ^δδ0/Vδ+∑t=2Tσt−γσt−1−γ/η2 is the mean.Then, using the approximate Gaussian densityqδ|σ,γ,η2∼NIDμ^δ,Λ^δ with δ<1 as proposal density, one can implement an independence-chain Metropolis–Hastings in which the proposal δ∗ is accepted in the feasible region with probability(46)min1,gδ∗|σ1,γ,η2gδ|σ1,γ,η2,where(47)gδ|σ1,γ,η2=−σ1−γ21−δ22η2+log1−δ22−δ−δ022Vδ+c11,where c11 is independent of δ.


## 4.3. Posterior Simulators for the Stateϕ and the Parameters θ=α,β,κ2 of MMt Models (19)–(22)

Sincepϕ|r,σ,λ,θ∝pr|σ,ϕ,λpϕ|θ, pr|σ,ϕ,λ, and pϕ|θ are all Gaussian distributions in ϕ, one immediately obtains a Gaussian distribution pϕ|r,σ,λ,θ. Now, write the return equation (17) and exceed demand equation (18) of ϕt in matrix form(48)r=Arϕ+εr,εr∼NID0,Λr,Aϕϕ=μϕ+εϕ,εϕ∼NID0,Λϕ,where Ar=diageσ1,eσ2,…,eσTT, Λr=diagλ1eσ1,λ2eσ2,…,λTeσTT, μϕ=α,1−βα,…,1−βαT, Λϕ=diagκ2/1−β2,κ2,…,κ2, and(49)Aϕ=100…0−β10…00−β1…0⋮⋱⋱⋱⋮00…−β1.Then,pϕ|θ∼NIDAϕ−1μϕ,AϕTΛϕ−1Aϕ−1 and pr|σ,ϕ,λ∼NIDArϕ,Λr with log-density(50)logpr|σ,ϕ,λ=−ϕTArTΛr−1Arϕ2+ϕTArTΛr−1r+c12,(51)logpϕ|θ=−ϕTAϕTΛϕ−1Aϕϕ2+ϕTAϕTΛϕ−1μϕ+c13,where c12 and c13 are independent of ϕ. Combining (49) and (50), one has(52)logpϕ|r,σ,λ,θ=logpr|σ,ϕ,λ+logpϕ|θ,=−ϕTArTΛr−1Ar+AϕTΛϕ−1Aϕϕ2+ϕTArTΛr−1r+AϕTΛϕ−1μϕ+c14,where c14 is independent of ϕ. As this log-density is quadratic in ϕ, it is Gaussian and one has(53)pϕ|r,σ,λ,θ∼NIDμ^ϕ,Λ^ϕ,where Λ^ϕ−1=ArTΛr−1Ar+AϕTΛϕ−1Aϕ is the precision matrix or inverse of variance and μ^ϕ=Λ^ϕArTΛr−1r+AϕTΛϕ−1μϕ is the mean.Lastly, following equations (18) and (19), α,β, and κ2 are sampled the same as those of γ,δ, and η2, respectively.Forκ2, it has inverse-gamma posterior density(54)pκ2|ϕ,α,β∼IGv^ϕ,S^ϕ,where v^ϕ=vϕ+T/2, S^ϕ=Sϕ+ϕ1−α21−β2+∑t=2Tϕt−α−βϕt−1−α2/2.Forα, one gets the Gaussian posterior density(55)pα|σ,β,κ2∼NIDμ^α,Λ^α,where(56)Λ^α−1=1Vα+1−β2κ2+∑t=2T1−β2κ2,μ^α=Λ^αα0Vα+ϕ11−β2κ2+∑t=2Tϕt−βϕt−11−βκ2,where Λ^α−1 is the precision matrix or inverse of variance and μ^α is the mean.Forβ, it has approximate Gaussian posterior density(57)qβ|σ,α,κ2∼NIDμ^β,Λ^β,where Λ^β−1=1/Vβ+∑t=2Tϕt−1−α2/κ2 is the precision matrix or inverse of variance and μ^β=Λ^ββ0/Vβ+∑t=2Tϕt−αϕt−1−α/κ2 is the mean.Then, one implements an independence-chain Metropolis–Hastings in which the proposalβ∗ is accepted in the feasible region with probability(58)min1,gβ∗|ϕ1,α,κ2gβ|ϕ1,α,κ2,where gβ|ϕ1,α,κ2=−ϕ1−α21−β2/2κ2+log1−β2/2−β−β02/2Vβ.Now an overall description of Algorithm2 is given below.Algorithm 2: The MCMC algorithm for estimation of MMt models (19)–(22) [8, 10, 24, 56].
Step 1. Initialize α,β,κ2,γ,δ,η2,v and σ,ϕ,λ;Step 2. Sample α,β,κ2,ϕ|σ,λ,r by(a)
Samplingϕ|α,β,κ2,σ,λ,r by expression (52),(b)
Samplingκ2|ϕ,α,β by expression (53),(c)
Samplingα|ϕ,β,κ2 by expression (54),(d)
Samplingβ|ϕ,α,κ2 by expression (56);Step 3. Sample γ,δ,η2,σ|ϕ,λ,r by(a)
Samplingσ|γ,δ,η2,ϕ,λ,r by expression (30),(b)
Samplingη2|σ,γ,δ by expression (39),(c)
Samplingγ|σ,δ,η2 by expression (41),(d)
Samplingδ|σ,γ,η2 by expression (44);Step 4. Sample v,λ|σ,ϕ,r by(a)
Samplingλ|v,σ,ϕ,r by expression (34),(b)
Samplingv|λ by expression (37);Step 5. Go to Step 2.


## 4.4. Bayesian Estimation for Other Stochastic Volatility Models

### 4.4.1. For Standard Stochastic Volatility Models (1) and (2)

We firstly give an overall description of Algorithm3.Algorithm 3: The MCMC algorithm for estimation of SV models (1) and (2) [1, 24, 30, 31, 56].
Step 1. Initialize μ,γ,δ,η2 and σ;Step 2. Sample μ|σ,r;Step 3. Sample γ,δ,η2,σ|μ,r by(a)
Samplingσ|γ,δ,η2,μ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2 above, samplingpμ|r,σ can be done easily, as the conditional distribution is standard. Since the prior for μ∼NIDμ0,Vμ and pμ|r,σ∝pr|σ,μpμ, following (1), one has(59)logprt|σt,μ=−rt−μ2/2eσt+c15,and(60)logpμ|r,σ=logpμ+∑t=1Tlogprt|σt,μ,=−μ−μ022Vμ−∑t=1Trt−μ22eσt+c16,=−μ21/Vμ+∑t=1Te−σt2+μμ0Vμ+∑t=1Trte−σt+c17,where c15, c16, and c17 are independent of μ. As this log-density is quadratic in μ, it is Gaussian and one has(61)pμ|r,σ∼NIDμ^μ,Λ^μ,where Λ^μ−1=1/Vμ+∑t=1Te−σt is the precision matrix or inverse of variance and μ^μ=Λ^μμ0/Vμ+∑t=1Trte−σt is the mean.In Step 3(a), the joint conditional densitypσ|r,μ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,μ follows equation (1) instead of equation (17).(62)logprt|σt,μ=−σt2−rt−μ22eσt+c18,where c18 is independent of σt. Then, starting from equation (61) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


### 4.4.2. For Stochastic Volatility in Mean

We firstly give an overall description of Algorithm4.Algorithm 4: The MCMC algorithm for estimation of SVM model (3) and (4) [24, 32, 56].
Step 1. Initialize ϕ,γ,δ,η2 and σ;Step 2. Sample ϕ|σ,r;Step 3. Sample γ,δ,η2,σ|ϕ,r by(a)
Samplingσ|γ,δ,η2,ϕ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2, Samplingϕ can be done easily, as the conditional distribution is standard. Since the prior for ϕ∼NIDϕ0,Vϕ and pϕ|r,σ∝pr|σ,ϕpϕ, following (3), one has(63)logprt|σt,ϕ=−rt−ϕeσt22eσt+c19,and(64)logpϕ|r,σ=logpϕ+logpr|σ,ϕ,=logpϕ+∑t=1Tlogprt|σt,ϕ,=−ϕ−ϕ022Vϕ−∑t=1Trt−ϕeσt22eσt+c20,=−ϕ21/Vϕ+∑t=1Teσt2+ϕϕ0Vϕ+∑t=1Trt+c21,where c19, c20, and c21 are independent of ϕ. As the log-density is quadratic in ϕ, it is Gaussian and one has(65)pϕ|r,σ∼NIDμ^ϕ,Λ^ϕ,where Λ^ϕ−1=1/Vϕ+∑t=1Teσt is the precision matrix or inverse of variance and μ^ϕ=Λ^ϕϕ0/Vϕ+∑t=1Trt is the mean.In Step 3(a), the joint conditional densitypσ|r,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,ϕ follows equation (6) instead of equation (17).(66)logprt|σt,ϕ=−σt2−rt−ϕeσt22eσt+c22,where c22 is independent of σt. Then, starting from equation (49) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


### 4.4.3. For Market Microstructure Model

We firstly give an overall description of Algorithm5.Algorithm 5: The MCMC algorithm for estimation of MM models (10)–(12) [5, 10, 24, 56].
Step 1. Initialize α,β,κ2,γ,δ,η2 and σ,ϕ;Step 2. Sample α,β,κ2,ϕ|σ,r by(a)
Samplingϕ|α,β,κ2,σ,r,(b)
Samplingκ2|ϕ,α,β,(c)
Samplingα|ϕ,β,κ2,(d)
Samplingβ|ϕ,α,κ2;Step 3. Sample γ,δ,η2,σ|ϕ,r by(a)
Samplingσ|γ,δ,η2,ϕ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2(a), the joint conditional densitypϕ|r,σ,γ,δ,η2 is similar to that of the MMt model in which pr|σ,ϕ follows equation (8) instead of equation (17). One can write return equation (8) in matrix form:(67)r=Arϕ+εr,εr∼NID0,Λr,where Ar=Λr=diageσ1,eσ2,…,eσTT. Then, starting from equation (66) instead of equation (47), ϕ can be sampled as that of the MMt model.In Step 3(a), the joint conditional densitypσ|r,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,ϕt follows equation (8) instead of equation (17).(68)logprt|σt,ϕt=−σt/2−rt−ϕteσt2/2eσt+c22,where c22 is independent of σt. Then, starting from equation (67) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


### 4.4.4. For Stochastic Volatility Model plus Student-t Heavy Tails

We firstly give an overall description of Algorithm6.Algorithm 6: The MCMC algorithm for estimation of SVt models (13)–(15) [8, 24, 43, 44, 56].
Step 1. Initialize μ,γ,δ,η2,v and σ,λ;Step 2. Sample μ|σ,λ,r;Step 3. Sample γ,δ,η2,σ|μ,λ,r by(a)
Samplingσ|γ,δ,η2,μ,λ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Sample v,λ|σ,μ,r by(a)
Samplingλ|v,σ,μ,r,(b)
Samplingv|λ;Step 5. Go to Step 2.In Step 2, samplingpμ|r,σ,λ is similar to that of the MMt model in which the log-density logprt|σt,λt,μ follows equation (11) instead of equation (1).(69)logprt|σt,λt,μ=−rt−μ2/2λteσt+c23,where c23 is independent of μ. Then, starting from equation (68) instead of equation (58), μ can be sampled as that of the SV model.In Step 3(a), the joint conditional densitypσ|r,λ,μ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,λt,μ follows equation (11) instead of equation (17).(70)logprt|σt,λt,μ=−σt2−rt−μ22λteσt+c24,where c24 is independent of σt. Then, starting from equation (69) instead of equation (21), σt can be sampled as that of the MMt model.In Step 4(a), to sampleλt is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:(71)pλt|rt,σt,μ∼IGυ+12,υ+d2,where d=e−σtrt−μ2.Lastly, other steps are the same as those of the MMt model.


### 4.4.5. For Stochastic Volatility in Mean plus Student-t Heavy Tails

We firstly give an overall description of Algorithm7.Algorithm 7: The MCMC algorithm for estimation of SVMt models (16)–(18) [8, 24, 43, 44, 56].
Step 1. Initialize ϕ,γ,δ,η2,v and σ,λ;Step 2. Sample ϕ|σ,λ,r byStep 3. Sample γ,δ,η2,σ|ϕ,λ,r by(a)
Samplingσ|γ,δ,η2,ϕ,λ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Sample v,λ|σ,ϕ,r by(a)
Samplingλ|v,σ,ϕ,r,(b)
Samplingv|λ;Step 5. Go to Step 2.In Step 2, Samplingϕ is similar to that of the SVM model in which the log-density logprt|σt,λt,ϕ follows equation (14) instead of equation (3).(72)logprt|σt,λt,ϕ=−rt−ϕeσt22λteσt+c25,where c25 is independent of ϕ. Then, starting from equation (71) instead of equation (64), ϕ can be sampled as that of the SVM model.In Step 3(a), the joint conditional densitypσ|r,λ,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,λt,ϕ follows equation (14) instead of equation (17).(73)logprt|σt,λt,μ=−σt2−rt−μ22λteσt+c26,where c26 is independent of σt. Then, starting from equation (73) instead of equation (19), σt can be sampled as that of the MMt model.In Step 4(a), to sampleλt is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:(74)pλt|rt,σt,ϕ∼IGυ+12,υ+d2,where d=e−σtrt−ϕeσt2.Lastly, other steps are the same as those of the MMt model.In summary, this paper studied the estimation problem of six typical SV models listed in Table1, and the difference between previous research and our research in this paper is briefly explained in Table 2.Table 2
Difference between previous research and our research in this paper.
MethodsMain characteristicsKalman filter and maximum likelihood method-based approach [5, 10, 19, 29–32]Needs a large computation burden. Not good for non-Gaussian non-linear data.MCMC method based on Metropolis–Hastings sampling [56]For high-dimensional data, it has two shortcomings: one is to need quite large computation burden and another is difficulty in solving the joint distribution of variables.MCMC method based on Gibbs sampling [8, 10]Gibbs sampling is the application of the Metropolis–Hastings sampling in special cases. It does not need to judge the state reception. It has fast convergence speed and high sampling efficiency.MCMC method based on auxiliary mixture sampler [23, 24, 27]Loses precision due to transforming the return equation to a linear equation.The proposed MCMC method in this paperA fast and efficient MCMC approach based on a precision-based simulation smoother algorithm and an acceptance-rejection Metropolis–Hastings algorithm for estimation of non-linear non-Gaussian SV type models.Uses the band diagonal structure and positive definition of Hessian matrix of the logarithmic density to achieve a small computation burden and high precision.


## 4.4.1. For Standard Stochastic Volatility Models (1) and (2)

We firstly give an overall description of Algorithm3.Algorithm 3: The MCMC algorithm for estimation of SV models (1) and (2) [1, 24, 30, 31, 56].
Step 1. Initialize μ,γ,δ,η2 and σ;Step 2. Sample μ|σ,r;Step 3. Sample γ,δ,η2,σ|μ,r by(a)
Samplingσ|γ,δ,η2,μ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2 above, samplingpμ|r,σ can be done easily, as the conditional distribution is standard. Since the prior for μ∼NIDμ0,Vμ and pμ|r,σ∝pr|σ,μpμ, following (1), one has(59)logprt|σt,μ=−rt−μ2/2eσt+c15,and(60)logpμ|r,σ=logpμ+∑t=1Tlogprt|σt,μ,=−μ−μ022Vμ−∑t=1Trt−μ22eσt+c16,=−μ21/Vμ+∑t=1Te−σt2+μμ0Vμ+∑t=1Trte−σt+c17,where c15, c16, and c17 are independent of μ. As this log-density is quadratic in μ, it is Gaussian and one has(61)pμ|r,σ∼NIDμ^μ,Λ^μ,where Λ^μ−1=1/Vμ+∑t=1Te−σt is the precision matrix or inverse of variance and μ^μ=Λ^μμ0/Vμ+∑t=1Trte−σt is the mean.In Step 3(a), the joint conditional densitypσ|r,μ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,μ follows equation (1) instead of equation (17).(62)logprt|σt,μ=−σt2−rt−μ22eσt+c18,where c18 is independent of σt. Then, starting from equation (61) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


## 4.4.2. For Stochastic Volatility in Mean

We firstly give an overall description of Algorithm4.Algorithm 4: The MCMC algorithm for estimation of SVM model (3) and (4) [24, 32, 56].
Step 1. Initialize ϕ,γ,δ,η2 and σ;Step 2. Sample ϕ|σ,r;Step 3. Sample γ,δ,η2,σ|ϕ,r by(a)
Samplingσ|γ,δ,η2,ϕ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2, Samplingϕ can be done easily, as the conditional distribution is standard. Since the prior for ϕ∼NIDϕ0,Vϕ and pϕ|r,σ∝pr|σ,ϕpϕ, following (3), one has(63)logprt|σt,ϕ=−rt−ϕeσt22eσt+c19,and(64)logpϕ|r,σ=logpϕ+logpr|σ,ϕ,=logpϕ+∑t=1Tlogprt|σt,ϕ,=−ϕ−ϕ022Vϕ−∑t=1Trt−ϕeσt22eσt+c20,=−ϕ21/Vϕ+∑t=1Teσt2+ϕϕ0Vϕ+∑t=1Trt+c21,where c19, c20, and c21 are independent of ϕ. As the log-density is quadratic in ϕ, it is Gaussian and one has(65)pϕ|r,σ∼NIDμ^ϕ,Λ^ϕ,where Λ^ϕ−1=1/Vϕ+∑t=1Teσt is the precision matrix or inverse of variance and μ^ϕ=Λ^ϕϕ0/Vϕ+∑t=1Trt is the mean.In Step 3(a), the joint conditional densitypσ|r,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,ϕ follows equation (6) instead of equation (17).(66)logprt|σt,ϕ=−σt2−rt−ϕeσt22eσt+c22,where c22 is independent of σt. Then, starting from equation (49) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


## 4.4.3. For Market Microstructure Model

We firstly give an overall description of Algorithm5.Algorithm 5: The MCMC algorithm for estimation of MM models (10)–(12) [5, 10, 24, 56].
Step 1. Initialize α,β,κ2,γ,δ,η2 and σ,ϕ;Step 2. Sample α,β,κ2,ϕ|σ,r by(a)
Samplingϕ|α,β,κ2,σ,r,(b)
Samplingκ2|ϕ,α,β,(c)
Samplingα|ϕ,β,κ2,(d)
Samplingβ|ϕ,α,κ2;Step 3. Sample γ,δ,η2,σ|ϕ,r by(a)
Samplingσ|γ,δ,η2,ϕ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Go to Step 2.In Step 2(a), the joint conditional densitypϕ|r,σ,γ,δ,η2 is similar to that of the MMt model in which pr|σ,ϕ follows equation (8) instead of equation (17). One can write return equation (8) in matrix form:(67)r=Arϕ+εr,εr∼NID0,Λr,where Ar=Λr=diageσ1,eσ2,…,eσTT. Then, starting from equation (66) instead of equation (47), ϕ can be sampled as that of the MMt model.In Step 3(a), the joint conditional densitypσ|r,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,ϕt follows equation (8) instead of equation (17).(68)logprt|σt,ϕt=−σt/2−rt−ϕteσt2/2eσt+c22,where c22 is independent of σt. Then, starting from equation (67) instead of equation (21), σt can be sampled as that of the MMt model.Lastly, other steps are the same as those of the MMt model.


## 4.4.4. For Stochastic Volatility Model plus Student-t Heavy Tails

We firstly give an overall description of Algorithm6.Algorithm 6: The MCMC algorithm for estimation of SVt models (13)–(15) [8, 24, 43, 44, 56].
Step 1. Initialize μ,γ,δ,η2,v and σ,λ;Step 2. Sample μ|σ,λ,r;Step 3. Sample γ,δ,η2,σ|μ,λ,r by(a)
Samplingσ|γ,δ,η2,μ,λ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Sample v,λ|σ,μ,r by(a)
Samplingλ|v,σ,μ,r,(b)
Samplingv|λ;Step 5. Go to Step 2.In Step 2, samplingpμ|r,σ,λ is similar to that of the MMt model in which the log-density logprt|σt,λt,μ follows equation (11) instead of equation (1).(69)logprt|σt,λt,μ=−rt−μ2/2λteσt+c23,where c23 is independent of μ. Then, starting from equation (68) instead of equation (58), μ can be sampled as that of the SV model.In Step 3(a), the joint conditional densitypσ|r,λ,μ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,λt,μ follows equation (11) instead of equation (17).(70)logprt|σt,λt,μ=−σt2−rt−μ22λteσt+c24,where c24 is independent of σt. Then, starting from equation (69) instead of equation (21), σt can be sampled as that of the MMt model.In Step 4(a), to sampleλt is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:(71)pλt|rt,σt,μ∼IGυ+12,υ+d2,where d=e−σtrt−μ2.Lastly, other steps are the same as those of the MMt model.


## 4.4.5. For Stochastic Volatility in Mean plus Student-t Heavy Tails

We firstly give an overall description of Algorithm7.Algorithm 7: The MCMC algorithm for estimation of SVMt models (16)–(18) [8, 24, 43, 44, 56].
Step 1. Initialize ϕ,γ,δ,η2,v and σ,λ;Step 2. Sample ϕ|σ,λ,r byStep 3. Sample γ,δ,η2,σ|ϕ,λ,r by(a)
Samplingσ|γ,δ,η2,ϕ,λ,r,(b)
Samplingη2|σ,γ,δ,(c)
Samplingγ|σ,δ,η2,(d)
Samplingδ|σ,γ,η2;Step 4. Sample v,λ|σ,ϕ,r by(a)
Samplingλ|v,σ,ϕ,r,(b)
Samplingv|λ;Step 5. Go to Step 2.In Step 2, Samplingϕ is similar to that of the SVM model in which the log-density logprt|σt,λt,ϕ follows equation (14) instead of equation (3).(72)logprt|σt,λt,ϕ=−rt−ϕeσt22λteσt+c25,where c25 is independent of ϕ. Then, starting from equation (71) instead of equation (64), ϕ can be sampled as that of the SVM model.In Step 3(a), the joint conditional densitypσ|r,λ,ϕ,γ,δ,η2 is similar to that of the MMt model in which the log-density logprt|σt,λt,ϕ follows equation (14) instead of equation (17).(73)logprt|σt,λt,μ=−σt2−rt−μ22λteσt+c26,where c26 is independent of σt. Then, starting from equation (73) instead of equation (19), σt can be sampled as that of the MMt model.In Step 4(a), to sampleλt is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:(74)pλt|rt,σt,ϕ∼IGυ+12,υ+d2,where d=e−σtrt−ϕeσt2.Lastly, other steps are the same as those of the MMt model.In summary, this paper studied the estimation problem of six typical SV models listed in Table1, and the difference between previous research and our research in this paper is briefly explained in Table 2.Table 2
Difference between previous research and our research in this paper.
MethodsMain characteristicsKalman filter and maximum likelihood method-based approach [5, 10, 19, 29–32]Needs a large computation burden. Not good for non-Gaussian non-linear data.MCMC method based on Metropolis–Hastings sampling [56]For high-dimensional data, it has two shortcomings: one is to need quite large computation burden and another is difficulty in solving the joint distribution of variables.MCMC method based on Gibbs sampling [8, 10]Gibbs sampling is the application of the Metropolis–Hastings sampling in special cases. It does not need to judge the state reception. It has fast convergence speed and high sampling efficiency.MCMC method based on auxiliary mixture sampler [23, 24, 27]Loses precision due to transforming the return equation to a linear equation.The proposed MCMC method in this paperA fast and efficient MCMC approach based on a precision-based simulation smoother algorithm and an acceptance-rejection Metropolis–Hastings algorithm for estimation of non-linear non-Gaussian SV type models.Uses the band diagonal structure and positive definition of Hessian matrix of the logarithmic density to achieve a small computation burden and high precision.


## 5. Simulation Research

With the purpose of evaluating the performance of the proposed precision-based simulation smoother algorithm and Markov Chain Monte Carlo method, we first utilize the simulated datasets to check the method. In this section, we only consider the SV model and the MM model; since the performance of the SV model and the MM model could compare to a voluminous literature, the comparison results can show the effectiveness of our method. In the simulation design, hundred datasets each composed ofT=1500 observations are generated from the two models. For each dataset, we estimate the two models by constructing Markov chain of length 10000 after a burn-in period of 10000.The parameter values are chosen to be comparable to those estimated from financial daily returns data (measured in percentage). They are also similar to those used in other simulation studies, like those in Kim et al. [19], Jacquier et al. [43], Omori et al. [20], and Pitt et al. [58]. In particular, we set μ=0 for the SV model and α=0 for the MM model. Parameters for the log-volatility or market liquidity transition are made to be γ=−0.25, δ=0.97, and η2=0.06 for both the SV model and the MM model. Parameters for the surplus demand transition are made to be β=0.6 and κ2=0.005 for the MM model. These parameters above are chosen by referring to the results given in Tables 3 and 4 in Section 6 for guaranteeing the reasonability of the chosen parameters. These parameters are used to build a SV model and MM model for generating observation data for simulation study so as to evaluate the performance of the proposed precision-based simulation smoother algorithm and MCMC method. In this section, we utilize the simulated datasets generated by the SV model and MM model with the chosen parameters to check the proposed method.Table 3
Parameter posterior means, standard deviations (in parentheses), inefficiency factors (in braces), and 95% intervals (in square brackets) for the SV and SVM models using Kalman filter-based simulation smoother algorithm.
SVSVMΜ0.1051 (0.01924) {2.149} [0.06650, 0.1431]ϕ0.041144 (0.022499) {1.172} [−0.002488, 0.084986]γ−0.2456 (0.2406) {1.383} [−0.7101, 0.2373]−0.2725 (0.2553) {1.542} [−0.7593, 0.2532]δ0.9691 (0.008464) {22.89} [0.9515, 0.9847]0.9742 (0.007807) {20.99} [0.9578, 0.9883]η20.06961 (0.01349) {59.73} [0.04616, 0.09785]0.05223 (0.01067) {87.73} [0.03515, 0.07667]ℚLB16.72 (2.212)12.74 (2.040)ℚML18.70 (4.926)18.17 (5.047)Note. ℚLB and ℚML are the Ljung–Box and McLeod–Li statistics of order 20 computed on the standard residuals and squared standard residuals, and the numerical standard mistakes are in parentheses. The 0.05 and 0.01 important values are 31.41 and 37.57, respectively.Table 4
Parameter posterior means, standard deviations (in parentheses), inefficiency factors (in braces), and 95% intervals (in square brackets) for the stochastic volatility models by using the precision-based simulation smoother method.
SVSVMMMΜ0.1036 (0.01878) {2.411} [0.06645, 0.1409]ϕ0.04071 (0.02261) {1.177} [−0.003155, 0.08544]γ−0.2772 (0.2450) {1.998} [−0.7547, 0.2248]−0.2358 (0.2660) {1.776} [−0.7216, 0.2865]−0.2418 (0.2663) {1.481} [−0.7388, 0.2973]δ0.97079 (0.008393) {26.80} [0.9531, 0.9861]0.9748 (0.007572) {19.72} [0.9588, 0.9885]0.9749 (0.007630) {23.73} [0.95901, 0.9885]η20.06342 (0.01289) {65.44} [0.04163, 0.09177]0.04991 (0.009473){53.76} [0.03414, 0.07103]0.04931 (0.01004) {75.65} [0.03220, 0.07193]α0.04437 (0.02477) {50.46} [−0.004412, 0.09382]β0.6375 (0.06925) {27.89} [0.4976, 0.7694]κ20.004579 (0.001315) {99.00} [0.002510, 0.007710]ℚLB16.78 (2.167)12.77 (1.966)13.69 (2.299)ℚML18.39 (4.904)17.84 (4.741)17.43 (4.809)Note. ℚLB and ℚML are, respectively, the Ljung–Box and McLeod–Li statistics of order 20 computed on the standard residuals and squared standard residuals, and the numerical standard mistakes are in parentheses. The 0.05 and 0.01 critical values are 31.41 and 37.57, respectively.The parameter prior distributions studied in Section4 are considered. We select the same hyperparameters for the parameters that are common across two models. Moreover, the hyperparameters are selected so that the implied prior means are similar to the estimates of typical financial daily return data. Following γ∼NIDγ0,Vγ, δ∼NIDδ0,Vδ, β∼NIDβ0,Vβ, η2∼IGvη2,Sη2=IGvσ,Sσ, and κ2∼IGvκ2,Sκ2=IGvϕ,Sϕ, we set γ0=0, Vγ=10, δ0=0.97, Vδ=0.01, vσ=11, Sσ=0.6, β0=0.7, Vβ=0.01, vϕ=11, and Sϕ=0.05. These parameters are the hyperparameters of the conjugate prior distributions for the model parameters, which describe the prior information of the model parameters. As a MCMC Bayesian method, our method uses the observation data and prior information of model parameters to infer the posterior information of model parameters. The hyperparameters of the prior distributions may be given intuitively. If you set enough long MCMC and long dataset, it should result in the same hyperparameters for the posterior distribution which is only contingent on the observation data and is independent of the hyperparameters for prior distributions.Empirically, in order to accelerate MCMC convergence, in this simulation study, as the regression coefficientδ of the built SV models (1) and (2) and MM models (10)–(12) is set as δ=0.97, for the hyperparameters δ0,Vδ, δ∼NIDδ0,Vδ, we set the mean δ0=0.97 and a small variance Vδ=0.01 for making the MCMC convergence fast.Because the regression coefficientβ of the built MM models (10)–(12) is set as β=0.6 in this simulation study, for the hyperparameters β0,Vβ, β∼NIDβ0,Vβ, we set the mean β0=0.7that is close to the regression coefficientβ and a small variance Vβ=0.01 for making the MCMC convergence fast.Because the regression coefficientγ of the built SV models (1) and (2) and MM models (10)–(12) is set as γ=−0.25 in this simulation study, for the hyperparameters γ0,Vγ, γ∼NIDγ0,Vγ, we set the mean γ0=0that is close to the regression coefficientγ and a simple large variance Vγ=10 since the mean of the volatility process, γ∼NIDγ0,Vγ, has large disparity for different financial data; usually, it is small for exchange rate and is large for stock.For the volatility of the state,η2∼IGvσ,Sσ=IG11,0.6 and κ2∼IGvϕ,Sϕ=IG11,0.05, in financial economic viewpoint for investment arbitrage, the smaller the volatility of the state, the better, so we select sufficient simple small mean S/v−1 and variance S2/v−12v−2 for the η2andκ2, and thus we set η2=Sσ/vσ−1=0.6/11−1=0.06 and κ2=Sϕ/vϕ−1=0.05/11−1=0.005 for the built SV models (1) and (2) and MM models (10)–(12) that are used to generate the observation data for simulation study.Figure1 depicts the posterior densities of the parameters for the SV model and the MM model, respectively, and it shows that the most posterior samples are close to the true values, and then it is clear that all parameters are estimated very precisely with small variations across simulated samples (see e.g., Kim et al. [19]; Chib et al. [56]; Omori et al. [20]; and Nakajima and Omori [21]). The accepting rate of the acceptance-rejection Metropolis–Hastings step to sample σ is 89.4% and 89.3% for the SV model and the MM model, respectively. This indicates that the Gaussian proposal density is quite similar to the conditional posterior density of σ, and the precision-based simulation smoothing algorithm is effective. The accepting rates of the Metropolis–Hastings step to sample δ and β are all larger than 90%, and it indicates the truncate Gaussian proposal density will be close to the conditional posterior density of δ and β.Figure 1
The posterior densities of parameters for SV and MM models.Tables5 and 6 show the estimated posterior means, the standard deviations, the 95% credible intervals, and inefficiency factors calculated in the simulated datasets. The posterior means are very close to the true values, and all true values are contained in the 95% credible intervals. The posterior means are computed based on 100 independent datasets, each of which is of length 10000 after a burn-in period of 10000. Using a total of one million draws for each model, the estimated parameters are receivable and accurate.Table 5
The estimated result for SV model using simulated datasets.
TruePriorPosteriorCredible intervalsInefficiencyMeanStdMeanStdΓ−0.250.003.162−0.29300.3290(−0.9744, 0.3179)1.305Δ0.970.970.1000.96900.01191(0.9417, 0.9885)21.27η20.060.060.0020.059790.01391(0.03739, 0.09141)60.33Table 6
The estimated result for MM model using simulated datasets.
TruePriorPosteriorCredible intervalsInefficiencyMeanStdMeanStdΓ−0.250.003.162−0.29350.3218(−0.9599, 0.3082)1.309Δ0.970.970.1000.96860.01209(0.9404, 0.9881)21.49η20.060.060.0020.059540.01438(0.03718, 0.09316)61.27Β0.700.070.1000.70090.07026(0.5629, 0.8383)33.13κ20.0050.0051.667e −30.0051380.001748(0.002767, 0.009502)148.0To measure how well the chain mixes, we calculate the inefficiency factors (the inverse of inefficiency factor is also known as numerical efficiency in Geweke [59]), which are given in Tables 5 and 6. The inefficient factor (equal to the autocorrelation time) is defined as(75)factor=1+2∑l=1Lρl,where ρl is the sample autocorrelation at lag l, which is calculated from the sampled values. L is chosen to be large enough so that the autocorrelation tapers off, which is calculated based on an autocorrelation rule of thumb. The inefficiency factor is also the ratio of the numerical variance of the posterior sample mean to the variance of the posterior sample mean from the hypothetical uncorrelated draws. Thus, it suggests the relative number of correlated draws necessary to attain the same variance of the posterior sample mean from the uncorrelated draws. In the ideal case where the posterior draws are independent, the corresponding inefficiency factor is 1. More generally, the inefficiency factor measures how many extra draws are needed to give results equivalent to this ideal case. For example, an inefficiency factor of 100 indicates that roughly 10000 posterior draws are required to give the same information as 100 independent draws.


## 6. Empirical Applications

### 6.1. The Dataset and Hyperparameters for the Parameters

In this section, we consider applications of the proposed method to a dataset that has been extensively analyzed in the finance literature. The data series for the study consists of daily continuously compounded returns,rt=100logpt/pt−1, in percentage, on the Standard and Poor’s 500 (S&P 500) index (computed without considering dividends) from January 2009 to December 2014, with a total of T=1509 observations. A plot of the data is given in Figure 2. The S&P 500 index time series is stable, and there is no obvious market jump which is the problem we need to study in the future. The S&P 500 index is a popular benchmark index of large-cap stocks in the United States, which is a price index, meaning it represents the stock prices of the companies within the index.Figure 2
The S&P 500 index from January 2009 to December 2014.We use the priors given in Section3.1 and set the same hyperparameters for the parameters that are common across models as follows: μ0=0, Vμ=10, ϕ0=0, Vϕ=10, γ0=0, Vγ=10, Vδ=0.01, vσ=11, α0=0, Vα=10, Vβ=0.01, vϕ=11, Sϕ=0.05, υ¯=2, and υ¯=100. The other hyperparameters for the parameters are as follows: δ0=0.97 for the SV model, δ0=0.975 for the SVM model, the MM model, and the SVt model, δ0=0.98 for the SVMt model and the MMt model, Sσ=0.65 for the SV model, Sσ=0.5 for the SVM model and the MM model, Sσ=0.45 for the SVt model, Sσ=0.35 for the SVMt model and the MMt model, β0=0.65 for the MM model, and β0=0.80 for the MMt model.


### 6.2. Comparison Results of Methods

To demonstrate the computation performance of the precision-based simulation smoother, a Kalman filter-based MCMC method for SV model and SVM model is also developed. The Kalman filter-based simulation smoother algorithm and the auxiliary mixture sampler (see, e.g., Kim et al. [19]) method are compared with our precision-based simulation smoother method on the SV model. For the SVM model, the details of the Kalman filter-based smoother and simulation smoother algorithm can be seen in Appendix. Note that all procedures are implemented in Matlab R2014b on Intel 4 Core i7-4790 @ 3.60 GHz CPU.In Figure3, we report the inefficiency factors corresponding to the posterior draws of σ=σ1,σ2,…,σTT in the SV model and the SVM model using the precision-based simulation smoother method and Kalman filter-based simulation smoother method, respectively. Note that each vector is of length T, so we have a total of 2T inefficiency factors. To present the information visually, boxplots are reported in Figure 3, where the middle line of the box denotes the median and the box denotes the median, while the lower and upper lines represent the 25 and the 75 percentiles, respectively. The whiskers extend to the maximum and minimum. For example, the boxplot associated with σ indicates that about 50% of the log-volatilities of the precision-based method have inefficiency factors between 1 and 2, which is significantly low compared to the Kalman filter-based method.Figure 3
The inefficiency factor boxplots of the posterior drawsσ in the SV model and the SVM model.Figures4 and 5 report sample paths of the parameters in the SV model and the SVM model of the precision-based simulation smoother method and the Kalman filter-based simulation smoother method, respectively, and the results are similar and support that two methods are practically stable. The estimation results of the Kalman filter-based MCMC method for the SV model and the SVM model are reported in Table 3, which are also similar to the estimation results of the precision-based method (see Table 4). Last but most important, the MCMC estimation times are 61.09 s, 1583 s, 62.76 s, and 4517 s, using the precision-based method for the SV model, the Kalman filter-based method for the SV model, the precision-based method for the SVM model, and the Kalman filter-based method for the SVM model, respectively, which are all of length 10000 after a burn-in period of 10000. These results are inspiring, since the computational cost of the proposed precision-based method is only 3.86% or 1.39% of that of the Kalman filter-based method for the SV model or the SVM model. This is consistent with the analyses in McCausland et al. [25].Figure 4
The sample paths of the parameters in the SV model.Figure 5
The sample paths of the parameters in the SVM model.


### 6.3. The Estimation Results

The evolution ofσt and ϕt for MMt models (19)–(22) estimated by the proposed method is depicted in Figure 6. It is clear that the estimates of surplus demand ϕt in the MMt model are much smoother than the return rt itself in Figure 2, so we may take the evolution of ϕt as the moving trend of the return rt. On the other hand, the MM model can capture the effect of the dynamics offset between supply and demand, and it can capture the “surplus demand” and obtain the surplus yields paid over the risk return. In this way, using the market microstructure model, one may use a more stable and predictable surplus demand process rather than the strongly random return process itself to assist the trading decision making in the financial market. It is therefore suggested that one may allocate assets by monitoring and making use of the surplus demand process. This methodology was implemented in, e.g., Peng et al. [5], Peng et al. [6], Peng et al. [7], and Qin et al. [9].Figure 6
The evolution of the estimated posterior meansσt and ϕt.Tables4 and 7 show the estimated parameters posterior means, standard deviations, inefficiency factor, and the 95% credible intervals estimates for all six models based on 10000 posterior draws after a burn-in period of 10000. We also present the results of two diagnostic tests. Specifically, we report the Ljung–Box and McLeod–Li statistics of order 20 computed on the standardized residuals and squared standardized residuals, respectively. All the Ljung–Box tests fail to reject the null hypothesis of no serial correlation in the standardized residuals at the 5% level. This diagnostic result also supports the SVM model and the MM model, since the Ljung–Box statistics value of them is significantly lower than that of the SV model. In addition, all the McLeod–Li tests also fail to reject the null hypothesis of no serial correlation in the squared standardized residuals at the 5% level, and it indicates that the models adequately capture the time-varying volatility of the data. All in all, the Ljung–Box and McLeod–Li statistics results indicate that the estimations are successful.Table 7
Parameter posterior means, standard deviations (in parentheses), inefficiency factors (in braces), and 95% intervals (in square brackets) for the stochastic volatility models using the precision-based simulation smoother method.
SVtSVMtMMtμ0.1054 (0.01852) {1.913} [0.06918, 0.1410]ϕ0.06721 (0.02995) {8.675} [0.009542, 0.1287]γ−0.5033 (0.2784) {6.866} [−1.038, 0.05783]−0.4352 (0.3285) {3.655} [−1.045, 0.2024]−0.4248 (0.3211) {4.445} [−1.020, 0.2520]δ0.9764 (0.007722) {30.70} [0.9602, 0.9903]0.9816 (0.006569) {20.07} [0.9676, 0.9935]0.9817 (0.006574) {26.02} [0.9677, 0.9935]η20.04615 (0.01081) {81.64} [0.02876, 0.06980]0.03495 (0.007486) {69.03} [0.02252, 0.05221]0.03429 (0.007483) {76.51} [0.02308, 0.05213]α0.06715 (0.03168) {27.18} [0.007429, 0.1329]β0.7885 (0.06333) {33.93} [0.6542, 0.8985]κ20.005679 (0.001688) {110.5} [0.003102, 0.009628]υ7.998 (2.378) {143.6} [4.900, 14.19]8.889 (2.634) {124.1} [5.465, 15.05]9.551 (2.988) {123.0} [5.742, 17.22]ℚLB16.58 (2.877)12.71(2.577)14.99 (3.103)ℚML18.56 (5.803)19.23 (5.977)18.50 (5.775)Note. ℚLB and ℚML are, respectively, the Ljung–Box and McLeod–Li statistics of order 20 computed on the standard residuals and squared standard residuals, and the numerical standard mistakes are in parentheses. The 0.05 and 0.01 critical values are 31.41 and 37.57, respectively.


### 6.4. The Prediction Results

In this part, a recursive out-of-sample forecasting exercise is performed to calculate the forecasting performance of the models listed in Table1 for the S&P 500 index returns at various horizons. Every one of the six models is used to make both point and density k-step-ahead iterated forecasts with k=1,2,5. Especially given the data up to time t, symbolized as r1:t=r1,r2,…,rtT, the MCMC sampler presented in Sections 3–4 is implemented to get the posterior draws given r1:t. Then, we calculate the predictive mean Ert+k|r1:t as the point forecast and the predictive density prt+k|r1:t as the density forecast. Next, we move one period in advance with building samples and repeat the whole exercise with data r1:t+1, and so on. The entire sample beginning with the start of the asset return series is utilized to generate samples, and it is assumed that all stochastic volatility models are the stable models. These forecasts are then evaluated for t=t0,…,T−k, where t0 is January 2, 2014. An autoregressive model as standard benchmark is also considered for comparison. It is the following special case of the vector autoregressive model (VAR) which is composed of a system of more than one stochastic difference equation.(76)rt=μ+βrt−1+σεtr,where μ, β, and σ are the constants to be estimated.In practice, neither the predictive mean nor the predictive density ofrt+k can be calculated analytically. Instead, they are obtained using predictive simulation. More precisely, at every MCMC iteration, given the model parameters and states up to time t, one can simulate future states from time t+1 to t+k using the relevant transition equations. The future errors eσs/2εsr for the SV type model or eσs/2λs1/2εsr for the SVt type model are also simulated where s=t+1,…,t+k. Given these draws, rt+k is a normal random variable as specified in (1), and one can easily produce the point and density forecasts for rt+k. Hence, one has a pair of forecasts (point and density) at every MCMC iteration. These forecasts are then averaged over all the posterior draws to produce estimates for Ert+k|r1:t and prt+k|r1:t. The whole exercise is then repeated using the data up to timet+1 to produce Ert+k+1|r1:t+1 and prt+k+1|r1:t+1, and so on.Letrt+ko symbolize the observed value of rt+k which is known at time t+k. The metric that is used to calculate the point forecasts is the root mean squared forecast mistake.RMSFE is defined as(77)RMSFE=∑t=t0T−krt+ko−Ert+k|r1:t2T−k−t0+1.To evaluate the density forecastprt+k|r1:t, one natural method is the predictive likelihood prt+k=rt+ko|r1:t, i.e., the predictive density of rt+k calculated at the observed value rt+ko. Obviously, if the actual result rt+ko is impossibly under the density forecast, the value of the predictive likelihood can be small, and vice versa (see, e.g., Geweke and Amisano [60] for a more detailed discussion of the predictive likelihood and its combination with the marginal likelihood). The density forecasts are evaluated by using the average log predictive likelihood scores (LPSs) as follows:(78)LPS=∑t=t0T−klogprt+k=rt+ko|r1:tT−k−t0+1.For this metric, a larger value shows prediction performance.Table8 gives the point forecast and the density prediction results for seven models. One can find overwhelming support for all six stochastic volatility models compared to the VAR model from the one-step and two-step LPS, since the time-varying volatility is more flexible than const volatility. One can also find overwhelming support for the MM model and SVM model compared to standard SV model with their RMSFE values and LPS, since the volatility feedback effect is essential feature for financial time series. In addition, the stochastic volatility models plus Student-t heavy tails model may significantly improve the density forecast and the point forecast accuracy as seen in Table 8.Table 8
The predictive mean and the predictive density of stochastic volatility models.
k = 1k = 2k = 5RMSFELPSRMSFELPSRMSFELPSVAR0.71801−1.272680.71803−1.272860.71812−1.27338SV0.71860−1.156410.71852−1.186260.71846−1.27845SVt0.71860−0.978350.71851−0.992800.71853−1.03235SVM0.71435−1.137770.71422−1.159350.71390−1.23931SVMt0.71377−0.994230.71360−1.000310.71311−1.03287MM0.71428−1.139870.71355−1.157570.71327−1.24026MMt0.71363−1.002820.71355−1.014300.71372−1.05082


## 6.1. The Dataset and Hyperparameters for the Parameters

In this section, we consider applications of the proposed method to a dataset that has been extensively analyzed in the finance literature. The data series for the study consists of daily continuously compounded returns,rt=100logpt/pt−1, in percentage, on the Standard and Poor’s 500 (S&P 500) index (computed without considering dividends) from January 2009 to December 2014, with a total of T=1509 observations. A plot of the data is given in Figure 2. The S&P 500 index time series is stable, and there is no obvious market jump which is the problem we need to study in the future. The S&P 500 index is a popular benchmark index of large-cap stocks in the United States, which is a price index, meaning it represents the stock prices of the companies within the index.Figure 2
The S&P 500 index from January 2009 to December 2014.We use the priors given in Section3.1 and set the same hyperparameters for the parameters that are common across models as follows: μ0=0, Vμ=10, ϕ0=0, Vϕ=10, γ0=0, Vγ=10, Vδ=0.01, vσ=11, α0=0, Vα=10, Vβ=0.01, vϕ=11, Sϕ=0.05, υ¯=2, and υ¯=100. The other hyperparameters for the parameters are as follows: δ0=0.97 for the SV model, δ0=0.975 for the SVM model, the MM model, and the SVt model, δ0=0.98 for the SVMt model and the MMt model, Sσ=0.65 for the SV model, Sσ=0.5 for the SVM model and the MM model, Sσ=0.45 for the SVt model, Sσ=0.35 for the SVMt model and the MMt model, β0=0.65 for the MM model, and β0=0.80 for the MMt model.


## 6.2. Comparison Results of Methods

To demonstrate the computation performance of the precision-based simulation smoother, a Kalman filter-based MCMC method for SV model and SVM model is also developed. The Kalman filter-based simulation smoother algorithm and the auxiliary mixture sampler (see, e.g., Kim et al. [19]) method are compared with our precision-based simulation smoother method on the SV model. For the SVM model, the details of the Kalman filter-based smoother and simulation smoother algorithm can be seen in Appendix. Note that all procedures are implemented in Matlab R2014b on Intel 4 Core i7-4790 @ 3.60 GHz CPU.In Figure3, we report the inefficiency factors corresponding to the posterior draws of σ=σ1,σ2,…,σTT in the SV model and the SVM model using the precision-based simulation smoother method and Kalman filter-based simulation smoother method, respectively. Note that each vector is of length T, so we have a total of 2T inefficiency factors. To present the information visually, boxplots are reported in Figure 3, where the middle line of the box denotes the median and the box denotes the median, while the lower and upper lines represent the 25 and the 75 percentiles, respectively. The whiskers extend to the maximum and minimum. For example, the boxplot associated with σ indicates that about 50% of the log-volatilities of the precision-based method have inefficiency factors between 1 and 2, which is significantly low compared to the Kalman filter-based method.Figure 3
The inefficiency factor boxplots of the posterior drawsσ in the SV model and the SVM model.Figures4 and 5 report sample paths of the parameters in the SV model and the SVM model of the precision-based simulation smoother method and the Kalman filter-based simulation smoother method, respectively, and the results are similar and support that two methods are practically stable. The estimation results of the Kalman filter-based MCMC method for the SV model and the SVM model are reported in Table 3, which are also similar to the estimation results of the precision-based method (see Table 4). Last but most important, the MCMC estimation times are 61.09 s, 1583 s, 62.76 s, and 4517 s, using the precision-based method for the SV model, the Kalman filter-based method for the SV model, the precision-based method for the SVM model, and the Kalman filter-based method for the SVM model, respectively, which are all of length 10000 after a burn-in period of 10000. These results are inspiring, since the computational cost of the proposed precision-based method is only 3.86% or 1.39% of that of the Kalman filter-based method for the SV model or the SVM model. This is consistent with the analyses in McCausland et al. [25].Figure 4
The sample paths of the parameters in the SV model.Figure 5
The sample paths of the parameters in the SVM model.


## 6.3. The Estimation Results

The evolution ofσt and ϕt for MMt models (19)–(22) estimated by the proposed method is depicted in Figure 6. It is clear that the estimates of surplus demand ϕt in the MMt model are much smoother than the return rt itself in Figure 2, so we may take the evolution of ϕt as the moving trend of the return rt. On the other hand, the MM model can capture the effect of the dynamics offset between supply and demand, and it can capture the “surplus demand” and obtain the surplus yields paid over the risk return. In this way, using the market microstructure model, one may use a more stable and predictable surplus demand process rather than the strongly random return process itself to assist the trading decision making in the financial market. It is therefore suggested that one may allocate assets by monitoring and making use of the surplus demand process. This methodology was implemented in, e.g., Peng et al. [5], Peng et al. [6], Peng et al. [7], and Qin et al. [9].Figure 6
The evolution of the estimated posterior meansσt and ϕt.Tables4 and 7 show the estimated parameters posterior means, standard deviations, inefficiency factor, and the 95% credible intervals estimates for all six models based on 10000 posterior draws after a burn-in period of 10000. We also present the results of two diagnostic tests. Specifically, we report the Ljung–Box and McLeod–Li statistics of order 20 computed on the standardized residuals and squared standardized residuals, respectively. All the Ljung–Box tests fail to reject the null hypothesis of no serial correlation in the standardized residuals at the 5% level. This diagnostic result also supports the SVM model and the MM model, since the Ljung–Box statistics value of them is significantly lower than that of the SV model. In addition, all the McLeod–Li tests also fail to reject the null hypothesis of no serial correlation in the squared standardized residuals at the 5% level, and it indicates that the models adequately capture the time-varying volatility of the data. All in all, the Ljung–Box and McLeod–Li statistics results indicate that the estimations are successful.Table 7
Parameter posterior means, standard deviations (in parentheses), inefficiency factors (in braces), and 95% intervals (in square brackets) for the stochastic volatility models using the precision-based simulation smoother method.
SVtSVMtMMtμ0.1054 (0.01852) {1.913} [0.06918, 0.1410]ϕ0.06721 (0.02995) {8.675} [0.009542, 0.1287]γ−0.5033 (0.2784) {6.866} [−1.038, 0.05783]−0.4352 (0.3285) {3.655} [−1.045, 0.2024]−0.4248 (0.3211) {4.445} [−1.020, 0.2520]δ0.9764 (0.007722) {30.70} [0.9602, 0.9903]0.9816 (0.006569) {20.07} [0.9676, 0.9935]0.9817 (0.006574) {26.02} [0.9677, 0.9935]η20.04615 (0.01081) {81.64} [0.02876, 0.06980]0.03495 (0.007486) {69.03} [0.02252, 0.05221]0.03429 (0.007483) {76.51} [0.02308, 0.05213]α0.06715 (0.03168) {27.18} [0.007429, 0.1329]β0.7885 (0.06333) {33.93} [0.6542, 0.8985]κ20.005679 (0.001688) {110.5} [0.003102, 0.009628]υ7.998 (2.378) {143.6} [4.900, 14.19]8.889 (2.634) {124.1} [5.465, 15.05]9.551 (2.988) {123.0} [5.742, 17.22]ℚLB16.58 (2.877)12.71(2.577)14.99 (3.103)ℚML18.56 (5.803)19.23 (5.977)18.50 (5.775)Note. ℚLB and ℚML are, respectively, the Ljung–Box and McLeod–Li statistics of order 20 computed on the standard residuals and squared standard residuals, and the numerical standard mistakes are in parentheses. The 0.05 and 0.01 critical values are 31.41 and 37.57, respectively.


## 6.4. The Prediction Results

In this part, a recursive out-of-sample forecasting exercise is performed to calculate the forecasting performance of the models listed in Table1 for the S&P 500 index returns at various horizons. Every one of the six models is used to make both point and density k-step-ahead iterated forecasts with k=1,2,5. Especially given the data up to time t, symbolized as r1:t=r1,r2,…,rtT, the MCMC sampler presented in Sections 3–4 is implemented to get the posterior draws given r1:t. Then, we calculate the predictive mean Ert+k|r1:t as the point forecast and the predictive density prt+k|r1:t as the density forecast. Next, we move one period in advance with building samples and repeat the whole exercise with data r1:t+1, and so on. The entire sample beginning with the start of the asset return series is utilized to generate samples, and it is assumed that all stochastic volatility models are the stable models. These forecasts are then evaluated for t=t0,…,T−k, where t0 is January 2, 2014. An autoregressive model as standard benchmark is also considered for comparison. It is the following special case of the vector autoregressive model (VAR) which is composed of a system of more than one stochastic difference equation.(76)rt=μ+βrt−1+σεtr,where μ, β, and σ are the constants to be estimated.In practice, neither the predictive mean nor the predictive density ofrt+k can be calculated analytically. Instead, they are obtained using predictive simulation. More precisely, at every MCMC iteration, given the model parameters and states up to time t, one can simulate future states from time t+1 to t+k using the relevant transition equations. The future errors eσs/2εsr for the SV type model or eσs/2λs1/2εsr for the SVt type model are also simulated where s=t+1,…,t+k. Given these draws, rt+k is a normal random variable as specified in (1), and one can easily produce the point and density forecasts for rt+k. Hence, one has a pair of forecasts (point and density) at every MCMC iteration. These forecasts are then averaged over all the posterior draws to produce estimates for Ert+k|r1:t and prt+k|r1:t. The whole exercise is then repeated using the data up to timet+1 to produce Ert+k+1|r1:t+1 and prt+k+1|r1:t+1, and so on.Letrt+ko symbolize the observed value of rt+k which is known at time t+k. The metric that is used to calculate the point forecasts is the root mean squared forecast mistake.RMSFE is defined as(77)RMSFE=∑t=t0T−krt+ko−Ert+k|r1:t2T−k−t0+1.To evaluate the density forecastprt+k|r1:t, one natural method is the predictive likelihood prt+k=rt+ko|r1:t, i.e., the predictive density of rt+k calculated at the observed value rt+ko. Obviously, if the actual result rt+ko is impossibly under the density forecast, the value of the predictive likelihood can be small, and vice versa (see, e.g., Geweke and Amisano [60] for a more detailed discussion of the predictive likelihood and its combination with the marginal likelihood). The density forecasts are evaluated by using the average log predictive likelihood scores (LPSs) as follows:(78)LPS=∑t=t0T−klogprt+k=rt+ko|r1:tT−k−t0+1.For this metric, a larger value shows prediction performance.Table8 gives the point forecast and the density prediction results for seven models. One can find overwhelming support for all six stochastic volatility models compared to the VAR model from the one-step and two-step LPS, since the time-varying volatility is more flexible than const volatility. One can also find overwhelming support for the MM model and SVM model compared to standard SV model with their RMSFE values and LPS, since the volatility feedback effect is essential feature for financial time series. In addition, the stochastic volatility models plus Student-t heavy tails model may significantly improve the density forecast and the point forecast accuracy as seen in Table 8.Table 8
The predictive mean and the predictive density of stochastic volatility models.
k = 1k = 2k = 5RMSFELPSRMSFELPSRMSFELPSVAR0.71801−1.272680.71803−1.272860.71812−1.27338SV0.71860−1.156410.71852−1.186260.71846−1.27845SVt0.71860−0.978350.71851−0.992800.71853−1.03235SVM0.71435−1.137770.71422−1.159350.71390−1.23931SVMt0.71377−0.994230.71360−1.000310.71311−1.03287MM0.71428−1.139870.71355−1.157570.71327−1.24026MMt0.71363−1.002820.71355−1.014300.71372−1.05082


## 7. Conclusions and Future Research

In this paper, we studied the more general MM and MMt models compared to the MM type models researched in previous works and proposed a new approach for efficiently estimating the general MM type models, which are non-linear non-Gaussian models. Because of the general applicability of the proposed approach, it proved its usefulness in a wide range of applications for stochastic volatility models. The precision-based efficient simulation smoother algorithm and the efficient acceptance-rejection Metropolis–Hastings procedure were extended to develop a fast and efficient MCMC method for six SV models. The presented methodology was demonstrated with a simulation study and the inefficiency factors and posterior densities of parameters and states. The model estimation results showed substantial highlight for the empirical relevance of the proposed extension. In order to demonstrate computation performance of the precision-based simulation smoother, a Kalman filter-based MCMC method for SVM model was also developed. Furthermore, in a recursive forecasting exercise, the market microstructure model plus Student-t heavy tails performed better than various standard benchmarks in both point forecast and density forecast.For future research, it would be interesting to investigate if the proposed market microstructure model variant (plus jump, leverage, and intertemporal volatility feedback) also fits other macroeconomic or financial time series better than the standard market microstructure model. Since it is difficult to obtain an accurate estimate of marginal likelihood and deviance information criterion (DIC) through existing algorithms for the market microstructure model, due to the marginal likelihood evaluation involving “integrating out” both types of states, i.e.,ϕt and σt, we will consider an adaptive importance sampling approach to estimate the marginal likelihood and DIC in the future, where the draws are typically costly to obtain and highly correlated in high-dimensional settings. Furthermore, in the recursive forecasting exercise, since the market microstructure model performed better than various standard benchmarks in both point forecast and density forecast, it would be interesting to do optimization portfolio and asset allocation using the market microstructure model and study the relationship between the excess demand ϕt and the trend of return rt.



---

*Source: 5523468-2021-10-27.xml*