Estimation and Properties of a Time-Varying GQARCH(1,1)-M Model

Abstract

Time-varying GARCH-M models are commonly used in econometrics and financial economics. Yet the recursive nature of the conditional variance makes exact likelihood analysis of these models computationally infeasible. This paper outlines the issues and suggests to employ a Markov chain Monte Carlo algorithm which allows the calculation of a classical estimator via the simulated EM algorithm or a simulated Bayesian solution in only 𝑂(𝑇) computational operations, where 𝑇 is the sample size. Furthermore, the theoretical dynamic properties of a time-varying GQARCH(1,1)-M are derived. We discuss them and apply the suggested Bayesian estimation to three major stock markets.

1. Introduction

Time series data, emerging from diverse fields appear to possess time-varying second conditional moments. Furthermore, theoretical results seem to postulate quite often, specific relationships between the second and the first conditional moment. For instance, in the stock market context, the first conditional moment of stock market’s excess returns, given some information set, is a possibly time-varying, linear function of volatility (see, e.g., Merton [1], Glosten et al. [2]). These have led to modifications and extensions of the initial ARCH model of Engle [3] and its generalization by Bollerslev [4], giving rise to a plethora of dynamic heteroscedasticity models. These models have been employed extensively to capture the time variation in the conditional variance of economic series, in general, and of financial time series, in particular (see Bollerslev et al. [5] for a survey).

Although the vast majority of the research in conditional heteroscedasticity is being processed aiming at the stylized facts of financial stock returns and of economic time series in general, Arvanitis and Demos [6] have shown that a family of time-varying GARCH-M models can in fact be consistent with the sample characteristics of time series describing the temporal evolution of velocity changes of turbulent fluid and gas molecules. Despite the fact that the latter statistical characteristics match in a considerable degree their financial analogues (e.g., leptokurtosis, volatility clustering, and quasi long-range dependence in the squares are common), there are also significant differences in the behavior of the before mentioned physical systems as opposed to financial markets (examples are the anticorrelation effect and asymmetry of velocity changes in contrast to zero autocorrelation and the leverage effect of financial returns) (see Barndorff-Nielsen and Shephard [7] as well as Mantegna and Stanley [8, 9]). It was shown that the above-mentioned family of models can even create anticorrelation in the means as far as an AR(1) time-varying parameter is introduced.

It is clear that from an econometric viewpoint it is important to study how to efficiently estimate models with partially unobserved GARCH processes. In this context, our main contribution is to show how to employ the method proposed in Fiorentini et al. [10] to achieve MCMC likelihood-based estimation of a time-varying GARCH-M model by means of feasible 𝑂(𝑇) algorithms, where 𝑇 is the sample size. The crucial idea is to transform the GARCH model in a first-order Markov’s model. However, in our model, the error term enters the in-mean equation multiplicatively and not additively as it does in the latent factor models of Fiorentini et al. [10]. Thus, we show that their method applies to more complicated models, as well.

We prefer to employ a GQARCH specification for the conditional variance (Engle [3] and Sentana [11]) since it encompasses all the existing restricted quadratic variance functions (e.g., augmented ARCH model), its properties are very similar to those of GARCH models (e.g., stationarity conditions) but avoids some of their criticisms (e.g., very easy to generalize to multivariate models). Moreover, many theories in finance involve an explicit tradeoff between the risk and the expected returns. For that matter, we use an in-mean model which is ideally suited to handling such questions in a time series context where the conditional variance may be time varying. However, a number of studies question the existence of a positive mean/variance ratio directly challenging the mean/variance paradigm. In Glosten et al. [2] when they explicitly include the nominal risk free rate in the conditioning information set, they obtain a negative ARCH-M parameter. For the above, we allow the conditional variance to affect the mean with a possibly time varying coefficient which we assume for simplicity that it follows an AR(1) process. Thus, our model is a time-varying GQARCH-M-AR(1) model.

As we shall see in Section 2.1, this model is able to capture the, so-called, stylized facts of excess stock returns. These are (i) the sample mean is positive and much smaller than the standard deviation, that is, high coefficient of variation, (ii) the autocorrelation of excess returns is insignificant with a possible exception of the 1st one, (iii) the distribution of returns is nonnormal mainly due to excess kurtosis and may be asymmetry (negative), (iv) there is strong volatility clustering, that is, significant positive autocorrelation of squared returns even for high lags, and (v) the so-called leverage effect; that is, negative errors increase future volatility more than positive ones of the same size.

The structure of the paper is as follows. In Section 2, we present the model and derive the theoretical properties the GQARCH(1,1)-M-AR(1) model. Next, we review Bayesian and classical likelihood approaches to inference for the time-varying GQARCH-M model. We show that the key task (in both cases) is to be able to produce consistent simulators and that the estimation problem arises from the existence of two unobserved processes, causing exact likelihood-based estimations computationally infeasible. Hence, we demonstrate that the method proposed by Fiorentini et al. [10] is needed to achieve a first-order Markov’s transformation of the model and thus, reducing the computations from 𝑂(𝑇²) to 𝑂(𝑇). A comparison of the efficient (𝑂(𝑇) calculations) and the inefficient (𝑂(𝑇²) ones) simulator is also given. An illustrative empirical application on weekly returns from three major stock markets is presented in Section 4, and we conclude in Section 5.

2. GQARCH(1,1)-M-AR(1) Model

The definition of our model is as follows.

Definition 2.1. The time-varying parameter GQARCH(1,1)-M-AR(1) model: 𝑟_𝑡=𝛿_𝑡ℎ_𝑡+𝜀_𝑡,𝜀_𝑡=𝑧_𝑡ℎ^𝑡1/2,(2.1) where ℎ𝛿_𝑡=(1−𝜑)𝛿+𝜑𝛿_𝑡−1+𝜑_𝑢𝑢_𝑡,(2.2)_𝑡=𝜔+𝛼𝜀_𝑡−1−𝛾²+𝛽ℎ_𝑡−1,(2.3)𝑧_𝑡∽i.i.d.𝑁(0,1), 𝑢_𝑡∽i.i.d.𝑁(0,1), 𝑢_𝑡, 𝑧_𝑡 are independent for all𝑡^𝑠, and where {𝑟_𝑡}^𝑇𝑡=1 are the observed excess returns, 𝑇 is the sample size, {𝛿_𝑡}^𝑇𝑡=1 is an unobserved AR(1) process independent (with 𝛿₀=𝛿) of {𝜀_𝑡}^𝑇𝑡=1, and {ℎ_𝑡}^𝑇𝑡=1 is the conditional variance (with ℎ₀ equal to the unconditional variance and 𝜀₀=0) which is supposed to follow a GQARCH(1,1). It is obvious that 𝛿_𝑡 is the market price of risk (see, e.g., Merton [1] Glosten at al. [2]). Let us call ℱ_𝑡−1 the sequence of natural filtrations generated by the past values of {𝜀_𝑡} and {𝑟_𝑡}.

Modelling the theoretical properties of this model has been a quite important issue. Specifically, it would be interesting to investigate whether this model can accommodate the main stylized facts of the financial markets. On other hand, the estimation of the model requires its transformation into a first-order Markov’s model to implement the method of Fiorentini et al. [10]. Let us start with the theoretical properties.

2.1. Theoretical Properties

Let us consider first the moments of the conditional variance ℎ_𝑡, needed for the moments of 𝑟_𝑡. The proof of the following lemma is based on raising ℎ_𝑡 to the appropriate power, in (2.3), and taking into account that 𝐸(𝑧^4𝑡)=3, 𝐸(𝑧^6𝑡)=15 and 𝐸(𝑧^8𝑡)=105.

Lemma 2.2. If 105𝛼⁴+60𝛼³𝛽+18𝛼²𝛽²+4𝛼𝛽³+𝛽⁴<1, the first four moments of the conditional variance of (2.3) are given by 𝐸ℎ_𝑡=𝜔+𝛼𝛾²𝐸ℎ1−(𝛼+𝛽)_2𝑡,=𝜔+𝛼𝛾²𝜔+𝛼𝛾²(1+𝛼+𝛽)+4𝛼²𝛾²𝐸ℎ1−3𝛼²−2𝛼𝛽−𝛽²(1−𝛼−𝛽)_3𝑡,=𝜔+𝛼𝛾²³+3𝜔+𝛼𝛾²𝜔+𝛼𝛾²(𝛼+𝛽)+4𝛼²𝛾²𝐸ℎ_𝑡1−𝛽³−15𝛼³−9𝛼²𝛽−3𝛼𝛽²+3𝜔+𝛼𝛾²3𝛼²+𝛽²+2𝛽𝛼+4𝛼²𝛾²(3𝛼+𝛽)𝐸ℎ^2𝑡1−𝛽³−15𝛼³−9𝛼²𝛽−3𝛼𝛽²𝐸ℎ_4𝑡,=𝜔+𝛼𝛾²²𝜔+𝛼𝛾²²+4𝜔+𝛼𝛾²(𝛼+𝛽)+6𝛼²𝛾²𝐸ℎ_𝑡1−105𝛼⁴−60𝛼³𝛽−18𝛼²𝛽²−4𝛼𝛽³−𝛽⁴+6𝜔+𝛼𝛾²²3𝛼²+𝛽²+2𝛽𝛼+8𝜔+𝛼𝛾²(3𝛼+𝛽)+𝛼²𝛾²𝛼²𝛾²1−105𝛼⁴−60𝛼³𝛽−18𝛼²𝛽²−4𝛼𝛽³−𝛽⁴+4𝜔+𝛼𝛾₂𝐸ℎ^2𝑡15𝛼³+9𝛼²𝛽+3𝛼𝛽²+𝛽³+6𝛼²𝛾²15𝛼²+6𝛼𝛽+𝛽²1−105𝛼⁴−60𝛼³𝛽−18𝛼²𝛽²−4𝛼𝛽³−𝛽⁴𝐸ℎ^3𝑡.(2.4)

