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.
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 Section 2, 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 to date is , where is the asset price. The one-period simple compound return is . For convenience, the log return is usually used in modeling; especially in case study, one often uses 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:where and 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, and the logarithmic volatility is initialized with , where 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:where the regression coefficient measures the volatility-in-mean effect or volatility feedback effect; when , 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:where is the log price of an asset, is the market surplus demand, and is the market liquidity. At any given instant of time, there is an instantaneous demand and an instantaneous supply for the asset. The dynamical equations (5) and (6) describe the impact of supply-demand balance, i.e., tends to push the price up (if ) or down (if ) 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 is driven by the surplus demand , and the range of price changes depends on the market liquidity .
However, in the original setting, models (5) and (6) only provide an abstract description to the dynamics of stock market because the surplus demand and the market liquidity cannot be easily measured directly. Iino and Ozaki [4] treated the market liquidity and the surplus demand 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:where , , and 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:where , , and , , and 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 and , respectively. Equations (9) and (10) describe the unobservable surplus demand and market liquidity 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 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:where are conditionally independent, given the model parameters (including ) and data, and 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:and the MM plus Student-t heavy tail (MMt) model:
All six models studied in this paper are summarized in Table 1. 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.
For 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 . However, this approach cannot be easily generalized to the time-varying surplus demand in MM models (10)–(12) and MMt models (19)–(22) due to the likelihood evaluation in the case that involves two kinds of status, and . Meanwhile, applying the Kalman filter-based smoother to obtain a smoothed estimate for 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 Table 1. 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 , rather than directly, information about the sign of the return is lost. Therefore, the distribution of or 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 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 given the return and the excess demand 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 . Using forward- and back-substitution to simulate smoother market liquidity , the Cholesky decomposition diagonal band matrix Hessian can be obtained in operations instead of operations for full matrices. Meanwhile, applying the precision-based smoother to obtain a smoothed estimate for market liquidity leads to a very simple and fast Newton–Raphson iterative procedure, since the Hessian of the log-conditional density of the market liquidity 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 appears in both mean term and variance term of return equation (17) as volatility and the model is non-linear in the market liquidity , one cannot directly sample the market liquidity 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, let , , , and .
Now, we discuss a quick way to get the Gaussian approximation density for the joint conditional density where . By Bayes’ theorem, one has , and the prior density that characterizes the state equation (19) is Gaussian density. By approximating the observation density , which characterizes return equation (17) via a Gaussian density in , one obtains a Gaussian approximation of at once. To this end, following equation (17), the log-density related to is given as follows:where is independent of . It is easy to examine that the first and second partial derivatives of the log-density with respect to are given by
As the first and second derivatives can be calculated quickly, we should note that when a point is given, one can estimate the log conditional likelihood by making use of a second-order Taylor expansion around to obtain
Let , which can be written in matrix form as follows:where is independent of , , , and . Next, rewrite state equation (19) in matrix form below:where , , and
By a change of variable, one haswith log-densitywhere and are independent of . Then, considering the expressions in (24) and (28), one haswhere is independent of . Since the expression in (29) is the log-kernel of in , where is the precision matrix or inverse of the variance and is the mean vector, can be calculated by the Gaussian density
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 of such as and obtain a random draw from efficiently such that and . A crucial result is that the Cholesky decomposition can be obtained in operations rather than 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 and can be quickly obtained by the numerical optimization. Using the Newton–Raphson method, one can quickly locate the maximum of , or equivalently that of , or to select the point as the mode of or smooth . It follows from (29) that negative Hessian of evaluated at is and the gradient at is . Hence, one can carry out the Newton–Raphson method as follows: initializing with for some lasting vector . For assume in the evaluation of and and compute
Note that since is Gaussian density, its mode equals its mean vector . Repeating this procedure until some convergence criterion is reached, e.g., we get the mode or when for some prefixed tolerance level . Also, note that the negative Hessian 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 of and we have the Gaussian approximation proposal density from which we generate candidate draws . In the classic accepting-rejecting sampling, an important requirement is that there exists a constant , which is selected by equation (34), such thatfor in support of . 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 setand let denote the complement code. Then, the ARMH procedure is as follows.
| ||||||||||||||||||
Adjusting the original Gaussian proposal density by the accept-reject step, a better approximation of the target density is achieved, i.e., the new proposed density is identical with the target density on the set (albeit different normalizing constants ), whereas on , the new proposed density is decreased to the original one.
Chib and Jeliazkov [46] presented a practical way to select the constant and the trade-off in such a choice, which is outlined here. Notice that if a bigger 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 setwhere is the mode of the conditional density and is, say, between 1 and 5. Such kind of choice would make sure that is sufficiently small to decrease the required number of draws from , 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 , , , , , , and , where 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 , which are conditionally independently, to give the parameters and data. Each follows an independent inverse-gamma distribution:where .
Next, since and the prior for , one can first derive the log-density following equation (20) to sample :where 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 bywhere and 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 distribution and maximize using numerical optimization based on the approximation of 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 , respectively. Then, the independence-chain Metropolis–Hastings step with proposal distribution is implemented as follows:and the proposed draw is accepted with probability
4.2. Posterior Simulators for the Parameters 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 following equation (19) as follows.
Sampling is straightforward, if we assume the prior for . It yieldswhere , .
Sampling can be done easily because the conditional distribution is standard. Especially following and the prior for , one haswhere and are independent of . As the log-density is quadratic in , it is Gaussian and one haswherewhere is the precision matrix or inverse of variance and is the mean.
Since is non-standard in , the conditional density is non-standard, and one can carry out an independence-chain Metropolis–Hastings step to draw . Firstly, following and the prior for , one haswhere and are independent of . As the approximate log-density is quadratic in , it is approximate Gaussian and one haswhere is the precision matrix or inverse of variance and is the mean.
Then, using the approximate Gaussian density with as proposal density, one can implement an independence-chain Metropolis–Hastings in which the proposal is accepted in the feasible region with probabilitywherewhere is independent of .
4.3. Posterior Simulators for the State and the Parameters of MMt Models (19)–(22)
Since , , and are all Gaussian distributions in , one immediately obtains a Gaussian distribution . Now, write the return equation (17) and exceed demand equation (18) of in matrix formwhere , , , , and
Then, and with log-densitywhere and are independent of . Combining (49) and (50), one haswhere is independent of . As this log-density is quadratic in , it is Gaussian and one haswhere is the precision matrix or inverse of variance and is the mean.
Lastly, following equations (18) and (19), , and are sampled the same as those of , and , respectively.
For , it has inverse-gamma posterior densitywhere , .
For , one gets the Gaussian posterior densitywherewhere is the precision matrix or inverse of variance and is the mean.
For , it has approximate Gaussian posterior densitywhere is the precision matrix or inverse of variance and is the mean.
Then, one implements an independence-chain Metropolis–Hastings in which the proposal is accepted in the feasible region with probabilitywhere .
Now an overall description of Algorithm 2 is given below.
| ||||||||||||||||||||||||||||||
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 Algorithm 3.
| ||||||||||||||||
In Step 2 above, sampling can be done easily, as the conditional distribution is standard. Since the prior for and , following (1), one hasandwhere , , and are independent of . As this log-density is quadratic in , it is Gaussian and one haswhere is the precision matrix or inverse of variance and is the mean.
In Step 3(a), the joint conditional density is similar to that of the MMt model in which the log-density follows equation (1) instead of equation (17).where is independent of . Then, starting from equation (61) instead of equation (21), 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 Algorithm 4.
| ||||||||||||||||
In Step 2, Sampling can be done easily, as the conditional distribution is standard. Since the prior for and , following (3), one hasandwhere , , and are independent of . As the log-density is quadratic in , it is Gaussian and one haswhere is the precision matrix or inverse of variance and is the mean.
In Step 3(a), the joint conditional density is similar to that of the MMt model in which the log-density follows equation (6) instead of equation (17).where is independent of . Then, starting from equation (49) instead of equation (21), 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 Algorithm 5.
| ||||||||||||||||||||||||
In Step 2(a), the joint conditional density is similar to that of the MMt model in which follows equation (8) instead of equation (17). One can write return equation (8) in matrix form:where . 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 density is similar to that of the MMt model in which the log-density follows equation (8) instead of equation (17).where is independent of . Then, starting from equation (67) instead of equation (21), 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 Algorithm 6.
| ||||||||||||||||||||||
In Step 2, sampling is similar to that of the MMt model in which the log-density follows equation (11) instead of equation (1).where 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 density is similar to that of the MMt model in which the log-density follows equation (11) instead of equation (17).where is independent of . Then, starting from equation (69) instead of equation (21), can be sampled as that of the MMt model.
In Step 4(a), to sample is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:where .
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 Algorithm 7.
| ||||||||||||||||||||||
In Step 2, Sampling is similar to that of the SVM model in which the log-density follows equation (14) instead of equation (3).where 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 density is similar to that of the MMt model in which the log-density follows equation (14) instead of equation (17).where is independent of . Then, starting from equation (73) instead of equation (19), can be sampled as that of the MMt model.
In Step 4(a), to sample is similar to that of the MMt model, and it follows an independent inverse-gamma distribution:where .
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 Table 1, and the difference between previous research and our research in this paper is briefly explained in Table 2.
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 of 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 for the SV model and for the MM model. Parameters for the log-volatility or market liquidity transition are made to be , , and for both the SV model and the MM model. Parameters for the surplus demand transition are made to be and 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.
The parameter prior distributions studied in Section 4 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 , , , , and , we set , , , , , , , , , and . 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 , for the hyperparameters , we set the mean and a small variance for making the MCMC convergence fast.
Because the regression coefficient of the built MM models (10)–(12) is set as in this simulation study, for the hyperparameters , we set the mean and a small variance 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 in this simulation study, for the hyperparameters , , we set the mean and a simple large variance since the mean of the volatility process, , 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, and , in financial economic viewpoint for investment arbitrage, the smaller the volatility of the state, the better, so we select sufficient simple small mean and variance for the , and thus we set and for the built SV models (1) and (2) and MM models (10)–(12) that are used to generate the observation data for simulation study.
Figure 1 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 .
Tables 5 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.
To 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 aswhere is the sample autocorrelation at lag , which is calculated from the sampled values. 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, , 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 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.
We use the priors given in Section 3.1 and set the same hyperparameters for the parameters that are common across models as follows: , , , , , , , , , , , , , , and . The other hyperparameters for the parameters are as follows: for the SV model, for the SVM model, the MM model, and the SVt model, for the SVMt model and the MMt model, for the SV model, for the SVM model and the MM model, for the SVt model, for the SVMt model and the MMt model, for the MM model, and 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 Figure 3, we report the inefficiency factors corresponding to the posterior draws of 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 , so we have a total of 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.
Figures 4 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].
6.3. The Estimation Results
The evolution of and for MMt models (19)–(22) estimated by the proposed method is depicted in Figure 6. It is clear that the estimates of surplus demand in the MMt model are much smoother than the return itself in Figure 2, so we may take the evolution of as the moving trend of the return . 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].
Tables 4 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.
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 Table 1 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 . Especially given the data up to time , symbolized as , the MCMC sampler presented in Sections 3–4 is implemented to get the posterior draws given . Then, we calculate the predictive mean as the point forecast and the predictive density as the density forecast. Next, we move one period in advance with building samples and repeat the whole exercise with data , 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 , where 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.where , , and are the constants to be estimated.
In practice, neither the predictive mean nor the predictive density of 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 , one can simulate future states from time to using the relevant transition equations. The future errors for the SV type model or for the SVt type model are also simulated where . Given these draws, is a normal random variable as specified in (1), and one can easily produce the point and density forecasts for . 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 and . The whole exercise is then repeated using the data up to time to produce and , and so on.
Let symbolize the observed value of which is known at time . The metric that is used to calculate the point forecasts is the root mean squared forecast mistake.
RMSFE is defined as
To evaluate the density forecast , one natural method is the predictive likelihood , i.e., the predictive density of calculated at the observed value . Obviously, if the actual result 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:
For this metric, a larger value shows prediction performance.
Table 8 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.
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., and , 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 and the trend of return .
Appendix
In this section, we present an augmented Kalman filter-based smoother and simulation smoother and consider a common state-space model:where and are uncorrelated, is fixed, has full column rank, and is non-singular unless . , , , , , and are the system matrices. The state-space models (A1) and (A2) are considered as a natural set for wide time series models to calculate parameters and to make prediction of future observations: the regression influence in time series models is put in ; we can find autoregressive moving average dynamics in and ; a model of additive outlier and structural change interventions can be constructed through using and ; unobservable elements are put in the state vector ; and the covariance construction is up to the matrices and .
There is a stationary initial condition, and , since . For notational convenience, let .
The Kalman filter is a forward recursion and run for .
, , , , , , .
and are the predictor and the mean squared error of , respectively, is the error of predictor , is covariance of , and is Kalman gain.
For state simulation smoother, set and , and the Kalman smoother is a backward recursion and run for .
, , , , , , , .
The disturbance smoothing is .
The simulation smoother for states is a forward recursion and run for .which is started off with .
For state smoother, set and , and the Kalman smoother is a backward recursion and run for .
The quick state smoother is a forward recursion and run for .which is started off with .
The detailed procedure may be seen, e.g., in De Jong [15], Koopman [16], De Jong and Shephard [17], Koopman et al. [61], and Durbin and Koopman [18].
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (71271215).