Consequently, the moments of 𝑟_𝑡 are given in the following theorem taken from Arvanitis and Demos [6].

Theorem 2.3. The first two moments of the model in (2.1), (2.2), and (2.3) are given by 𝐸𝑟_𝑡=𝛿𝐸ℎ_𝑡,𝐸𝑟^2𝑡=𝛿²+𝜑^2𝑢1−𝜑²𝐸ℎ^2𝑡+𝐸ℎ_𝑡,(2.5) whereas the skewness and kurtosis coefficients are Sk𝑟_𝑡=𝑆𝑟^𝑡Var^1.5𝑟_𝑡,kurt𝑟^𝑡=𝜅Var²𝑟_𝑡,(2.6) where 𝑆𝑟_𝑡=𝛿𝛿²+3𝜑^2𝑢1−𝜑²𝐸ℎ^3𝑡−3𝛿𝛿₂+3𝛿𝐸ℎ²+2𝛿³𝐸³ℎ_𝑡+𝜑^2𝑢1−𝜑²𝐸ℎ^2𝑡𝜅=𝛿₄+𝐸ℎ_𝑡𝐸ℎ_𝑡,+6𝛿²𝜑^2𝑢1−𝜑²+3𝜑^4𝑢1−𝜑²²𝐸ℎ^4𝑡+6𝛿₂+3𝛿²2−𝛿²𝐸ℎ_𝑡𝐸³ℎ_𝑡+𝜑^2𝑢1−𝜑²𝐸ℎ^3𝑡−4𝛿²𝛿²+3𝜑^2𝑢1−𝜑²𝐸ℎ^3𝑡+6𝛿₂𝐸ℎ_𝑡𝛿²+𝜑^2𝑢1−𝜑²𝐸(ℎ)−2𝐸(ℎ)+3𝐸ℎ^2𝑡,(2.7) and 𝐸(ℎ_𝑡), 𝐸(ℎ^2𝑡), 𝐸(ℎ^3𝑡),and 𝐸(ℎ^4𝑡) are given in Lemma 2.2.

In terms of stationarity, the process {𝑟_𝑡} is 4th-order stationary if and only if ||𝜑||<1,105𝛼⁴+60𝛼³𝛽+18𝛼²𝛽²+4𝛼𝛽³+𝛽⁴<1.(2.8) These conditions are the same as in Arvanitis and Demos [6], indicating that the presence of the asymmetry parameter, 𝛾, does not affect the stationarity conditions (see also Bollerslev [4] and Sentana [11]). Furthermore, the 4th-order stationarity is needed as we would like to measure the autocorrelation of the squared 𝑟_𝑡’s (volatility clustering), as well as the correlation of 𝑟^2𝑡 and 𝑟_𝑡−1 (leverage effect). The dynamic moments of the conditional variance and those between the conditional variance ℎ_𝑡 and the error 𝜀_𝑡 are given in the following two lemmas (for a proof see Appendix B).

Lemma 2.4. Under the assumption of Lemma 2.2, one has that CCoovvℎℎ_𝑡2𝑡,ℎ_𝑡−𝑘=(𝛼+𝛽)^𝑘𝑉ℎ_𝑡,,ℎ_𝑡−𝑘=𝐴^𝑘𝐸ℎ^3𝑡−𝐸ℎ^2𝑡𝐸ℎ_𝑡+(𝛼+𝛽)^𝑘−𝐴^𝑘Covℎ_𝑡(𝛼+𝛽)−𝐴𝐵𝑉ℎ^𝑡,Covℎ_2𝑡,ℎ^2𝑡−𝑘=(𝛼+𝛽)^𝑘𝐸ℎ^3𝑡−1−𝐸ℎ^2𝑡−1𝐸ℎ_𝑡−1,,ℎ^2𝑡−𝑘=𝐴^𝑘𝑉ℎ^2𝑡+𝐵(𝛼+𝛽)^𝑘−𝐴^𝑘(𝛼+𝛽)−𝐴𝐸ℎ^3𝑡−𝐸ℎ^2𝑡𝐸ℎ_𝑡,(2.9) where 𝐴=3𝛼²+𝛽²+2𝛼𝛽 and 𝐵=2[2𝛼²𝛾²+(𝜔+𝛼𝛾²)(𝛼+𝛽)].

Lemma 2.5. CovCℎov_𝑡ℎ_𝑡,𝜀_𝑡−𝑘=𝐸ℎ_𝑡𝜀_𝑡−𝑘=−2𝛼𝛾(𝛼+𝛽)^𝑘−1𝐸ℎ_𝑡,Covℎℎ_{𝑡−𝑘𝑡},𝜀_𝑡−𝑘=𝐸ℎ_𝑡ℎ_𝑡−𝑘𝜀_𝑡−𝑘=−2𝛼𝛾(𝛼+𝛽)^𝑘−1𝐸ℎ^2𝑡−1,Covℎ_2𝑡,𝜀^2𝑡−𝑘=(𝛼+𝛽)^𝑘𝑉ℎ_𝑡+2𝛼(𝛼+𝛽)^𝑘−1𝐸ℎ^2𝑡,,𝜀_𝑡−𝑘=−𝐸4𝛼ℎ𝛾^2𝑡𝜀𝜔_𝑡−𝑘+=𝛼−𝛾₂4𝛼𝛾(3𝛼+𝛽)𝐴^𝑘−1𝐸ℎ^2𝑡(𝛼+𝛽)^𝑘−𝐴^𝑘(𝛼+𝛽)−𝐴+2𝛼²𝛾²(𝛼+𝛽)^𝑘−1−𝐴^𝑘−1Covℎ_2𝑡(𝛼+𝛽)−𝐴𝐸ℎ^𝑡,,𝜀^2𝑡−𝑘=𝐴^𝑘𝐸ℎ^3𝑡−𝐸ℎ^2𝑡𝐸ℎ_𝑡+𝐵(𝛼+𝛽)^𝑘−𝐴^𝑘(𝛼+𝛽)−𝐴𝑉ℎ^𝑡+4𝛼(3𝛼+𝛽)𝐴^𝑘−1+2𝛼𝐵(𝛼+𝛽)_𝑘−1𝐸ℎ^3𝑡−𝐴^𝑘−1(𝛼+𝛽)−𝐴+42𝛼²𝛾²+𝜔+𝛼𝛾²𝐴^𝑘−1Covℎ_2𝑡𝐸ℎ^2𝑡,ℎ_𝑡−𝑘,𝜀_𝑡−𝑘=−𝐸2𝛼ℎ𝛾^2𝑡𝐵ℎ(𝛼_𝑡−𝑘𝜀+𝛽_𝑡−𝑘)=_𝑘−1−4𝛼𝛾𝐴^𝑘−1(3𝛼+𝛽)𝐸ℎ^3𝑡+𝜔+𝛼𝛾²𝐸ℎ^2𝑡−𝐴^𝑘−1(𝛼+𝛽)−𝐴𝐸ℎ^2𝑡,(2.10) where 𝐴 and 𝐵 are given in Lemma 2.4.

Furthermore, from Arvanitis and Demos [6] we know that the following results hold.

Theorem 2.6. The autocovariance of returns for the model in (2.1)–(2.3) is given by 𝛾_𝑘=Cov𝑟_𝑡,𝑟_𝑡−𝑘=𝛿²Covℎ_𝑡,ℎ_𝑡−𝑘+𝛿𝐸ℎ_𝑡𝜀_𝑡−𝑘+𝜑^𝑘𝜑^2𝑢1−𝜑²𝐸ℎ_𝑡ℎ_𝑡−𝑘,(2.11) and the covariance of squares’ levels and the autocovariance of squares are Cov𝑟^2𝑡,𝑟_𝑡−𝑘=+𝐸𝐸𝛿^2𝑡𝛿𝛿_{𝑡−𝑘𝑡−𝑘}Covℎ^2𝑡,ℎ_𝑡−𝑘+Cov𝛿^2𝑡,𝛿_𝑡−𝑘𝐸ℎ^2𝑡𝐸ℎ_𝑡−𝑘Cov𝑟_2𝑡Covℎ_𝑡,ℎ_𝑡−𝑘+𝐸𝛿^2𝑡𝐸ℎ^2𝑡𝜀_𝑡−𝑘+𝐸ℎ_𝑡𝜀_𝑡−𝑘,,𝑟^2𝑡−𝑘=+𝐸𝐸²𝛿𝛿^2𝑡_2𝑡Covℎ^2𝑡,ℎ^2𝑡−𝑘+Cov𝛿^2𝑡,𝛿^2𝑡−𝑘𝐸ℎ^2𝑡ℎ^2𝑡−𝑘+𝐸𝛿_2𝑡Covℎ^2𝑡,𝜀^2𝑡−𝑘+2𝐸𝛿^2𝑡𝛿_𝑡−𝑘𝐸ℎ^2𝑡ℎ_𝑡−𝑘𝜀_𝑡−𝑘Covℎ_𝑡,ℎ^2𝑡−𝑘+Covℎ_𝑡,𝜀^2𝑡−𝑘+2𝛿𝐸ℎ_𝑡ℎ_𝑡−𝑘𝜀_𝑡−𝑘,(2.12) where all needed covariances and expectations of the right-hand sizes are given in Lemmas 2.4 and 2.5, Cov𝛿^2𝑡,𝛿_𝑡−𝑘=Cov𝛿_𝑡,𝛿^2𝑡−𝑘=2𝜑^𝑘𝛿𝜑^2𝑢1−𝜑²Cov𝛿_2𝑡,,𝛿^2𝑡−𝑘=4𝜑^𝑘𝛿²𝜑^2𝑢1−𝜑²+2𝜑^2𝑘𝜑^4𝑢1−𝜑²².(2.13)

From the above theorems and lemmas it is obvious that our model can accommodate all stylized facts. For example, negative asymmetry is possible; volatility clustering and the leverage effect (negative Cov(ℎ_𝑡,𝜀_𝑡−𝑘)) can be accommodated, and so forth. Furthermore, the model can accommodate negative autocorrelations, 𝛾_𝑘, something that is not possible for the GARCH-M model (see Fiorentini and Sentana [12]). Finally, another interesting issue is the diffusion limit of the time-varying GQARCH-M process. As already presented by Arvanitis [13], the weak convergence of the Time-varying GQARCH(1,1)-M coincides with the general conclusions presented elsewhere in the literature. These are that weak limits of the endogenous volatility models are exogenous (stochastic volatility) continuous-time processes. Moreover, Arvanitis [13] suggests that there is a distributional relation between the GQARCH model and the continuous-time Ornstein-Uhlenbeck models with respect to appropriate nonnegative Levy’s processes.

Let us turn our attention to the estimation of our model. We will show that estimating our model is a hard task and the use of well-known methods such as the EM-algorithm cannot handle the problem due to the huge computational load that such methods require.

3. Likelihood-Inference: EM and Bayesian Approaches

The purpose of this section is the estimation of the time-varying GQARCH(1,1)-M model. Since our model involves two unobserved components (one from the time-varying in-mean parameter and one from the error term), the estimation method required is an EM and more specifically a simulated EM (SEM), as the expectation terms at the 𝐸 step cannot be computed. The main modern way of carrying out likelihood inference in such situations is via a Markov chain Monte Carlo (MCMC) algorithm (see Chib [14] for an extensive review). This simulation procedure can be used either to carry out Bayesian inference or to classically estimate the parameters by means of a simulated EM algorithm.

The idea behind the MCMC methods is that in order to sample a given probability distribution, which is referred to as the target distribution, a suitable Markov chain is constructed (using a Metropolis-Hasting (M-H) algorithm or a Gibbs sampling method) with the property that its limiting, invariant distribution is the target distribution. In most problems, the target distribution is absolutely continuous, and as a result the theory of MCMC methods is based on that of the Markov chains on continuous state spaces [15]. This means that by simulating the Markov chain a large number of times and recording its values a sample of (correlated) draws from the target distribution can be obtained. It should be noted that the Markov chain samplers are invariant by construction, and, therefore, the existence of the invariant distribution does not have to be checked in any particular application of MCMC method.

The Metropolis-Hasting algorithm (M-H) is a general MCMC method to produce sample variates from a given multivariate distribution. It is based on a candidate generating density that is used to supply a proposal value that is accepted with probability given as the ratio of the target density times the ratio of the proposal density. There are a number of choices of the proposal density (e.g., random walk M-H chain, independence M-H chain, tailored M-H chain) and the components may be revised either in one block or in several blocks. Another MCMC method, which is special case of the multiple block M-H method with acceptance rate always equal to one, is called the Gibbs sampling method and was brought into statistical prominence by Gelfand and Smith [16]. In this algorithm, the parameters are grouped into blocks, and each block is sampled according to the full conditional distribution denoted as 𝜋(𝜙_𝑡/𝜙_/𝑡). By Bayes' theorem, we have 𝜋(𝜙_𝑡/𝜙_/𝑡)∝𝜋(𝜙_𝑡𝜙_/𝑡), the joint distribution of all blocks, and so full conditional distributions are usually quite simply derived. One cycle of the Gibbs sampling algorithm is completed by simulating {𝜙_𝑡}^𝑝𝑡=1, where 𝑝 is the number of blocks, from the full conditional distributions, recursively updating the conditioning variables as one moves through each distribution. Under some general conditions, it is verified that the Markov chain generated by the M-H or the Gibbs sampling algorithm converges to the target density as the number of iterations becomes large.

Within the Bayesian framework, MCMC methods have proved very popular, and the posterior distribution of the parameters is the target density (see [17]). Another application of the MCMC is the analysis of hidden Markov’s models where the approach relies on augmenting the parameter space to include the unobserved states and simulate the target distribution via the conditional distributions (this procedure is called data augmentation and was pioneered by Tanner and Wong [18]). Kim et al. [19] discuss an MCMC algorithm of the stochastic volatility (SV) model which is an example of a state space model in which the state variable ℎ_𝑡 (log-volatility) appears non-linearly in the observation equation. The idea is to approximate the model by a conditionally Gaussian state space model with the introduction of multinomial random variables that follow a seven-point discrete distribution.

The analysis of a time-varying GQARCH-M model becomes substantially complicated since the log-likelihood of the observed variables can no longer be written in closed form. In this paper, we focus on both the Bayesian and the classical estimation of the model. Unfortunately, the non-Markovian nature of the GARCH process implies that each time we simulate one error we implicitly change all future conditional variances. As pointed out by Shephard [20], a regrettable consequence of this path dependence in volatility is that standard MCMC algorithms will evolve in 𝑂(𝑇²) computational load (see [21]). Since this cost has to be borne for each parameter value, such procedures are generally infeasible for large financial datasets that we see in practice.

3.1. Estimation Problem: Simulated EM Algorithm

As mentioned already, the estimation problem arises because of the fact that we have two unobserved processes. More specifically, we cannot write down the likelihood function in closed form since we do not observe both 𝜀_𝑡 and 𝛿_𝑡. On the other hand, the conditional log-likelihood function of our model assuming that 𝛿_𝑡 were observed would be the following: ℓ𝐫,𝜹∣𝜙,ℱ₀==l−n𝑇𝑝ln𝐫2∣𝜋𝜹−,𝜙1,ℱ₀+ln𝑝𝜹∣𝜙,ℱ₀2^𝑇_𝑡=1lnℎ_𝑡−12^𝑇_𝑡=1𝜀_𝑡²ℎ_𝑡−𝑇2ln𝜑^2𝑢−12^𝑇_𝑡=1𝛿_𝑡−𝛿(1−𝜑)−𝜑𝛿_𝑡−1²𝜑^2𝑢,(3.1) where 𝐫=(𝑟₁,…,𝑟_𝑇)^,𝜹=(𝛿₁,…,𝛿_𝑇)^, and 𝐡=(ℎ₁,…,ℎ_𝑇)^.

However, the 𝛿_𝑡’s are unobserved, and, thus, to classically estimate the model, we have to rely on an EM algorithm [22] to obtain estimates as close to the optimum as desired. At each iteration, the EM algorithm obtains 𝜙^(𝑛+1), where 𝜙 is the parameter vector, by maximizing the expectation of the log-likelihood conditional on the data and the current parameter values, that is, 𝐸(ℓ(⋅)∣𝐫,𝜙^(𝑛),ℱ₀) with respect to 𝜙 keeping 𝜙^(𝑛) fixed.

The 𝐸 step, thus, requires the expectation of the complete log-likelihood. For our model, this is given by: 𝐸ℓ(⋅)∣𝐫,𝜙^(𝑛),ℱ₀=−𝑇ln2𝜋−𝑇2ln𝜑^2𝑢−12^𝑇_𝑡=1−1𝐸lnℎ_𝑡∣𝐫,𝜙^(𝑛),ℱ₀2^𝑇_𝑡=1𝐸𝜀^𝑡²ℎ_𝑡∣𝐫,𝜙^(𝑛),ℱ₀−12^𝑇_𝑡=1𝐸𝛿^𝑡−𝛿(1−𝜑)−𝜑𝛿_𝑡−1²𝜑^2𝑢∣𝐫,𝜙^(𝑛),ℱ₀.(3.2)

It is obvious that we cannot compute such quantities. For that matter, we may rely on a simulated EM where the expectation terms are replaced by averages over simulations, and so we will have an SEM or a simulated score. The SEM log-likelihood is: SEMℓ=−𝑇ln2𝜋−𝑇2ln𝜑^2𝑢−12𝑀1^𝑀_𝑇𝑖=1_𝑡=1lnℎ^𝑡(𝑖)−12𝑀1^𝑀_𝑇𝑖=1_𝑡=1𝜀^𝑡(𝑖)²ℎ^𝑡(𝑖)−𝑇2(1−𝜑)²𝛿²𝜑^2𝑢−12𝜑1^2𝑢1𝑀^𝑀_𝑇𝑖=1_𝑡=1𝛿^𝑡(𝑖)²+(1−𝜑)𝛿𝜑^2𝑢1𝑀^𝑀_𝑇𝑖=1_𝑡=1𝛿^𝑡(𝑖)+𝜑𝜑^2𝑢1𝑀^𝑀_𝑇𝑖=1_𝑡=1𝛿^𝑡(𝑖)𝛿^{(𝑖)𝑡−1}+(1−𝜑)𝜑𝛿𝜑^2𝑢1𝑀^𝑀_𝑇𝑖=1_𝑡=1𝛿^{(𝑖)𝑡−1}−𝜑²2𝜑^2𝑢1𝑀^𝑀_𝑇𝑖=1_𝑡=1𝛿^{(𝑖)𝑡−1}².(3.3) Consequently, we need to obtain the following quantities: (1/𝑀)∑^𝑀𝑖=1∑^𝑇𝑡=1lnℎ^𝑡(𝑖), (1/𝑀)∑^𝑀𝑖=1∑^𝑇𝑡=1((𝜀^𝑡(𝑖))²/ℎ^𝑡(𝑖)), (1/𝑀)∑^𝑀𝑖=1∑^𝑇𝑡=1𝛿^𝑡(𝑖), (1/𝑀)∑^𝑀𝑖=1∑^𝑇𝑡=1𝛿^{(𝑖)𝑡−1},  (1/𝑀)∑^𝑀𝑖=1∑^𝑇𝑡=1𝛿^𝑡(𝑖)𝛿^{(𝑖)𝑡−1} and (1/𝑀)∑^𝑀𝑖=1∑^𝑇𝑡=1(𝛿^𝑡(𝑖))², (1/𝑀)∑^𝑀𝑖=1∑^𝑇𝑡=1𝛿^{2(𝑖)𝑡−1}, where 𝑀 is the number of simulations.

Thus, to classically estimate our model by using an SEM algorithm, the basic problem is to sample from 𝐡∣𝜙,𝐫,ℱ₀ where 𝜙 is the vector of the unknown parameters and also sample from 𝜹∣𝜙,𝐫,ℱ₀.

In terms of identification, the model is not, up to second moment, identified (see Corollary 1 in Sentana and Fiorentini [23]). The reason is that we can transfer unconditional variance from the error, 𝜀_𝑡, to the price of risk, 𝛿_𝑡, and vice versa. One possible solution is to fix 𝜔 such that 𝐸(ℎ_𝑡) is 1 or to set 𝜑_𝑢 to a specific value. In fact in an earlier version of the paper, we fixed 𝜑_𝑢 to be 1 (see Anyfantaki and Demos [24]). Nevertheless, from a Bayesian viewpoint, the lack of identification is not too much of a problem, as the parameters are identified through their proper priors (see Poirier [25]).

Next, we will exploit the Bayesian estimation of the model, and, since we need to resort to simulations, we will show that the key task is again to simulate from 𝜹∣𝜙,𝐫,ℱ₀.

3.2. Simulation-Based Bayesian Inference

In our problem, the key issue is that the likelihood function of the sample 𝑝(𝐫∣𝜙,ℱ₀) is intractable which precludes the direct analysis of the posterior density 𝑝(𝜙∣𝐫,ℱ₀). This problem may be overcome by focusing instead on the posterior density of the model using Bayes’ rule:𝑝(𝜙,𝜹∣𝐫)∝𝑝(𝜙,𝜹)𝑝(𝐫∣𝜙,𝜹)∝𝑝(𝜹∣𝜙)𝑝(𝜙)𝑝(𝐫∣𝜙,𝜹),(3.4) where𝜙=𝛿,𝜑,𝜑_𝑢,𝛼,𝛽,𝛾,𝜔^.(3.5)

Now,𝑝(𝜹∣𝜙)=^𝑇_𝑡=1𝑝𝛿^𝑡𝛿,𝜑,𝜑^2𝑢=^𝑇_𝑡=11√2𝜋𝜑^2𝑢exp−𝛿^𝑡−𝛿(1−𝜑)−𝜑𝛿_𝑡−1²2𝜑^2𝑢.(3.6)

On the other hand,𝑝(𝐫∣𝜙,𝜹)=^𝑇_𝑡=1𝑝𝑟^𝑡𝑟_𝑡−1,𝜹,𝜙=^𝑇_𝑡=11√2𝜋ℎ_𝑡exp−𝜀^2𝑡2ℎ_𝑡(3.7) is the full-information likelihood. Once we have the posterior density, we get the parameters' marginal posterior density by integrating the posterior density. MCMC is one way of numerical integration.

The Hammersley-Clifford theorem (see Clifford [26]) says that a joint distribution can be characterized by its complete conditional distribution. Hence, given initial values {𝛿_𝑡}⁽⁰⁾,𝜙⁽⁰⁾, we draw {𝛿_𝑡}⁽¹⁾ from 𝑝({𝛿_𝑡}⁽¹⁾∣𝐫,𝜙⁽⁰⁾) and then 𝜙⁽¹⁾ from 𝑝(𝜙⁽¹⁾∣{𝛿_𝑡}⁽¹⁾,𝐫). Iterating these steps, we finally get ({𝛿_𝑡}^(𝑖),𝜙^(𝑖))^𝑀𝑖=1, and under mild conditions it is shown that the distribution of the sequence converges to the joint posterior distribution 𝑝(𝜙,𝜹∣𝐫).

The above simulation procedure may be carried out by first dividing the parameters into two blocks:𝜙𝜙₁₂=𝛿,𝜑,𝜑^2𝑢,=(𝛼,𝛽,𝛾,𝜔),(3.8) Then the algorithm is described as follows.(1)Initialize 𝜙.(2)Draw from 𝑝(𝛿_𝑡∣𝛿_≠𝑡,𝐫,𝜙).(3)Draw from 𝑝(𝜙∣𝜹,𝐫) in the following blocks:(i)draw from 𝑝(𝜙₁∣𝜹,𝐫) using the Gibbs sampling. This is updated in one block;(ii)draw from 𝑝(𝜙₂∣𝐫) by M-H. This is updated in a second block.(4)Go to  (2).

We review the implementation of each step.

3.2.1. Gibbs Sampling

The task of simulating from an AR model has been already discussed. Here, we will follow the approach of Chib [27], but we do not have any MA terms which makes inference simpler.

Suppose that the prior distribution of (𝛿,𝜑^2𝑢,𝜑) is given by:𝑝𝛿,𝜑^2𝑢,𝜑=𝑝𝛿∣𝜑^2𝑢𝑝𝜑^2𝑢𝑝(𝜑),(3.9) which means that 𝛿,𝜑^2𝑢 is a priory independent of 𝜑.

Also the following holds for the prior distributions of the parameter subvector 𝜙₁: 𝑝𝛿𝑝∣𝜑𝜑^2𝑢_2𝑢∼𝑁𝛿_𝑝𝑟,𝜑^2𝑢𝜎^2𝛿_𝑝𝑟,𝑝(𝜑)∼∼𝑁IG𝜑𝑣₀2⁰,𝑑2⁰,,𝜎^2𝜑₀𝐼_𝜑,(3.10) where 𝐼_𝜑 ensures that 𝜑 lies outside the unit circle, IG is the inverted gamma distribution, and the hyperparameters 𝑣₀,𝑑₀,𝛿_𝑝𝑟,𝜎^2𝛿_𝑝𝑟,𝜑₀,𝜎^2𝜑₀ have to be defined.

Now, the joint posterior is proportional to 𝑝𝛿,𝜑,𝜑^2𝑢∣𝐫,𝜹∝^𝑇_𝑡=11√2𝜋𝜑^2𝑢exp−𝛿^𝑡−(1−𝜑)𝛿−𝜑𝛿_𝑡−1²×𝑁𝛿_𝑝𝑟2𝜑^2𝑢,𝜑^2𝑢𝜎^2𝛿_𝑝𝑟×IG𝑣⁰2,𝑑2⁰×𝑁𝜑⁰,𝜎^2𝜑₀𝐼_𝜑.(3.11) From a Bayesian viewpoint, the right-hand side of the above equation is equal to the “augmented” prior, that is, the prior augmented by the latent 𝜹 (We would like to thank the associate editor for bringing this to our attention.) We proceed to the generation of these parameters.

Generation of 𝛿
First we see how to generate 𝛿. Following again Chib [27], we may write 𝛿^∗𝑡=𝛿_𝑡−𝜑𝛿_𝑡−1,𝛿^∗𝑡∣ℱ_𝑡−1∼𝑁(1−𝜑)𝛿,𝜑^2𝑢,(3.12) or, otherwise, 𝛿^∗𝑡=(1−𝜑)𝛿+𝑣_𝑡,𝑣_𝑡∼𝑁0,𝜑^2𝑢.(3.13) Under the above and using Chib's (1993) notation, we have that the proposal distribution is the following Gaussian distribution (see Chib [27] for a proof).Proposition 3.1. The proposal distribution of 𝛿 is 𝛿∣𝜹,𝜙,𝜑^2𝑢∼𝑁𝛿̃,𝜑^2𝑢𝜎^2𝛿,(3.14) where 𝛿̃=𝜎^2𝛿𝜎𝛿^{𝑝𝑟2𝛿}_𝑝𝑟+(1−𝜑)^𝑇_𝑡=1𝛿^∗𝑡𝜎_2𝛿,=1𝜎^2𝛿_𝑝𝑟+(1−𝜑)²⁻¹.(3.15)
Hence, the generation of 𝛿 is completed, and we may turn on the generation of the other parameters.

Generation of 𝜑^2𝑢
For the generation of 𝜑^2𝑢 and using [27] notation, we have the following.
Proposition 3.2. The proposal distribution of 𝜑^2𝑢 is 𝜑^2𝑢∣𝜹,𝜙,𝛿∼IG𝑇−𝑣⁰2,𝑑⁰+𝑄+𝑑2,(3.16) where 𝑑=𝑄_𝑇=𝛿−𝛿_p𝑟²𝜎^𝛿−2_𝑝𝑟,_𝑡=2𝛿^∗𝑡−𝛿(1−𝜑)².(3.17)
Finally, we turn on the generation of 𝜑.

Generation of 𝜑
For the generation of 𝜑, we follow again Chib [27] and write 𝛿_𝑡=(1−𝜑)𝛿−𝜑𝛿_𝑡−1+𝑣_𝑡,𝑣_𝑡∼𝑁0,𝜑^2𝑢.(3.18) We may now state the following proposition (see Chib [27] for a proof).
Proposition 3.3. The proposal distribution of 𝜑 is 𝜑²∣𝜹,𝛿,𝜑^2𝑢∼𝑁𝜑,𝜎^2𝜑,(3.19) where 𝜑=𝜎^2𝜑𝜎^𝜑−20𝜑₀+𝜑^𝑢𝑇−2_𝑡=1𝛿_𝑡−1𝜎_𝜑−2−𝛿𝛿_𝑡−𝛿,=𝜎^𝜑−2₀+𝜑^𝑢𝑇−2_𝑡=1𝛿_𝑡−1−𝛿².(3.20)
The Gibbs sampling scheme has been completed, and the next step of the algorithm requires the generation of the conditional variance parameters via an M-H algorithm which is now presented.

3.2.2. Metropolis-Hasting

Step (3)-(ii) is the task of simulating from the posterior of the parameters of a GQARCH-M process. This has been already addressed by Kim et al. [19], Bauwens and Lubrano [28], Nakatsuma [29], Ardia [30] and others.

First, we need to decide on the priors. For the parameters 𝛼,𝛽,𝛾,𝜔, we use normal densities as priors:𝑝𝑝((𝜶𝛽))∼∼𝑁𝑁𝜇𝜇_𝛼𝛽,Σ_𝜶𝐼_𝜶,,𝜎^2𝛽𝐼^𝛽(3.21) and similarly 𝑝(𝛾)∼𝑁𝜇_𝛾,𝜎^2𝛾,(3.22) where 𝜶=(𝜔,𝛼)^,𝐼_𝜶,𝐼_𝛽 are the indicators ensuring the constraints 𝜶>0,𝛼+𝛽<1 and 𝛽>0,𝛼+𝛽<1, respectively. 𝜇_.,𝜎^2. are the hyperparameters.

We form the joint prior by assuming prior independence between 𝜶,𝛽,𝛾, and the joint posterior is then obtained by combining the joint prior and likelihood function by Bayes' rule:𝑝(𝜶,𝛽,𝛾∣𝐫)∝^𝑇_𝑡=11√2𝜋ℎ_𝑡exp−𝜀^2𝑡2ℎ_𝑡×𝑁𝜇^𝛼,Σ_𝜶𝐼_𝜶×𝑁𝜇_𝛽,𝜎^2𝛽𝐼^𝛽×𝑁𝜇_𝛾,𝜎^2𝛾.(3.23)

For the M-H algorithm, we use the following approximated GARCH model as in Nakatsuma [29] which is derived by the well-known property of GARCH models [4]: 𝜀^2𝑡=𝜔+𝛼𝜀_𝑡−1−𝛾²+𝛽𝜀^2𝑡−1+𝑤_𝑡−𝛽𝑤_𝑡−1,(3.24) where 𝑤_𝑡=𝜀^2𝑡−ℎ_𝑡 with 𝑤_𝑡∼𝑁(0,2ℎ^2𝑡).

Then the corresponding approximated likelihood is written as 𝑝𝜀²∣𝜹,𝐫,𝜙₂=^𝑇_𝑡=112ℎ_𝑡√𝜋exp⎡⎢⎢⎢⎣−𝜀^2𝑡−𝜔−𝛼𝜀_𝑡−1−𝛾²−𝛽𝜀^2𝑡−1+𝛽𝑤_𝑡−1²4ℎ^2𝑡⎤⎥⎥⎥⎦,(3.25) and the generation of 𝜶,𝛽,𝛾 is based on the above likelihood where we update {ℎ_𝑡} each time after the corresponding parameters are updated. The generation of the four variance parameters is given.

Generation of 𝛼
For the generation of 𝜶, we first note that 𝑤_𝑡 in (3.32), below, can be written as a linear function of 𝜶: 𝑤_𝑡=𝜀^2𝑡−𝜁_𝑡𝜶,(3.26) where 𝜁_𝑡=[𝜄_𝑡̃,𝜀̂^2𝑡] with 𝜀𝜀̃^2𝑡_2𝑡=𝜀̃^2𝑡−𝛽𝜀̃^2𝑡−1,𝜀̂_2𝑡=𝜀^2𝑡+𝛽𝜀̃^2𝑡−1,𝜄_𝑡̃=𝜀_𝑡−1−𝛾²+𝛽𝜀̂^2𝑡−1,=1+𝛽𝜄̃_𝑡−1.(3.27)
Now, let the two following vectors be𝑋𝑌_𝜶𝜶=𝜀²¹,…,𝜀^2𝑇^,=𝜁^1,…,𝜁^𝑇^.(3.28)
Then the likelihood function of the approximated model is rewritten as 𝑝𝜀²∣𝐫,𝜹,𝜙₂=^𝑇_𝑡=112𝜋2ℎ^2𝑡exp⎡⎢⎢⎢⎣−𝜀^2𝑡−𝜁_𝑡𝜶²22ℎ^2𝑡⎤⎥⎥⎥⎦.(3.29)
Using this we have the following proposal distribution of 𝜶 (see Nakatsuma [29] or Ardia [30] for a proof).
Proposition 3.4. The proposal distribution of 𝜶 is 𝜶∣𝑌,𝑋,Σ,𝜙_2−𝜶∼𝑁𝜇_𝜶,Σ_𝜶𝐼^𝜶,(3.30) where 𝜇_𝜶=Σ_𝜶(𝑋^𝜶Λ⁻¹𝑌_𝜶+Σ^𝜶−1𝜇_𝜶), Σ_𝜶=(𝑋^𝜶Λ⁻¹𝑋_𝜶+Σ^𝜶−1)⁻¹, and Λ=diag(2ℎ²¹,…,2ℎ^2𝑇). 𝐼_𝜶 imposes the restriction that 𝜶>0 and 𝛼+𝛽<1.
Hence a candidate 𝜶 is sampled from this proposal density and accepted with probability: min𝑝𝜶,𝛽,𝛾∣𝜹,𝐫𝑞𝜶^∗∣𝜶,𝛽,𝛾,𝜹,𝐫𝑝(𝜶^∗,𝛽,𝛾∣𝜹,𝐫)𝑞𝜶∣𝜶^∗,𝛽,𝛾,𝜹,𝐫,1,(3.31) where 𝜶^∗ is the previous draw.
Similar procedure is used for the generation of 𝛽 and 𝛾.

Generation of 𝛽
Following Nakatsuma [29], we linearize 𝑤_𝑡 by the first-order Taylor expansion 𝑤_𝑡(𝛽)≈𝑤_𝑡𝛽^∗+𝜉_𝑡𝛽^∗𝛽−𝛽^∗,(3.32) where 𝜉_𝑡 is the first-order derivative of 𝑤_𝑡(𝛽) evaluated at 𝛽^∗ the previous draw of the M-H sampler.
Define as 𝑟_𝑡=𝑤_𝑡𝛽^∗+𝑔_𝑡𝛽^∗𝛽^∗,(3.33) where 𝑔_𝑡(𝛽^∗)=−𝜉_𝑡(𝛽^∗) which is computed by the recursion: 𝑔_𝑡=𝜀^2𝑡−1−𝑤_𝑡−1+𝛽^∗𝑔_𝑡−1,(3.34)𝜉_𝑡=0 for 𝑡≤0 [30]. Then, 𝑤_𝑡(𝛽)≈𝑟_𝑡−𝑔_𝑡(𝛽)𝛽.(3.35)
Let the following two vectors be 𝑋𝑌_𝛽𝛽=𝑟₁,…,𝑟_𝑇^,=𝑔₁,…,𝑔_𝑇^.(3.36)
Then, the likelihood function of the approximated model is rewritten as 𝑝𝜀²∣𝜹,𝐫,𝜙₂=^𝑇_𝑡=112𝜋2ℎ^2𝑡exp−𝑤^𝑡𝛽^∗+𝜉_𝑡𝛽^∗𝛽−𝛽^∗²22ℎ^2𝑡.(3.37)
We have the following proposal distribution for 𝛽 (for a proof see Nakatsuma [29] or Ardia [30]).

Proposition 3.5. The proposal distribution for 𝛽 is 𝛽∣𝑌,𝑋,𝜎^2𝛽,𝜙_2−𝛽∼𝑁𝜇_𝛽,𝜎^2𝛽𝐼^𝛽,(3.38) where 𝜇_𝛽=𝜎^2𝛽(𝑋^𝛽Λ⁻¹𝑌_𝛽+(𝜇_𝛽/𝜎^2𝛽)), 𝜎^2𝛽=(𝑋^𝛽Λ⁻¹𝑋_𝛽+(1/𝜎^2𝛽))⁻¹, and Λ=diag(2ℎ⁴¹,…,2ℎ^4𝑇). 𝐼_𝛽 imposes the restriction that 𝛽>0 and 𝛼+𝛽<1.
Hence, a candidate 𝛽̃ is sampled from this proposal density and accepted with probability: min𝑝𝛽̃,𝛼,𝛾∣𝜹,𝐫𝑞𝛽^∗∣𝛽̃,𝛼,𝛾,𝜹,𝐫𝑝𝛽^∗,𝛼,𝛾∣𝜹,𝐫𝑞𝛽̃∣𝛽^∗,𝛼,𝛾,𝜹,𝐫,1(3.39)

Finally, we explain the generation of 𝛾.

Generation of 𝛾
As with 𝛽, we linearize 𝑤_𝑡 by a first-order Taylor, expansion at a point 𝛾^∗ the previous draw in the M-H sampler. In this case, 𝑟_𝑡=𝑤_𝑡𝛾^∗−𝑔_𝑡𝛾^∗𝛾^∗,(3.40) where 𝑔_𝑡(𝛾^∗)=−𝜉_𝑡(𝛾^∗) which is computed by the recursion: 𝑔_𝑡=−2𝛼𝜀_𝑡−1−𝛾^∗+𝛽𝑔_𝑡−1,(3.41) and 𝑔_𝑡=0 for 𝑡≤0.
Then 𝑤_𝑡(𝛾)≈𝑟_𝑡−𝑔_𝑡𝛾.(3.42)
Let again 𝑋𝑌_𝛾𝛾=𝑟₁,…,𝑟_𝑇^,=𝑔₁,…,𝑔_𝑇^(3.43) and the likelihood function of the approximated model is rewritten as 𝑝𝜀²∣𝜹,𝐫,𝜙₂=^𝑇_𝑡=112𝜋2ℎ^2𝑡exp−𝑤^𝑡𝛾^∗−𝑔_𝑡𝛾^∗𝛾^∗²22ℎ^2𝑡.(3.44)
Thus, we have the following proposal distribution for 𝛾 (for proof see Nakatsuma [29] and Ardia [30]).

Proposition 3.6. The proposal distribution of 𝛾 is 𝛾∣𝑌,𝑋,𝜎^2𝛾,𝜙_2−𝛾∼𝑁𝜇_𝛾,𝜎^2𝛾,(3.45) where 𝜇_𝛾=𝜎^2𝛾(𝑋^𝛾Λ⁻¹𝑌_𝛾+(𝜇_𝛾/𝜎^2𝛾)), 𝜎^2𝛾=(𝑋^𝛾Λ⁻¹𝑋_𝛾+(1/𝜎^2𝛾))⁻¹, and Λ=diag(2ℎ⁴¹,…,2ℎ^4𝑇). A candidate 𝛾̃ is sampled from this proposal density and accepted with probability: min𝑝(𝛾̃,𝛼,𝛽∣𝜹,𝐫)𝑞𝛾^∗∣𝛾̃,𝛼,𝛽,𝜹,𝐫𝑝(𝛾^∗,𝛼,𝛽∣𝜹,𝐫)𝑞(𝛾̃∣𝛾^∗,𝛼,𝛽,𝜹,𝐫),1.(3.46)

The algorithm described above is a special case of a MCMC algorithm, which converges as it iterates, to draws from the required density 𝑝(𝜙,𝜹∣𝐫). Posterior moments and marginal densities can be estimated (simulation consistently) by averaging the relevant function of interest over the sample variates. The posterior mean of 𝜙 is simply estimated by the sample mean of the simulated 𝜙 values. These estimated values can be made arbitrarily accurate by increasing the simulation sample size. However, it should be remembered that sample variates from an MCMC algorithm are a high dimensional (correlated) sample from the target density, and sometimes the serial correlation can be quite high for badly behaved algorithms.

All that remains, therefore, is step (2). Thus, from the above, it is seen that the main task is again as with the classical estimation of the model, to simulate from 𝜹∣𝜙,𝐫,ℱ₀.

3.2.3. MCMC Simulation of 𝜀∣𝜙,𝐫,ℱ₀

For a given set of parameter values and initial conditions, it is generally simpler to simulate {𝜀_𝑡} for 𝑡=1,…,𝑇 and then compute {𝛿_𝑡}^𝑇𝑡=1 than to simulate {𝛿_𝑡}^𝑇𝑡=1 directly. For that matter, we concentrate on simulators of 𝜀_𝑡 given 𝐫 and 𝜙. We set the mean and the variance of 𝜀₀ equal to their unconditional values, and, given that ℎ_𝑡 is a sufficient statistic for ℱ_𝑡−1 and the unconditional variance is a deterministic function of 𝜙, ℱ₀ can be eliminated from the information set without any information loss.

Now sampling from 𝑝(𝜀∣𝐫,𝜙)∝𝑝(𝐫∣𝜀,𝜙)𝑝(𝜀∣𝜙) is feasible by using an M-H algorithm where we update each time only one 𝜀_𝑡 leaving all the other unchanged [20]. In particular, let us write the nth iteration of a Markov chain as 𝜀^𝑛𝑡. Then we generate a potential new value of the Markov chain 𝜀^new_𝑡 by proposing from some candidate density 𝑔(𝜀_𝑡∣𝜀^𝑛⧵𝑡,𝐫,𝜙) where 𝜀^𝑛⧵𝑡={𝜀^1𝑛+1,…,𝜀^{𝑛+1𝑡−1},𝜀^𝑛𝑡+1,…,𝜀^𝑛𝑇} which we accept with probabilitymin⎡⎢⎢⎣1,𝑝𝜀^new𝑡∣𝜀^𝑛⧵𝑡,𝐫,𝜙𝑔𝜀^new_𝑡∣𝜀^𝑛⧵𝑡,𝐫,𝜙𝑝𝜀^𝑛𝑡∣𝜀^𝑛⧵𝑡,𝐫,𝜙𝑔𝜀^𝑛𝑡∣𝜀^𝑛⧵𝑡,𝐫,𝜙⎤⎥⎥⎦.(3.47)

If it is accepted then, we set 𝜀^𝑡𝑛+1=𝜀^new_𝑡 and otherwise we keep 𝜀^𝑡𝑛+1=𝜀^𝑛𝑡. Although the proposal is much better since it is only in a single dimension, each time we consider modifying a single error we have to compute: 𝑝𝜀^new_𝑡∣𝜀^𝑛⧵𝑡,𝐫,𝜙𝑝𝜀^𝑛𝑡∣𝜀^𝑛⧵𝑡,𝐫,𝜙=𝑝𝑟^𝑡∣𝜀^new_𝑡,ℎ^new𝑡,𝑡,𝜙𝑝𝜀^new_𝑡∣ℎ^new𝑡,𝑡,𝜙𝑝𝑟_𝑡∣ℎ^{𝑡𝑛,𝑡},𝜙∗𝑝𝑟_𝑡∣ℎ^new𝑡,𝑡,𝜙𝑝𝑟_𝑡∣𝜀^𝑛𝑡,ℎ^{𝑡𝑛,𝑡},𝜙𝑝𝜀^new_𝑡∣ℎ^{𝑡𝑛,𝑡},𝜙_𝑇_𝑠_=𝑡+1𝑝𝑟_𝑠∣𝜀^𝑟𝑠,ℎ^new𝑠,𝑡,𝜙𝑝𝜀^𝑟𝑠∣ℎ^new𝑠,𝑡,𝜙𝑝𝑟_𝑠∣ℎ^{𝑠𝑛,𝑡},𝜙=𝑝𝑟𝑝𝑟_𝑠∣ℎ^new𝑠,𝑡,𝜙𝑝𝑟_𝑠∣𝜀^𝑛𝑠,ℎ^{𝑠𝑛,𝑡},𝜙𝑝𝜀^𝑛𝑠∣ℎ^{𝑠𝑛,𝑡},𝜙_𝑡∣𝜀^new_𝑡,ℎ^new𝑡,𝑡,𝜙𝑝𝜀^new_𝑡∣ℎ^new𝑡,𝑡,𝜙∗𝑝𝑟_𝑡∣𝜀^𝑛𝑡,ℎ^{𝑡𝑛,𝑡},𝜙𝑝𝜀^new_𝑡∣ℎ^{𝑡𝑛,𝑡},𝜙_𝑇_𝑠_=𝑡+1𝑝𝑟_𝑠∣𝜀^𝑛𝑠,ℎ^new𝑠,𝑡,𝜙𝑝𝜀^𝑛𝑠∣ℎ^new𝑠,𝑡,𝜙𝑝𝑟_𝑠∣𝜀^𝑛𝑠,ℎ^{𝑠𝑛,𝑡},𝜙𝑝𝜀^𝑛𝑠∣ℎ^{𝑠𝑛,𝑡},𝜙,(3.48) where for 𝑠=𝑡+1,…,𝑇,ℎℎ^new𝑠,𝑡_{𝑠𝑛,𝑡}=𝑉𝜀_𝑠∣𝜀^𝑛𝑠−1,𝜀^𝑛𝑠−2,…,𝜀^𝑛𝑡+1,𝜀^new_𝑡,𝜀^{𝑛+1𝑡−1},…,𝜀^1𝑛+1,=𝑉𝜀_𝑠∣𝜀^𝑛𝑠−1,𝜀^𝑛𝑠−2,…,𝜀^𝑛𝑡+1,𝜀^𝑛𝑡,𝜀^{𝑛+1𝑡−1},…,𝜀^1𝑛+1(3.49) whileℎ^new𝑡,𝑡=ℎ^{𝑡𝑛,𝑡}.(3.50)

Nevertheless, each time we revise one 𝜀_𝑡, we have also to revise 𝑇−𝑡 conditional variances because of the recursive nature of the GARCH model which makes ℎ^new𝑠,𝑡 depend upon 𝜀^new_𝑡 for 𝑠=𝑡+1,…,𝑇. And since 𝑡=1,…,𝑇, it is obvious that we need to calculate 𝑇² normal densities, and so this algorithm is 𝑂(𝑇²). And this should be done for every 𝜙. To avoid this huge computational load, we show how to use the method proposed by Fiorentini et al. [10] and so do MCMC with only 𝑂(𝑇) calculations. The method is described in the following subsection.

3.3. Estimation Method Proposed: Classical and Bayesian Estimation

The method proposed by Fiorentini et al. [10] is to transform the GARCH model into a first-order Markov’s model and so do MCMC with only 𝑂(𝑇) calculations. Following their transformation, we augment the state vector with the variables ℎ_𝑡+1 and then sample the joint Markov process {ℎ_𝑡+1,𝑠_𝑡}∣𝐫,𝜙∈ℱ_𝑡 where 𝑠_𝑡=sign𝜀_𝑡−𝛾,(3.51) so that 𝑠_𝑡=±1 with probability one. The mapping is one to one and has no singularities. More specifically, if we know {ℎ_𝑡+1} and 𝜑, then we know the value of𝜀_𝑡−𝛾²=ℎ^𝑡+1−𝜔−𝛽ℎ_𝑡𝛼∀𝑡≥1.(3.52)

Hence the additional knowledge of the signs of (𝜀_𝑡−𝛾) would reveal the entire path of {𝜀_𝑡} so long as ℎ₀ (which equals the unconditional value in our case) is known, and, thus, we may now reveal also the unobserved random variable {𝛿_𝑡}∣𝐫,𝜙,{ℎ_𝑡+1}.

Now we have to sample from 𝑝𝑠_𝑡,ℎ_𝑡+1∣𝐫,𝜙∝^𝑇_𝑡=1𝑝𝑠_𝑡∣ℎ_𝑡+1,ℎ_𝑡,𝜙𝑝ℎ_𝑡+1∣ℎ_𝑡,𝜙𝑝𝑟_𝑡∣𝑠_𝑡,ℎ_𝑡,ℎ_𝑡+1,𝜙,(3.53) where the second and the third term come from the model, and the first comes from the fact that 𝜀_𝑡∣ℱ_𝑡−1∼𝑁(0,ℎ_𝑡) but 𝜀_𝑡∣{ℎ_𝑡+1},ℱ_𝑡−1 takes values𝜀_𝑡=𝛾±𝑑_𝑡,(3.54) where 𝑑_𝑡=ℎ_𝑡+1−𝜔−𝛽ℎ_𝑡𝛼.(3.55)

From the above, it is seen that we should first simulate {ℎ_𝑡+1}∣𝐫,𝜙 since we do not alter the volatility process when we flip from 𝑠_𝑡=−1 to 𝑠_𝑡=1 (implying that the signs do not cause the volatility process), but we do alter 𝜀_𝑡 and then simulate {𝑠_𝑡}∣{ℎ_𝑡+1},𝐫,𝜙. The second step is a Gibbs sampling scheme whose acceptance rate is always one and also conditional on {ℎ_𝑡+1},𝐫,𝜙 the elements of {𝑠_𝑡} are independent which further simplifies the calculations. We prefer to review first the Gibbs sampling scheme and then the simulation of the conditional variance.

3.3.1. Simulations of {𝑠_𝑡}∣{ℎ_𝑡+1},𝐫, 𝜙

First, we see how to sample from {𝑠_𝑡}∣{ℎ_𝑡+1},𝐫,𝜙. To obtain the required conditionally Bernoulli distribution, we establish first some notation. We have the following (see Appendix A): 𝑐_𝑡=1√𝑣_{𝑡∣𝑟𝑡,ℎ𝑡}𝜑𝛾+𝑑^𝑡−𝜀_{𝑡∣𝑟𝑡,ℎ𝑡}√𝑣_{𝑡∣𝑟𝑡,ℎ𝑡}+𝜑𝛾−𝑑^𝑡−𝜀_{𝑡∣𝑟𝑡,ℎ𝑡}𝜀√𝑣_{𝑡∣𝑟𝑡,ℎ𝑡}_{𝑡∣𝑟𝑡,ℎ𝑡},=𝐸𝜀_𝑡∣𝑟_𝑡,ℎ_𝑡=1−𝜑²𝑟_𝑡−𝛿ℎ_𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑²,𝑣_{𝑡∣𝑟𝑡,ℎ𝑡}=Var𝜀_𝑡∣𝑟_𝑡,ℎ_𝑡=𝜑^2𝑢ℎ^2𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑².(3.56)

Using the above notation, we see that the probability of drawing 𝑠_𝑡=1 conditional on {ℎ_𝑡+1} is equal to the probability of drawing 𝜀_𝑡=𝛾+𝑑_𝑡 conditional on ℎ_𝑡+1,ℎ_𝑡,𝑟_𝑡,𝜙, where 𝑑_𝑡 is given by (3.70), which is given by𝑝𝑠_𝑡=1∣ℎ_𝑡+1,𝐫,𝜙==𝑝𝜀_𝑡1=𝛾+𝑑_𝑡∣ℎ_𝑡+1,ℎ_𝑡,𝑟_𝑡,𝜙𝑐_𝑡√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}𝜑𝛾+𝑑^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}.(3.57)

Similarly for the probability of drawing 𝑠_𝑡=−1. Both of these quantities are easy to compute; for example,𝜑𝛾+𝑑^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}=√12𝜋exp⎧⎪⎨⎪⎩−21𝛾+𝑑^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}²⎫⎪⎬⎪⎭(3.58) and so we may simulate {𝑠_𝑡}∣{ℎ_𝑡+1},𝐫,𝜙 using a Gibbs sampling scheme. Specifically, since conditional on {ℎ_𝑡+1},𝐫,𝜙 the elements of {𝑠_𝑡} are independent, we actually draw from the marginal distribution, and the acceptance rate for this algorithm is always one.

The Gibbs sampling algorithm for drawing {𝑠_𝑡}∣{ℎ_𝑡+1},𝐫,𝜙 may be described as below.(1)Specify an initial value 𝑠⁽⁰⁾=(𝑠¹⁽⁰⁾,…,𝑠^𝑇(0)).(2)Repeat for 𝑘=1,…,𝑀.(a)Repeat for 𝑡=0,…,𝑇−1.(i)Draw 𝑠^(𝑘)=1 with probability, 1𝑐_𝑡√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}𝜑𝛾+𝑑^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡},(3.59) and 𝑠^(𝑘)=−1 with probability, 1−1𝑐_𝑡√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}𝜑𝛾+𝑑^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}.(3.60)(3)Return the values {𝑠⁽¹⁾,…,𝑠^(𝑀)}.

3.3.2. Simulations of {ℎ_𝑡+1}∣𝐫, 𝜙 (Single Move Samplers)

On the other hand, the first step involves simulating from {ℎ_𝑡+1}∣𝐫,𝜙. To avoid large dependence in the chain, we use an M-H algorithm where we simulate one ℎ_𝑡+1 at a time leaving the others unchanged [20, 31]. So if (ℎ_𝑡+1)^𝑛 is the current value of the 𝑛th iteration of a Markov chain, then we draw a candidate value of the Markov chain ℎ^new_𝑡+1 by proposing it from a candidate density (proposal density) 𝑔(ℎ_𝑡+1∣(ℎ)^𝑛/𝑡+1,𝐫,𝜙) where (ℎ)^𝑛/𝑡+1={ℎ^1𝑛+1,ℎ^2𝑛+1,…,ℎ^𝑡𝑛+1,ℎ^𝑛𝑡+2,…,ℎ^𝑛𝑇+1}. We set (ℎ_𝑡+1)^𝑛+1=(ℎ_𝑡+1)^new with acceptance probabilitymin1,𝑝ℎ^new𝑡+1∣(ℎ)^𝑛/𝑡+1,𝐫,𝜙𝑔ℎ^𝑛𝑡+1∣(ℎ)^𝑛/𝑡+1,𝐫,𝜙𝑝ℎ^𝑛𝑡+1∣(ℎ)^𝑛/𝑡+1,𝐫,𝜙𝑔ℎ^new_𝑡+1∣(ℎ)^𝑛/𝑡+1,𝐫,𝜙,(3.61) where we have used the fact that𝑝(𝐡∣𝐫,𝜙)=𝑝(ℎ)_/𝑡∣𝐫,𝜙𝑝ℎ_𝑡∣(ℎ)_/𝑡,𝐫,𝜙.(3.62)

However, we may simplify further the acceptance rate. More specifically, we have that𝑝ℎ_𝑡+1∣(ℎ)_/𝑡+1,𝐫,𝜙∝𝑝ℎ_𝑡+2∣ℎ_𝑡+1,𝜙𝑝ℎ_𝑡+1∣ℎ_𝑡,𝜙𝑝𝑟_𝑡+1∣ℎ_𝑡+2,ℎ_𝑡+1,𝜙𝑝𝑟_𝑡∣ℎ_𝑡+1,ℎ_𝑡,𝜙.(3.63)

Now, since the following should hold: ℎ_𝑡+1≥𝜔+𝛽ℎ_𝑡(3.64) and similarly ℎ_𝑡+1≤𝛽⁻¹ℎ_𝑡+2−𝜔,(3.65) we have the support of the conditional distribution of ℎ_𝑡+1 given that ℎ_𝑡 is bounded from below by 𝜔+𝛽ℎ_𝑡, and the same applies to the distribution of ℎ_𝑡+2 given ℎ_𝑡+1 (lower limit corresponds to 𝑑_𝑡=0 and the upper limit to 𝑑_𝑡+1=0). This means that the range of values of ℎ_𝑡+1 compatible with ℎ_𝑡 and ℎ_𝑡+2 in the GQARCH case is bounded from above and below; that is:ℎ_𝑡+1∈𝜔+𝛽ℎ_𝑡,𝛽⁻¹ℎ_𝑡+2−𝜔.(3.66)

From the above, we understand that it makes sense to make the proposal to obey the support of the density, and so it is seen that we can simplify the acceptance rate by setting 𝑔ℎ_𝑡+1∣(ℎ)_/𝑡+1,𝐫,𝜙=𝑝ℎ_𝑡+1∣ℎ_𝑡,𝜙(3.67) appropriately truncated from above (since the truncation from below will automatically be satisfied). But the above proposal density ignores the information contained in 𝑟_𝑡+1, and so according to Fiorentini et al. [10] we can achieve a substantially higher acceptance rate if we propose from 𝑔ℎ_𝑡+1∣(ℎ)_/𝑡+1,𝐫,𝜙=𝑝ℎ_𝑡+1∣𝑟_𝑡,ℎ_𝑡,𝜙.(3.68)

A numerically efficient way to simulate ℎ_𝑡+1 from 𝑝(ℎ_𝑡+1∣𝑟_𝑡,ℎ_𝑡,𝜙) is to sample an underlying Gaussian random variable doubly truncated by using an inverse transform method. More specifically, we may draw𝜀_𝑡∣𝑟_𝑡,ℎ_𝑡,𝜙∼𝑁1−𝜑²𝑟_𝑡−𝛿ℎ_𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑²,𝜑^2𝑢ℎ^2𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑²(3.69) doubly truncated so that it remains within the following bounds:𝜀^new_𝑡∈𝛾−𝑙_𝑡,𝛾+𝑙_𝑡,(3.70) where 𝑙_𝑡=ℎ_𝑡+2−𝜔−𝛽𝜔−𝛽²ℎ_𝑡𝛽𝛼,(3.71) using an inverse transform method and then computeℎ^new_𝑡+1=𝜔+𝛼𝜀^new_𝑡−𝛾²+𝛽ℎ_𝑡,(3.72) which in turn implies a real value for 𝑑^new_𝑡+1=(ℎ_𝑡+2−𝜔−𝛽ℎ^new_𝑡+1)/𝛼 and so guarantees that ℎ^new_𝑡+1 lies within the acceptance bounds.

The inverse transform method to draw the doubly truncated Gaussian random variable first draws a uniform random number𝑢∼𝑈(0,1)(3.73) and then computes the following:𝑢=(1−𝑢)Φ⎛⎜⎜⎜⎝𝛾−𝑙^𝑡−1−𝜑²𝑟_𝑡−𝛿ℎ_𝑡/𝜑^2𝑢ℎ_𝑡+1−𝜑²+𝑢Φ⎛⎜⎜⎜⎝𝛾+𝑙_𝑡𝜑^2𝑢ℎ^2𝑡/𝜑^2𝑢ℎ_𝑡+1−𝜑²⎞⎟⎟⎟⎠−1−𝜑²𝑟_𝑡−𝛿ℎ_𝑡/𝜑^2𝑢ℎ_𝑡+1−𝜑²𝜑^2𝑢ℎ^2𝑡/𝜑^2𝑢ℎ_𝑡+1−𝜑²⎞⎟⎟⎟⎠.(3.74)

A draw is then given by𝜀^new_𝑡=Φ⁻¹𝑢.(3.75)

However, if the bounds are close to each other (the degree of truncation is small) the extra computations involved make this method unnecessarily slow, and so we prefer to use the accept-reject method where we draw 𝜀^new_𝑡∣𝑟_𝑡,ℎ_𝑡,𝜙∼𝑁1−𝜑²𝑟_𝑡−𝛿ℎ_𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑²,𝜑^2𝑢ℎ^2𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑²(3.76) and accept the draw if 𝛾−𝑙_𝑡≤𝜀^new_𝑡≤𝛾+𝑙_𝑡, and otherwise we repeat the drawing (this method is inefficient if the truncation lies in the tails of the distribution). It may be worth assessing the degree of truncation first, and, depending on its tightness, choose one simulation method or the other.

The conditional density of 𝜀^new_𝑡 will be given according to the definition of a truncated normal distribution: 𝑝𝜀^new_𝑡∣||𝜀^new_𝑡−𝛾||≤𝑙_𝑡,𝑟_𝑡,ℎ_𝑡,𝜙=1√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}𝜑𝜀^new𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}×Φ𝛾+𝑙_𝑡√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}−Φ𝛾−𝑙^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}⁻¹,(3.77) where Φ(⋅) is the cdf of the standard normal.

By using the change of variable formula, we have that the density of ℎ^new_𝑡+1 will be𝑝ℎ=^new_𝑡+1∣ℎ𝑐^new_𝑡+1∈𝜔+𝛽ℎ_𝑡,𝛽⁻¹ℎ_𝑡+2−𝜔,𝑟_𝑡,ℎ_𝑡,𝜙_new𝑡||2𝛼𝑑^new_𝑡||×Φ𝛾+𝑙^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}−Φ𝛾−𝑙^𝑡−𝜀_{𝑡/𝑟𝑟,ℎ𝑡}√𝑣_{𝑡/𝑟𝑡,ℎ𝑡}⁻¹.(3.78)

Using Bayes theorem we have that the acceptance probability will bemin1,𝑝ℎ^𝑡+2∣ℎ^new_𝑡+1,𝑟_𝑡+1,𝜙𝑝𝑟_𝑡+1∣ℎ^new_𝑡+1,𝜙𝑝𝑟_𝑡+1∣ℎ^𝑛𝑡+1,𝜙𝑝ℎ_𝑡+2∣ℎ^𝑛𝑡+1,𝑟_𝑡+1,𝜙.(3.79)

Since the degree of truncation is the same for old and new, the acceptance probability will be min1,𝑝𝑟^𝑡+1∣ℎ^new_𝑡+1𝑝𝑟_𝑡+1∣ℎ^𝑛𝑡+1𝑐𝑐^{new𝑡+1𝑛𝑡+1}𝑑^𝑛𝑡+1𝑑^new_𝑡+1,(3.80) where 𝑝(𝑟_𝑡+1∣ℎ_𝑡+1) is a mixture of two univariate normal densities so𝑟_𝑡+1∣ℎ_𝑡+1∼𝑁𝛿ℎ_𝑡+1,𝜑^2𝑢1−𝜑²ℎ_𝑡+1+1ℎ_𝑡+1.(3.81)

Hence,𝑝𝑟_𝑡+1∣ℎ^𝑛𝑡+1=12𝜋𝜑^2𝑢/1−𝜑²ℎ^𝑛𝑡+1+1ℎ^𝑛𝑡+1exp−𝑟^𝑡+1−𝛿ℎ^𝑛𝑡+1²2𝜑^2𝑢/1−𝜑²ℎ^𝑛𝑡+1+1ℎ^𝑛𝑡+1,(3.82) and the acceptance probability becomesmin⎡⎢⎢⎢⎣1,ℎ^𝑛𝑡+1ℎ^new_𝑡+1^3/2√ℎ^𝑛𝑡+2−𝜔−𝛽ℎ^𝑛𝑡+1ℎ^𝑛𝑡+2−𝜔−𝛽ℎ^new_𝑡+1𝜅ℎ^new_𝑡+1𝜅ℎ^𝑛𝑡+1⎤⎥⎥⎥⎦,(3.83) where 𝜅ℎ^𝑖𝑡+1=exp⎡⎢⎢⎢⎢⎢⎢⎣−𝑟^𝑡+1−𝛿ℎ^𝑖𝑡+1²2𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1+1ℎ^𝑖𝑡+1−𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1+12𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1²⎛⎜⎜⎜⎜⎜⎜⎝⎛⎜⎜⎜⎜⎝𝛾−𝑟_𝑡+1−𝛿ℎ^𝑖𝑡+1𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1+1⎞⎟⎟⎟⎟⎠²+ℎ^𝑟𝑡+2−𝜔−𝛽ℎ^𝑖𝑡+1×⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣1+exp⎡⎢⎢⎢⎢⎣2𝜑_2𝑢𝛼⎞⎟⎟⎟⎟⎟⎟⎠⎤⎥⎥⎥⎥⎥⎥⎦1−𝜑²ℎ^𝑖𝑡+1+1𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1²⎛⎜⎜⎜⎜⎝𝛾−𝑟_𝑡+1−𝛿ℎ^𝑖𝑡+1𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1+1⎞⎟⎟⎟⎟⎠ℎ^𝑛𝑡+2−𝜔−𝛽ℎ^𝑖𝑡+1𝛼⎤⎥⎥⎥⎥⎦exp⎡⎢⎢⎢⎢⎣𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1+1𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1²⎛⎜⎜⎜⎜⎝𝛾−𝑟_𝑡+1−𝛿ℎ^𝑖𝑡+1𝜑^2𝑢1−𝜑²ℎ^𝑖𝑡+1+1⎞⎟⎟⎟⎟⎠ℎ^𝑛𝑡+2−𝜔−𝛽ℎ^𝑖𝑡+1𝛼⎤⎥⎥⎥⎥⎦⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦.(3.84)

Overall the MCMC of {ℎ_𝑡+1}∣𝐫,𝜙 includes the following steps.(1)Specify an initial value {ℎ⁽⁰⁾}.(2)Repeat for 𝑛=1,…,𝑀.(a)Repeat for 𝑡=0,…,𝑇−1.(i)Use an inverse transform method to simulate 𝜀^new_𝑡∣𝑟_𝑡,ℎ_𝑡,𝜙∼𝑁1−𝜑²𝑟_𝑡−𝛿ℎ_𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑²,𝜑^2𝑢ℎ^2𝑡𝜑^2𝑢ℎ_𝑡+1−𝜑²(3.85) doubly truncated.(ii)Calculate ℎ^new_𝑡+1=𝜔+𝛼𝜀^𝑐𝑡−𝛾²+𝛽ℎ_𝑡.(3.86) Steps (2)(a)(i) and (2)(a)(ii) are equivalent to draw ℎ_𝑡+1^new∼𝑝ℎ^new_𝑡+1∣𝑟_𝑡,ℎ_𝑡,𝜙(3.87) appropriately truncated.(iii)Calculate 𝛼_𝑟=min1,𝑝𝑟^𝑡+1∣ℎ^new_𝑡+1𝑝𝑟_𝑡+1∣ℎ^𝑛𝑡+1𝑐𝑐^{new𝑡+1𝑛𝑡+1}𝑑^𝑛𝑡+1𝑑^new_𝑡+1.(3.88)(iv)Set ℎ_𝑡+1^𝑛+1=ℎ^𝑡+1^newifUnif(0,1)≤𝛼_𝑟ℎ_𝑡+1^𝑛o