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

Authors: Sofia Anyfantaki; Antonis Demos
Journal: Journal of Probability and Statistics
 (2011)

---

## 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 onlyOT computational operations, where T 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.


---

## Body

## 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 O(T) algorithms, where T 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 Section2.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 Section2, 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 O(T2) to O(T). A comparison of the efficient (O(T) calculations) and the inefficient (O(T2) 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:(2.1)rt=δtht+ɛt,ɛt=ztht1/2,
where
(2.2)δt=(1-φ)δ+φδt-1+φuut,(2.3)ht=ω+α(ɛt-1-γ)2+βht-1,zt∽i.i.d.N(0,1), ut∽i.i.d.N(0,1), ut, zt are independent for allt′s, and where {rt}t=1T are the observed excess returns, T is the sample size, {δt}t=1T is an unobserved AR(1) process independent (with δ0=δ) of {ɛt}t=1T, and {ht}t=1T is the conditional variance (with h0 equal to the unconditional variance and ɛ0=0) which is supposed to follow a GQARCH(1,1). It is obvious that δt is the market price of risk (see, e.g., Merton [1] Glosten at al. [2]). Let us call ℱt-1 the sequence of natural filtrations generated by the past values of {ɛt} and {rt}.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 varianceht, needed for the moments of rt. The proof of the following lemma is based on raising ht to the appropriate power, in (2.3), and taking into account that E(zt4)=3, E(zt6)=15 and E(zt8)=105.Lemma 2.2.
If105α4+60α3β+18α2β2+4αβ3+β4<1, the first four moments of the conditional variance of (2.3) are given by
(2.4)E(ht)=ω+αγ21-(α+β),E(ht2)=(ω+αγ2)(ω+αγ2)(1+α+β)+4α2γ2(1-3α2-2αβ-β2)(1-α-β),E(ht3)=(ω+αγ2)3+3(ω+αγ2)[(ω+αγ2)(α+β)+4α2γ2]E(ht)1-β3-15α3-9α2β-3αβ2+3[(ω+αγ2)(3α2+β2+2βα)+4α2γ2(3α+β)]E(ht2)1-β3-15α3-9α2β-3αβ2,E(ht4)=(ω+αγ2)2(ω+αγ2)2+4[(ω+αγ2)(α+β)+6α2γ2]E(ht)1-105α4-60α3β-18α2β2-4αβ3-β4+6(ω+αγ2)2(3α2+β2+2βα)+8[(ω+αγ2)(3α+β)+α2γ2]α2γ21-105α4-60α3β-18α2β2-4αβ3-β4E(ht2)+4(ω+αγ2)(15α3+9α2β+3αβ2+β3)+6α2γ2(15α2+6αβ+β2)1-105α4-60α3β-18α2β2-4αβ3-β4E(ht3).Consequently, the moments ofrt 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.5)E(rt)=δE(ht),E(rt2)=(δ2+φu21-φ2)E(ht2)+E(ht),
whereas the skewness and kurtosis coefficients are
(2.6)Sk(rt)=S(rt)Var⁡1.5(rt),kurt(rt)=κVar⁡2(rt),
where
(2.7)S(rt)=δ(δ2+3φu21-φ2)E(ht3)+3δE(h2)+2δ3E3(ht)-3δ[(δ2+φu21-φ2)E(ht2)+E(ht)]E(ht),κ=(δ4+6δ2φu21-φ2+3φu4(1-φ2)2)E(ht4)+3δ2[2-δ2E(ht)]E3(ht)+6(δ2+φu21-φ2)E(ht3)-4δ2(δ2+3φu21-φ2)E(ht3)E(ht)+{6δ2[(δ2+φu21-φ2)E(h)-2]E(h)+3}E(ht2),
and E(ht), E(ht2), E(ht3),and E(ht4) are given in Lemma 2.2.In terms of stationarity, the process{rt} is 4th-order stationary if and only if
(2.8)|φ|<1,105α4+60α3β+18α2β2+4αβ3+β4<1.
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 rt’s (volatility clustering), as well as the correlation of rt2 and rt-1 (leverage effect). The dynamic moments of the conditional variance and those between the conditional variance ht and the error ɛt are given in the following two lemmas (for a proof see Appendix B).Lemma 2.4.
Under the assumption of Lemma2.2, one has that
(2.9)Cov(ht,ht-k)=(α+β)kV(ht),Cov(ht2,ht-k)=Ak[E(ht3)-E(ht2)E(ht)]+(α+β)k-Ak(α+β)-ABV(ht),Cov(ht,ht-k2)=(α+β)k[E(ht-13)-E(ht-12)E(ht-1)],Cov(ht2,ht-k2)=AkV(ht2)+B(α+β)k-Ak(α+β)-A[E(ht3)-E(ht2)E(ht)],
where A=3α2+β2+2αβ and B=2[2α2γ2+(ω+αγ2)(α+β)].Lemma 2.5.
(2.10)Cov(ht,ɛt-k)=E(htɛt-k)=-2αγ(α+β)k-1E(ht),Cov(htht-k,ɛt-k)=E(htht-kɛt-k)=-2αγ(α+β)k-1E(ht-12),Cov(ht,ɛt-k2)=(α+β)kV(ht)+2α(α+β)k-1E(ht2),Cov(ht2,ɛt-k)=E(ht2ɛt-k)=-4αγ(3α+β)Ak-1E(ht2)-4αγ[(ω+αγ2)(α+β)k-Ak(α+β)-A+2α2γ2(α+β)k-1-Ak-1(α+β)-A]E(ht),Cov(ht2,ɛt-k2)=Ak[E(ht3)-E(ht2)E(ht)]+B(α+β)k-Ak(α+β)-AV(ht)+4α(3α+β)Ak-1E(ht3)+[2αB(α+β)k-1-Ak-1(α+β)-A+4(2α2γ2+(ω+αγ2))Ak-1]E(ht2),Cov(ht2ht-k,ɛt-k)=E(ht2ht-kɛt-k)=-4αγAk-1[(3α+β)E(ht3)+(ω+αγ2)E(ht2)]-2αγB(α+β)k-1-Ak-1(α+β)-AE(ht2),
where A and B 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
(2.11)γk=Cov(rt,rt-k)=δ2Cov(ht,ht-k)+δE(htɛt-k)+φkφu21-φ2E(htht-k),
and the covariance of squares’ levels and the autocovariance of squares are
(2.12)Cov(rt2,rt-k)=E(δt2δt-k)Cov(ht2,ht-k)+Cov(δt2,δt-k)E(ht2)E(ht-k)+E(δt-k)Cov(ht,ht-k)+E(δt2)E(ht2ɛt-k)+E(htɛt-k),Cov(rt2,rt-k2)=E2(δt2)Cov(ht2,ht-k2)+Cov(δt2,δt-k2)E(ht2ht-k2)+E(δt2)Cov(ht2,ɛt-k2)+2E(δt2δt-k)E(ht2ht-kɛt-k)+E(δt2)Cov(ht,ht-k2)+Cov(ht,ɛt-k2)+2δE(htht-kɛt-k),
where all needed covariances and expectations of the right-hand sizes are given in Lemmas 2.4 and 2.5,
(2.13)Cov(δt2,δt-k)=Cov(δt,δt-k2)=2φkδφu21-φ2,Cov(δt2,δt-k2)=4φkδ2φu21-φ2+2φ2kφu4(1-φ2)2.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 (negativeCov(ht,ɛt-k)) can be accommodated, and so forth. Furthermore, the model can accommodate negative autocorrelations, γk, 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.


## 2.1. Theoretical Properties

Let us consider first the moments of the conditional varianceht, needed for the moments of rt. The proof of the following lemma is based on raising ht to the appropriate power, in (2.3), and taking into account that E(zt4)=3, E(zt6)=15 and E(zt8)=105.Lemma 2.2.
If105α4+60α3β+18α2β2+4αβ3+β4<1, the first four moments of the conditional variance of (2.3) are given by
(2.4)E(ht)=ω+αγ21-(α+β),E(ht2)=(ω+αγ2)(ω+αγ2)(1+α+β)+4α2γ2(1-3α2-2αβ-β2)(1-α-β),E(ht3)=(ω+αγ2)3+3(ω+αγ2)[(ω+αγ2)(α+β)+4α2γ2]E(ht)1-β3-15α3-9α2β-3αβ2+3[(ω+αγ2)(3α2+β2+2βα)+4α2γ2(3α+β)]E(ht2)1-β3-15α3-9α2β-3αβ2,E(ht4)=(ω+αγ2)2(ω+αγ2)2+4[(ω+αγ2)(α+β)+6α2γ2]E(ht)1-105α4-60α3β-18α2β2-4αβ3-β4+6(ω+αγ2)2(3α2+β2+2βα)+8[(ω+αγ2)(3α+β)+α2γ2]α2γ21-105α4-60α3β-18α2β2-4αβ3-β4E(ht2)+4(ω+αγ2)(15α3+9α2β+3αβ2+β3)+6α2γ2(15α2+6αβ+β2)1-105α4-60α3β-18α2β2-4αβ3-β4E(ht3).Consequently, the moments ofrt 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.5)E(rt)=δE(ht),E(rt2)=(δ2+φu21-φ2)E(ht2)+E(ht),
whereas the skewness and kurtosis coefficients are
(2.6)Sk(rt)=S(rt)Var⁡1.5(rt),kurt(rt)=κVar⁡2(rt),
where
(2.7)S(rt)=δ(δ2+3φu21-φ2)E(ht3)+3δE(h2)+2δ3E3(ht)-3δ[(δ2+φu21-φ2)E(ht2)+E(ht)]E(ht),κ=(δ4+6δ2φu21-φ2+3φu4(1-φ2)2)E(ht4)+3δ2[2-δ2E(ht)]E3(ht)+6(δ2+φu21-φ2)E(ht3)-4δ2(δ2+3φu21-φ2)E(ht3)E(ht)+{6δ2[(δ2+φu21-φ2)E(h)-2]E(h)+3}E(ht2),
and E(ht), E(ht2), E(ht3),and E(ht4) are given in Lemma 2.2.In terms of stationarity, the process{rt} is 4th-order stationary if and only if
(2.8)|φ|<1,105α4+60α3β+18α2β2+4αβ3+β4<1.
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 rt’s (volatility clustering), as well as the correlation of rt2 and rt-1 (leverage effect). The dynamic moments of the conditional variance and those between the conditional variance ht and the error ɛt are given in the following two lemmas (for a proof see Appendix B).Lemma 2.4.
Under the assumption of Lemma2.2, one has that
(2.9)Cov(ht,ht-k)=(α+β)kV(ht),Cov(ht2,ht-k)=Ak[E(ht3)-E(ht2)E(ht)]+(α+β)k-Ak(α+β)-ABV(ht),Cov(ht,ht-k2)=(α+β)k[E(ht-13)-E(ht-12)E(ht-1)],Cov(ht2,ht-k2)=AkV(ht2)+B(α+β)k-Ak(α+β)-A[E(ht3)-E(ht2)E(ht)],
where A=3α2+β2+2αβ and B=2[2α2γ2+(ω+αγ2)(α+β)].Lemma 2.5.
(2.10)Cov(ht,ɛt-k)=E(htɛt-k)=-2αγ(α+β)k-1E(ht),Cov(htht-k,ɛt-k)=E(htht-kɛt-k)=-2αγ(α+β)k-1E(ht-12),Cov(ht,ɛt-k2)=(α+β)kV(ht)+2α(α+β)k-1E(ht2),Cov(ht2,ɛt-k)=E(ht2ɛt-k)=-4αγ(3α+β)Ak-1E(ht2)-4αγ[(ω+αγ2)(α+β)k-Ak(α+β)-A+2α2γ2(α+β)k-1-Ak-1(α+β)-A]E(ht),Cov(ht2,ɛt-k2)=Ak[E(ht3)-E(ht2)E(ht)]+B(α+β)k-Ak(α+β)-AV(ht)+4α(3α+β)Ak-1E(ht3)+[2αB(α+β)k-1-Ak-1(α+β)-A+4(2α2γ2+(ω+αγ2))Ak-1]E(ht2),Cov(ht2ht-k,ɛt-k)=E(ht2ht-kɛt-k)=-4αγAk-1[(3α+β)E(ht3)+(ω+αγ2)E(ht2)]-2αγB(α+β)k-1-Ak-1(α+β)-AE(ht2),
where A and B 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
(2.11)γk=Cov(rt,rt-k)=δ2Cov(ht,ht-k)+δE(htɛt-k)+φkφu21-φ2E(htht-k),
and the covariance of squares’ levels and the autocovariance of squares are
(2.12)Cov(rt2,rt-k)=E(δt2δt-k)Cov(ht2,ht-k)+Cov(δt2,δt-k)E(ht2)E(ht-k)+E(δt-k)Cov(ht,ht-k)+E(δt2)E(ht2ɛt-k)+E(htɛt-k),Cov(rt2,rt-k2)=E2(δt2)Cov(ht2,ht-k2)+Cov(δt2,δt-k2)E(ht2ht-k2)+E(δt2)Cov(ht2,ɛt-k2)+2E(δt2δt-k)E(ht2ht-kɛt-k)+E(δt2)Cov(ht,ht-k2)+Cov(ht,ɛt-k2)+2δE(htht-kɛt-k),
where all needed covariances and expectations of the right-hand sizes are given in Lemmas 2.4 and 2.5,
(2.13)Cov(δt2,δt-k)=Cov(δt,δt-k2)=2φkδφu21-φ2,Cov(δt2,δt-k2)=4φkδ2φu21-φ2+2φ2kφu4(1-φ2)2.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 (negativeCov(ht,ɛt-k)) can be accommodated, and so forth. Furthermore, the model can accommodate negative autocorrelations, γk, 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 theE 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 π(ϕt/ϕ/t). By Bayes' theorem, we have π(ϕt/ϕ/t)∝π(ϕtϕ/t), 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 {ϕt}t=1p, where p 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 ht (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 O(T2) 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ɛt and δt. On the other hand, the conditional log-likelihood function of our model assuming that δt were observed would be the following:
(3.1)l(r,δ∣ϕ,F0)=ln⁡p(r∣δ,ϕ,F0)+ln⁡p(δ∣ϕ,F0)=-Tln⁡2π-12∑t=1Tln⁡ht-12∑t=1T(ɛt)2ht-T2ln⁡(φu2)-12∑t=1T(δt-δ(1-φ)-φδt-1)2φu2,
where r=(r1,…,rT)′,δ=(δ1,…,δT)′, and h=(h1,…,hT)′.However, theδt’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 ϕ(n+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, E(ℓ(·)∣r,ϕ(n),ℱ0) with respect to ϕ keeping ϕ(n) fixed.TheE step, thus, requires the expectation of the complete log-likelihood. For our model, this is given by: (3.2)E(l(⋅)∣r,ϕ(n),F0)=-Tln⁡2π-T2ln⁡(φu2)-12∑t=1TE(ln⁡ht∣r,ϕ(n),F0)-12∑t=1TE((ɛt)2ht∣r,ϕ(n),F0)-12∑t=1TE((δt-δ(1-φ)-φδt-1)2φu2∣r,ϕ(n),F0).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:(3.3)SEMl=-Tln⁡2π-T2ln⁡(φu2)-121M∑i=1M∑t=1Tln⁡ht(i)-121M∑i=1M∑t=1T(ɛt(i))2ht(i)-T2(1-φ)2δ2φu2-121φu21M∑i=1M∑t=1T(δt(i))2+(1-φ)δφu21M∑i=1M∑t=1Tδt(i)+φφu21M∑i=1M∑t=1Tδt(i)δt-1(i)+(1-φ)φδφu21M∑i=1M∑t=1Tδt-1(i)-φ22φu21M∑i=1M∑t=1T(δt-1(i))2.
Consequently, we need to obtain the following quantities: (1/M)∑i=1M∑t=1Tlnht(i),  (1/M)∑i=1M∑t=1T((ɛt(i))2/ht(i)), (1/M)∑i=1M∑t=1Tδt(i), (1/M)∑i=1M∑t=1Tδt-1(i),  (1/M)∑i=1M∑t=1Tδt(i)δt-1(i) and (1/M)∑i=1M∑t=1T(δt(i))2, (1/M)∑i=1M∑t=1Tδt-12(i), where M is the number of simulations.Thus, to classically estimate our model by using an SEM algorithm, the basic problem is to sample fromh∣ϕ,r,ℱ0 where ϕ is the vector of the unknown parameters and also sample from δ∣ϕ,r,ℱ0.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, ɛt, to the price of risk, δt, and vice versa. One possible solution is to fix ω such that E(ht) is 1 or to set φu to a specific value. In fact in an earlier version of the paper, we fixed φu 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δ∣ϕ,r,ℱ0.


### 3.2. Simulation-Based Bayesian Inference

In our problem, the key issue is that the likelihood function of the samplep(r∣ϕ,ℱ0) is intractable which precludes the direct analysis of the posterior density p(ϕ∣r,ℱ0). This problem may be overcome by focusing instead on the posterior density of the model using Bayes’ rule:(3.4)p(ϕ,δ∣r)∝p(ϕ,δ)p(r∣ϕ,δ)∝p(δ∣ϕ)p(ϕ)p(r∣ϕ,δ),
where(3.5)ϕ=(δ,φ,φu,α,β,γ,ω)′.Now,(3.6)p(δ∣ϕ)=∏t=1Tp(δtδ,φ,φu2)=∏t=1T12πφu2exp⁡(-(δt-δ(1-φ)-φδt-1)22φu2).On the other hand,(3.7)p(r∣ϕ,δ)=∏t=1Tp(rt{rt-1},δ,ϕ)=∏t=1T12πhtexp⁡(-ɛt22ht)
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 {δt}(0),ϕ(0), we draw {δt}(1) from p({δt}(1)∣r,ϕ(0)) and then ϕ(1) from p(ϕ(1)∣{δt}(1),r). Iterating these steps, we finally get ({δt}(i),ϕ(i))i=1M, and under mild conditions it is shown that the distribution of the sequence converges to the joint posterior distribution p(ϕ,δ∣r).The above simulation procedure may be carried out by first dividing the parameters into two blocks:(3.8)ϕ1=(δ,φ,φu2),ϕ2=(α,β,γ,ω),
Then the algorithm is described as follows.(1)
Initializeϕ.(2)
Draw fromp(δt∣δ≠t,r,ϕ).(3)
Draw fromp(ϕ∣δ,r) in the following blocks:
(i)
draw fromp(ϕ1∣δ,r) using the Gibbs sampling. This is updated in one block;(ii)
draw fromp(ϕ2∣r) 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(δ,φu2,φ) is given by:(3.9)p(δ,φu2,φ)=p(δ∣φu2)p(φu2)p(φ),
which means that δ,φu2 is a priory independent of φ.Also the following holds for the prior distributions of the parameter subvectorϕ1:
(3.10)p(δ∣φu2)~N(δpr,φu2σδpr2),p(φu2)~IG(v02,d02),p(φ)~N(φ0,σφ02)Iφ,
where Iφ ensures that φ lies outside the unit circle, IG is the inverted gamma distribution, and the hyperparameters v0,d0,δpr,σδpr2,φ0,σφ02 have to be defined.Now, the joint posterior is proportional to(3.11)p(δ,φ,φu2∣r,δ)∝∏t=1T12πφu2exp⁡{-(δt-(1-φ)δ-φδt-1)22φu2}×N(δpr,φu2σδpr2)×IG(v02,d02)×N(φ0,σφ02)Iφ.
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
(3.12)δt*=δt-φδt-1,δt*∣Ft-1~N((1-φ)δ,φu2),
or, otherwise,
(3.13)δt*=(1-φ)δ+vt,vt~N(0,φu2).
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
(3.14)δ∣δ,ϕ,φu2~N(δ̃,φu2σ̃δ2),
where
(3.15)δ̃=σ̃δ2(δprσδpr2+(1-φ)∑t=1Tδt*),σ̃δ2=(1σδpr2+(1-φ)2)-1.
Hence, the generation ofδ is completed, and we may turn on the generation of the other parameters.Generation ofφu2
For the generation ofφu2 and using [27] notation, we have the following.
Proposition 3.2.
The proposal distribution ofφu2 is
(3.16)φu2∣δ,ϕ,δ~IG(T-v02,d0+Q+d2),
where
(3.17)Q=(δ-δpr)2σδpr-2,d=∑t=2T[δt*-δ(1-φ)]2.
Finally, we turn on the generation ofφ.Generation ofφ
For the generation ofφ, we follow again Chib [27] and write
(3.18)δt=(1-φ)δ-φδt-1+vt,vt~N(0,φu2).
We may now state the following proposition (see Chib [27] for a proof).
Proposition 3.3.
The proposal distribution ofφ is
(3.19)φ2∣δ,δ,φu2~N(φ̃,σ̃φ2),
where
(3.20)φ̃=σ̃φ2(σφ0-2φ0+φu-2∑t=1T(δt-1-δ)(δt-δ)),σ̃φ-2=σφ0-2+φu-2∑t=1T(δt-1-δ)2.
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:(3.21)p(α)~N(μα,Σα)Iα,p(β)~N(μβ,σβ2)Iβ
and similarly
(3.22)p(γ)~N(μγ,σγ2),
where α=(ω,α)′,Iα,Iβ 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:(3.23)p(α,β,γ∣r)∝∏t=1T12πhtexp⁡(-ɛt22ht)×N(μα,Σα)Iα×N(μβ,σβ2)Iβ×N(μγ,σγ2).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]:
(3.24)ɛt2=ω+α(ɛt-1-γ)2+βɛt-12+wt-βwt-1,
where wt=ɛt2-ht with wt~N(0,2ht2).Then the corresponding approximated likelihood is written as(3.25)p(ɛ2∣δ,r,ϕ2)=∏t=1T12htπexp⁡[-(ɛt2-ω-α(ɛt-1-γ)2-βɛt-12+βwt-1)24ht2],
and the generation of α,β,γ is based on the above likelihood where we update {ht} 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 wt in (3.32), below, can be written as a linear function of α:
(3.26)wt=ɛ¯t2-ζtα,
where ζt=[ιt̃,ɛ̂t2] with
(3.27)ɛ¯t2=ɛ̃t2-βɛ̃t-12,ɛ̃t2=ɛt2+βɛ̃t-12,ɛ̂t2=(ɛt-1-γ)2+βɛ̂t-12,ιt̃=1+βι̃t-1.
Now, let the two following vectors be(3.28)Yα=[ɛ¯12,…,ɛ¯T2]′,Xα=[ζ1′,…,ζT′]′.
Then the likelihood function of the approximated model is rewritten as(3.29)p(ɛ2∣r,δ,ϕ2)=∏t=1T12π(2ht2)exp⁡[-(ɛ¯t2-ζtα)22(2ht2)].
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
(3.30)α∣Y,X,Σ,ϕ2-α~N(μ̂α,Σ̂α)Iα,
where μ̂α=Σ̂α(Xα′Λ-1Yα+Σα-1μα), Σ̂α=(Xα′Λ-1Xα+Σα-1)-1, and Λ=diag(2h12,…,2hT2). Iα imposes the restriction that α>0 and α+β<1.
Hence a candidateα̃ is sampled from this proposal density and accepted with probability:
(3.31)min⁡{p(α̃,β,γ∣δ,r)q(α*∣α̃,β,γ,δ,r)p(α*,β,γ∣δ,r)q(α̃∣α*,β,γ,δ,r),1},
where α* is the previous draw.
Similar procedure is used for the generation ofβ and γ.Generation ofβ
Following Nakatsuma [29], we linearize wt by the first-order Taylor expansion
(3.32)wt(β)≈wt(β*)+ξt(β*)(β-β*),
where ξt is the first-order derivative of wt(β) evaluated at β* the previous draw of the M-H sampler.
Define as(3.33)rt=wt(β*)+gt(β*)β*,
where gt(β*)=-ξt(β*) which is computed by the recursion:
(3.34)gt=ɛt-12-wt-1+β*gt-1,ξt=0 for t≤0 [30]. Then,
(3.35)wt(β)≈rt-gt(β)β.
Let the following two vectors be(3.36)Yβ=[r1,…,rT]′,Xβ=[g1,…,gT]′.
Then, the likelihood function of the approximated model is rewritten as(3.37)p(ɛ2∣δ,r,ϕ2)=∏t=1T12π(2ht2)exp⁡[-{wt(β*)+ξt(β*)(β-β*)}22(2ht2)].
We have the following proposal distribution forβ (for a proof see Nakatsuma [29] or Ardia [30]).Proposition 3.5.
The proposal distribution forβ is
(3.38)β∣Y,X,σβ2,ϕ2-β~N(μ̂β,σ̂β2)Iβ,
where μ̂β=σ̂β2(Xβ′Λ-1Yβ+(μβ/σβ2)), σ̂β2=(Xβ′Λ-1Xβ+(1/σβ2))-1, and Λ=diag(2h14,…,2hT4). Iβ imposes the restriction that β>0 and α+β<1.
Hence, a candidateβ̃ is sampled from this proposal density and accepted with probability:
(3.39)min⁡{p(β̃,α,γ∣δ,r)q(β*∣β̃,α,γ,δ,r)p(β*,α,γ∣δ,r)q(β̃∣β*,α,γ,δ,r),1}Finally, we explain the generation ofγ.Generation ofγ
As withβ, we linearize wt by a first-order Taylor, expansion at a point γ* the previous draw in the M-H sampler. In this case,
(3.40)rt=wt(γ*)-gt(γ*)γ*,
where gt(γ*)=-ξt(γ*) which is computed by the recursion:
(3.41)gt=-2α(ɛt-1-γ*)+βgt-1,
and gt=0 for t≤0.
Then(3.42)wt(γ)≈rt-gtγ.
Let again(3.43)Yγ=[r1,…,rT]′,Xγ=[g1,…,gT]′
and the likelihood function of the approximated model is rewritten as
(3.44)p(ɛ2∣δ,r,ϕ2)=∏t=1T12π(2ht2)exp⁡[-{wt(γ*)-gt(γ*)γ*}22(2ht2)].
Thus, we have the following proposal distribution forγ (for proof see Nakatsuma [29] and Ardia [30]).Proposition 3.6.
The proposal distribution ofγ is
(3.45)γ∣Y,X,σγ2,ϕ2-γ~N(μ̂γ,σ̂γ2),
where μ̂γ=σ̂γ2(Xγ′Λ-1Yγ+(μγ/σγ2)), σ̂γ2=(Xγ′Λ-1Xγ+(1/σγ2))-1, and Λ=diag(2h14,…,2hT4). A candidate γ̃ is sampled from this proposal density and accepted with probability:
(3.46)min⁡{p(γ̃,α,β∣δ,r)q(γ*∣γ̃,α,β,δ,r)p(γ*,α,β∣δ,r)q(γ̃∣γ*,α,β,δ,r),1}.The algorithm described above is a special case of a MCMC algorithm, which converges as it iterates, to draws from the required densityp(ϕ,δ∣r). 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δ∣ϕ,r,ℱ0.


#### 3.2.3. MCMC Simulation ofɛ∣ϕ,r,ℱ0

For a given set of parameter values and initial conditions, it is generally simpler to simulate{ɛt} for t=1,…,T and then compute {δt}t=1T than to simulate {δt}t=1T directly. For that matter, we concentrate on simulators of ɛt given r and ϕ. We set the mean and the variance of ɛ0 equal to their unconditional values, and, given that ht is a sufficient statistic for ℱt-1 and the unconditional variance is a deterministic function of ϕ, ℱ0 can be eliminated from the information set without any information loss.Now sampling fromp(ɛ∣r,ϕ)∝p(r∣ɛ,ϕ)p(ɛ∣ϕ) is feasible by using an M-H algorithm where we update each time only one ɛt leaving all the other unchanged [20]. In particular, let us write the nth  iteration of a Markov chain as ɛtn. Then we generate a potential new value of the Markov chain ɛtnew by proposing from some candidate density g(ɛt∣ɛ∖tn,r,ϕ) where ɛ∖tn={ɛ1n+1,…,ɛt-1n+1,ɛt+1n,…,ɛTn} which we accept with probability(3.47)min⁡[1,p(ɛtnew∣ɛ∖tn,r,ϕ)g(ɛtnew∣ɛ∖tn,r,ϕ)p(ɛtn∣ɛ∖tn,r,ϕ)g(ɛtn∣ɛ∖tn,r,ϕ)].If it is accepted then, we setɛtn+1=ɛtnew and otherwise we keep ɛtn+1=ɛtn. 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:
(3.48)p(ɛtnew∣ɛ∖tn,r,ϕ)p(ɛtn∣ɛ∖tn,r,ϕ)=p(rt∣ɛtnew,htnew,t,ϕ)p(ɛtnew∣htnew,t,ϕ)p(rt∣htn,t,ϕ)p(rt∣htnew,t,ϕ)p(rt∣ɛtn,htn,t,ϕ)p(ɛtnew∣htn,t,ϕ)*∏s=t+1Tp(rs∣ɛsr,hsnew,t,ϕ)p(ɛsr∣hsnew,t,ϕ)p(rs∣hsn,t,ϕ)p(rs∣hsnew,t,ϕ)p(rs∣ɛsn,hsn,t,ϕ)p(ɛsn∣hsn,t,ϕ)=p(rt∣ɛtnew,htnew,t,ϕ)p(ɛtnew∣htnew,t,ϕ)p(rt∣ɛtn,htn,t,ϕ)p(ɛtnew∣htn,t,ϕ)*∏s=t+1Tp(rs∣ɛsn,hsnew,t,ϕ)p(ɛsn∣hsnew,t,ϕ)p(rs∣ɛsn,hsn,t,ϕ)p(ɛsn∣hsn,t,ϕ),
where for s=t+1,…,T,(3.49)hsnew,t=V(ɛs∣ɛs-1n,ɛs-2n,…,ɛt+1n,ɛtnew,ɛt-1n+1,…,ɛ1n+1),hsn,t=V(ɛs∣ɛs-1n,ɛs-2n,…,ɛt+1n,ɛtn,ɛt-1n+1,…,ɛ1n+1)
while(3.50)htnew,t=htn,t.Nevertheless, each time we revise oneɛt, we have also to revise T-t conditional variances because of the recursive nature of the GARCH model which makes hsnew,t depend upon ɛtnew for s=t+1,…,T. And since t=1,…,T, it is obvious that we need to calculate T2 normal densities, and so this algorithm is O(T2). 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 O(T) 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 O(T) calculations. Following their transformation, we augment the state vector with the variables ht+1 and then sample the joint Markov process {ht+1,st}∣r,ϕ∈ℱt where
(3.51)st=sign⁡(ɛt-γ),
so that st=±1 with probability one. The mapping is one to one and has no singularities. More specifically, if we know {ht+1} and φ, then we know the value of(3.52)(ɛt-γ)2=ht+1-ω-βhtα∀t≥1.Hence the additional knowledge of the signs of(ɛt-γ) would reveal the entire path of {ɛt} so long as h0 (which equals the unconditional value in our case) is known, and, thus, we may now reveal also the unobserved random variable {δt}∣r,ϕ,{ht+1}.Now we have to sample from(3.53)p({st,ht+1}∣r,ϕ)∝∏t=1Tp(st∣ht+1,ht,ϕ)p(ht+1∣ht,ϕ)p(rt∣st,ht,ht+1,ϕ),
where the second and the third term come from the model, and the first comes from the fact that ɛt∣ℱt-1~N(0,ht) but ɛt∣{ht+1},ℱt-1 takes values(3.54)ɛt=γ±dt,
where
(3.55)dt=ht+1-ω-βhtα.From the above, it is seen that we should first simulate{ht+1}∣r,ϕ since we do not alter the volatility process when we flip from st=-1 to st=1 (implying that the signs do not cause the volatility process), but we do alter ɛt and then simulate {st}∣{ht+1},r,ϕ. The second step is a Gibbs sampling scheme whose acceptance rate is always one and also conditional on {ht+1},r,ϕ the elements of {st} 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{st}∣{ht+1},r, ϕ

First, we see how to sample from{st}∣{ht+1},r,ϕ. To obtain the required conditionally Bernoulli distribution, we establish first some notation. We have the following (see Appendix A):
(3.56)ct=1vt∣rt,ht[φ(γ+dt-ɛt∣rt,htvt∣rt,ht)+φ(γ-dt-ɛt∣rt,htvt∣rt,ht)],ɛt∣rt,ht=E(ɛt∣rt,ht)=(1-φ2)(rt-δht)φu2ht+1-φ2,vt∣rt,ht=Var⁡(ɛt∣rt,ht)=φu2ht2φu2ht+1-φ2.Using the above notation, we see that the probability of drawingst=1 conditional on {ht+1} is equal to the probability of drawing ɛt=γ+dt conditional on ht+1,ht,rt,ϕ, where dt is given by (3.70), which is given by(3.57)p(st=1∣{ht+1},r,ϕ)=p(ɛt=γ+dt∣ht+1,ht,rt,ϕ)=1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht).Similarly for the probability of drawingst=-1. Both of these quantities are easy to compute; for example,(3.58)φ(γ+dt-ɛt/rr,htvt/rt,ht)=12πexp⁡{-12(γ+dt-ɛt/rr,htvt/rt,ht)2}
and so we may simulate {st}∣{ht+1},r,ϕ using a Gibbs sampling scheme. Specifically, since conditional on {ht+1},r,ϕ the elements of {st} 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{st}∣{ht+1},r,ϕ may be described as below.(1)
Specify an initial values(0)=(s1(0),…,sT(0)).(2)
Repeat fork=1,…,M.
(a)
Repeat fort=0,…,T-1.
(i)
Draws(k)=1 with probability,
(3.59)1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht),
and s(k)=-1 with probability,
(3.60)1-1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht).(3)
Return the values{s(1),…,s(M)}.


#### 3.3.2. Simulations of{ht+1}∣r, ϕ (Single Move Samplers)

On the other hand, the first step involves simulating from{ht+1}∣r,ϕ. To avoid large dependence in the chain, we use an M-H algorithm where we simulate one ht+1 at a time leaving the others unchanged [20, 31]. So if (ht+1)n is the current value of the nth iteration of a Markov chain, then we draw a candidate value of the Markov chain ht+1new by proposing it from a candidate density (proposal density) g(ht+1∣(h)/t+1n,r,ϕ) where (h)/t+1n={h1n+1,h2n+1,…,htn+1,ht+2n,…,hT+1n}. We set (ht+1)n+1=(ht+1)new with acceptance probability(3.61)min⁡[1,p(ht+1new∣(h)/t+1n,r,ϕ)g(ht+1n∣(h)/t+1n,r,ϕ)p(ht+1n∣(h)/t+1n,r,ϕ)g(ht+1new∣(h)/t+1n,r,ϕ)],
where we have used the fact that(3.62)p(h∣r,ϕ)=p((h)/t∣r,ϕ)p(ht∣(h)/t,r,ϕ).However, we may simplify further the acceptance rate. More specifically, we have that(3.63)p(ht+1∣(h)/t+1,r,ϕ)∝p(ht+2∣ht+1,ϕ)p(ht+1∣ht,ϕ)p(rt+1∣ht+2,ht+1,ϕ)p(rt∣ht+1,ht,ϕ).Now, since the following should hold:(3.64)ht+1≥ω+βht
and similarly
(3.65)ht+1≤β-1(ht+2-ω),
we have the support of the conditional distribution of ht+1 given that ht is bounded from below by ω+βht, and the same applies to the distribution of ht+2 given ht+1 (lower limit corresponds to dt=0 and the upper limit to dt+1=0). This means that the range of values of ht+1 compatible with ht and ht+2 in the GQARCH case is bounded from above and below; that is:(3.66)ht+1∈[ω+βht,β-1(ht+2-ω)].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(3.67)g(ht+1∣(h)/t+1,r,ϕ)=p(ht+1∣ht,ϕ)
appropriately truncated from above (since the truncation from below will automatically be satisfied). But the above proposal density ignores the information contained in rt+1, and so according to Fiorentini et al. [10] we can achieve a substantially higher acceptance rate if we propose from
(3.68)g(ht+1∣(h)/t+1,r,ϕ)=p(ht+1∣rt,ht,ϕ).A numerically efficient way to simulateht+1 from p(ht+1∣rt,ht,ϕ) is to sample an underlying Gaussian random variable doubly truncated by using an inverse transform method. More specifically, we may draw(3.69)ɛt∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
doubly truncated so that it remains within the following bounds:(3.70)ɛtnew∈[γ-lt,γ+lt],
where
(3.71)lt=ht+2-ω-βω-β2htβα,
using an inverse transform method and then compute(3.72)ht+1new=ω+α(ɛtnew-γ)2+βht,
which in turn implies a real value for dt+1new=(ht+2-ω-βht+1new)/α and so guarantees that ht+1new lies within the acceptance bounds.The inverse transform method to draw the doubly truncated Gaussian random variable first draws a uniform random number(3.73)u~U(0,1)
and then computes the following:(3.74)u¯=(1-u)Φ(γ-lt-((1-φ2)(rt-δht)/(φu2ht+1-φ2))φu2ht2/(φu2ht+1-φ2))+uΦ(γ+lt-((1-φ2)(rt-δht)/(φu2ht+1-φ2))φu2ht2/(φu2ht+1-φ2)).A draw is then given by(3.75)ɛtnew=Φ-1(u¯).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(3.76)ɛtnew∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
and accept the draw if γ-lt≤ɛtnew≤γ+lt, 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ɛtnew will be given according to the definition of a truncated normal distribution: (3.77)p(ɛtnew∣|ɛtnew-γ|≤lt,rt,ht,ϕ)=1vt/rt,htφ(ɛtnew-ɛt/rr,htvt/rt,ht)×[Φ(γ+lt-ɛt/rr,htvt/rt,ht)-Φ(γ-lt-ɛt/rr,htvt/rt,ht)]-1,
where Φ(·) is the cdf of the standard normal.By using the change of variable formula, we have that the density ofht+1new will be(3.78)p(ht+1new∣ht+1new∈[ω+βht,β-1(ht+2-ω)],rt,ht,ϕ)=ctnew|2αdtnew|×[Φ(γ+lt-ɛt/rr,htvt/rt,ht)-Φ(γ-lt-ɛt/rr,htvt/rt,ht)]-1.Using Bayes theorem we have that the acceptance probability will be(3.79)min⁡(1,p(ht+2∣ht+1new,rt+1,ϕ)p(rt+1∣ht+1new,ϕ)p(rt+1∣ht+1n,ϕ)p(ht+2∣ht+1n,rt+1,ϕ)).Since the degree of truncation is the same for old  and new, the acceptance probability will be(3.80)min⁡(1,p(rt+1∣ht+1new)p(rt+1∣ht+1n)ct+1newct+1ndt+1ndt+1new),
where p(rt+1∣ht+1) is a mixture of two univariate normal densities so(3.81)rt+1∣ht+1~N(δht+1,(φu21-φ2ht+1+1)ht+1).Hence,(3.82)p(rt+1∣ht+1n)=12π((φu2/(1-φ2))ht+1n+1)ht+1nexp⁡(-(rt+1-δht+1n)22((φu2/(1-φ2))ht+1n+1)ht+1n),
and the acceptance probability becomes(3.83)min⁡[1,(ht+1nht+1new)3/2ht+2n-ω-βht+1nht+2n-ω-βht+1newκ(ht+1new)κ(ht+1n)],
where
(3.84)κ(ht+1i)=exp⁡[-(rt+1-δht+1i)22(φu21-φ2ht+1i+1)ht+1i-φu21-φ2ht+1i+12φu21-φ2(ht+1i)2((γ-rt+1-δht+1iφu21-φ2ht+1i+1)2+ht+2r-ω-βht+1iα)]×[1+exp⁡[2φu21-φ2ht+1i+1φu21-φ2(ht+1i)2(γ-rt+1-δht+1iφu21-φ2ht+1i+1)ht+2n-ω-βht+1iα]exp⁡[φu21-φ2ht+1i+1φu21-φ2(ht+1i)2(γ-rt+1-δht+1iφu21-φ2ht+1i+1)ht+2n-ω-βht+1iα]].Overall the MCMC of{ht+1}∣r,ϕ includes the following steps.(1)
Specify an initial value{h(0)}.(2)
Repeat forn=1,…,M.
(a)
Repeat fort=0,…,T-1.
(i)
Use an inverse transform method to simulate(3.85)ɛtnew∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
doubly truncated.(ii)
Calculate(3.86)ht+1new=ω+α(ɛtc-γ)2+βht.
Steps (2)(a)(i) and (2)(a)(ii) are equivalent to draw
(3.87)(ht+1)new~p(ht+1new∣rt,ht,ϕ)
appropriately truncated.(iii)
Calculate(3.88)αr=min⁡[1,p(rt+1∣ht+1new)p(rt+1∣ht+1n)ct+1newct+1ndt+1ndt+1new].(iv)
Set(3.89)(ht+1)n+1={(ht+1)newifUnif(0,1)≤αr(ht+1)notherwiseRemark 3.7.
Every time we changeht+1, we calculate only one normal density since the transformation is Markovian, and, since t=0,…,T-1, we need O(T) calculations.Notice that if we retainht+1new, then ɛtnew is retained and we will not need to simulate st at a later stage. In fact we only need to simulate st at t=T since we need to know ɛT. The final step involves computing(3.90)δt+1(n)=rt+1-ɛt+1(n)ht+1(n),t=0,…,T-1,n=1,…,M.Using all the above simulated values, we may now take average of simulations and compute the quantities needed for the SEM algortihm. As for the Bayesian inference, having completed Step (2) we may now proceed to the Gibbs sampling and M-H steps to obtain draws from the required posterior density. Thus, the first-order Markov transformation of the model made feasible an MCMC algorithm which allows the calculation of a classical estimator via the simulated EM algorithm and a simulation-based Bayesian inference inO(T) computational operations.


### 3.4. A Comparison of the Simulators

In order to compare the performance of the inefficient and the efficient MCMC sampler introduced in the previous subsection, we have generated realizations of sizeT=240 for the simple GQARCH(1,1)-M-AR(1) model with parameters δ=0.1, φ=0.85, α=0.084, β=0.688, γ=0.314 (which are centered around typical values that we tend to see in the empirical literature). We first examine the increase in the variance of the sample mean of ɛt across 500,000 simulations due to the autocorrelation in the drawings relative to an ideal but infeasible independent sampler.We do so by recording the inefficient ratios for the observationst=80 and t=160 using standard spectral density techniques. In addition, we record the mean acceptance probabilities over all observations and the average CPU time needed to simulate one complete drawing. The behavior of the two simulators is summarized in Table 1 and is very much as one would expect. The computationally inefficient sampler shows high serial correlation for both t=80 and t=160 and a low acceptance rate for each individual t. Moreover, it is extremely time consuming to compute even though our sample size is fairly small. In fact, when we increase T from 240 to 2400, and 24000 the average CPU time increases by a factor of 100 and 10000, respectively, as opposed to 10 and 100 for the other one (the efficient), which makes it impossible to implement in most cases of practical interest. On the other hand, the single-move efficient sampler produces results much faster, with a reasonably high acceptance rate but more autocorrelation in the drawings for t=160.Table 1
Comparison between the efficient and inefficient MCMC simulator.
MAPIRTime (CPU effort)t=80t=160T=240T=2400T=24000Inefficient0.1633845.174847.36020.0312504.046875382.05532Single-move0.605858.74229110.2450.0010000.0156250.1083750Note: MAP denotes mean acceptance probability over the whole sample, while IR refers to inefficiency ratio of the MCMC drawings at observations 80 and 160. Time refers to the total CPU time (in seconds) taken to simulate one complete drawing.


## 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ɛt and δt. On the other hand, the conditional log-likelihood function of our model assuming that δt were observed would be the following:
(3.1)l(r,δ∣ϕ,F0)=ln⁡p(r∣δ,ϕ,F0)+ln⁡p(δ∣ϕ,F0)=-Tln⁡2π-12∑t=1Tln⁡ht-12∑t=1T(ɛt)2ht-T2ln⁡(φu2)-12∑t=1T(δt-δ(1-φ)-φδt-1)2φu2,
where r=(r1,…,rT)′,δ=(δ1,…,δT)′, and h=(h1,…,hT)′.However, theδt’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 ϕ(n+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, E(ℓ(·)∣r,ϕ(n),ℱ0) with respect to ϕ keeping ϕ(n) fixed.TheE step, thus, requires the expectation of the complete log-likelihood. For our model, this is given by: (3.2)E(l(⋅)∣r,ϕ(n),F0)=-Tln⁡2π-T2ln⁡(φu2)-12∑t=1TE(ln⁡ht∣r,ϕ(n),F0)-12∑t=1TE((ɛt)2ht∣r,ϕ(n),F0)-12∑t=1TE((δt-δ(1-φ)-φδt-1)2φu2∣r,ϕ(n),F0).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:(3.3)SEMl=-Tln⁡2π-T2ln⁡(φu2)-121M∑i=1M∑t=1Tln⁡ht(i)-121M∑i=1M∑t=1T(ɛt(i))2ht(i)-T2(1-φ)2δ2φu2-121φu21M∑i=1M∑t=1T(δt(i))2+(1-φ)δφu21M∑i=1M∑t=1Tδt(i)+φφu21M∑i=1M∑t=1Tδt(i)δt-1(i)+(1-φ)φδφu21M∑i=1M∑t=1Tδt-1(i)-φ22φu21M∑i=1M∑t=1T(δt-1(i))2.
Consequently, we need to obtain the following quantities: (1/M)∑i=1M∑t=1Tlnht(i),  (1/M)∑i=1M∑t=1T((ɛt(i))2/ht(i)), (1/M)∑i=1M∑t=1Tδt(i), (1/M)∑i=1M∑t=1Tδt-1(i),  (1/M)∑i=1M∑t=1Tδt(i)δt-1(i) and (1/M)∑i=1M∑t=1T(δt(i))2, (1/M)∑i=1M∑t=1Tδt-12(i), where M is the number of simulations.Thus, to classically estimate our model by using an SEM algorithm, the basic problem is to sample fromh∣ϕ,r,ℱ0 where ϕ is the vector of the unknown parameters and also sample from δ∣ϕ,r,ℱ0.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, ɛt, to the price of risk, δt, and vice versa. One possible solution is to fix ω such that E(ht) is 1 or to set φu to a specific value. In fact in an earlier version of the paper, we fixed φu 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δ∣ϕ,r,ℱ0.


## 3.2. Simulation-Based Bayesian Inference

In our problem, the key issue is that the likelihood function of the samplep(r∣ϕ,ℱ0) is intractable which precludes the direct analysis of the posterior density p(ϕ∣r,ℱ0). This problem may be overcome by focusing instead on the posterior density of the model using Bayes’ rule:(3.4)p(ϕ,δ∣r)∝p(ϕ,δ)p(r∣ϕ,δ)∝p(δ∣ϕ)p(ϕ)p(r∣ϕ,δ),
where(3.5)ϕ=(δ,φ,φu,α,β,γ,ω)′.Now,(3.6)p(δ∣ϕ)=∏t=1Tp(δtδ,φ,φu2)=∏t=1T12πφu2exp⁡(-(δt-δ(1-φ)-φδt-1)22φu2).On the other hand,(3.7)p(r∣ϕ,δ)=∏t=1Tp(rt{rt-1},δ,ϕ)=∏t=1T12πhtexp⁡(-ɛt22ht)
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 {δt}(0),ϕ(0), we draw {δt}(1) from p({δt}(1)∣r,ϕ(0)) and then ϕ(1) from p(ϕ(1)∣{δt}(1),r). Iterating these steps, we finally get ({δt}(i),ϕ(i))i=1M, and under mild conditions it is shown that the distribution of the sequence converges to the joint posterior distribution p(ϕ,δ∣r).The above simulation procedure may be carried out by first dividing the parameters into two blocks:(3.8)ϕ1=(δ,φ,φu2),ϕ2=(α,β,γ,ω),
Then the algorithm is described as follows.(1)
Initializeϕ.(2)
Draw fromp(δt∣δ≠t,r,ϕ).(3)
Draw fromp(ϕ∣δ,r) in the following blocks:
(i)
draw fromp(ϕ1∣δ,r) using the Gibbs sampling. This is updated in one block;(ii)
draw fromp(ϕ2∣r) 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(δ,φu2,φ) is given by:(3.9)p(δ,φu2,φ)=p(δ∣φu2)p(φu2)p(φ),
which means that δ,φu2 is a priory independent of φ.Also the following holds for the prior distributions of the parameter subvectorϕ1:
(3.10)p(δ∣φu2)~N(δpr,φu2σδpr2),p(φu2)~IG(v02,d02),p(φ)~N(φ0,σφ02)Iφ,
where Iφ ensures that φ lies outside the unit circle, IG is the inverted gamma distribution, and the hyperparameters v0,d0,δpr,σδpr2,φ0,σφ02 have to be defined.Now, the joint posterior is proportional to(3.11)p(δ,φ,φu2∣r,δ)∝∏t=1T12πφu2exp⁡{-(δt-(1-φ)δ-φδt-1)22φu2}×N(δpr,φu2σδpr2)×IG(v02,d02)×N(φ0,σφ02)Iφ.
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
(3.12)δt*=δt-φδt-1,δt*∣Ft-1~N((1-φ)δ,φu2),
or, otherwise,
(3.13)δt*=(1-φ)δ+vt,vt~N(0,φu2).
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
(3.14)δ∣δ,ϕ,φu2~N(δ̃,φu2σ̃δ2),
where
(3.15)δ̃=σ̃δ2(δprσδpr2+(1-φ)∑t=1Tδt*),σ̃δ2=(1σδpr2+(1-φ)2)-1.
Hence, the generation ofδ is completed, and we may turn on the generation of the other parameters.Generation ofφu2
For the generation ofφu2 and using [27] notation, we have the following.
Proposition 3.2.
The proposal distribution ofφu2 is
(3.16)φu2∣δ,ϕ,δ~IG(T-v02,d0+Q+d2),
where
(3.17)Q=(δ-δpr)2σδpr-2,d=∑t=2T[δt*-δ(1-φ)]2.
Finally, we turn on the generation ofφ.Generation ofφ
For the generation ofφ, we follow again Chib [27] and write
(3.18)δt=(1-φ)δ-φδt-1+vt,vt~N(0,φu2).
We may now state the following proposition (see Chib [27] for a proof).
Proposition 3.3.
The proposal distribution ofφ is
(3.19)φ2∣δ,δ,φu2~N(φ̃,σ̃φ2),
where
(3.20)φ̃=σ̃φ2(σφ0-2φ0+φu-2∑t=1T(δt-1-δ)(δt-δ)),σ̃φ-2=σφ0-2+φu-2∑t=1T(δt-1-δ)2.
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:(3.21)p(α)~N(μα,Σα)Iα,p(β)~N(μβ,σβ2)Iβ
and similarly
(3.22)p(γ)~N(μγ,σγ2),
where α=(ω,α)′,Iα,Iβ 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:(3.23)p(α,β,γ∣r)∝∏t=1T12πhtexp⁡(-ɛt22ht)×N(μα,Σα)Iα×N(μβ,σβ2)Iβ×N(μγ,σγ2).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]:
(3.24)ɛt2=ω+α(ɛt-1-γ)2+βɛt-12+wt-βwt-1,
where wt=ɛt2-ht with wt~N(0,2ht2).Then the corresponding approximated likelihood is written as(3.25)p(ɛ2∣δ,r,ϕ2)=∏t=1T12htπexp⁡[-(ɛt2-ω-α(ɛt-1-γ)2-βɛt-12+βwt-1)24ht2],
and the generation of α,β,γ is based on the above likelihood where we update {ht} 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 wt in (3.32), below, can be written as a linear function of α:
(3.26)wt=ɛ¯t2-ζtα,
where ζt=[ιt̃,ɛ̂t2] with
(3.27)ɛ¯t2=ɛ̃t2-βɛ̃t-12,ɛ̃t2=ɛt2+βɛ̃t-12,ɛ̂t2=(ɛt-1-γ)2+βɛ̂t-12,ιt̃=1+βι̃t-1.
Now, let the two following vectors be(3.28)Yα=[ɛ¯12,…,ɛ¯T2]′,Xα=[ζ1′,…,ζT′]′.
Then the likelihood function of the approximated model is rewritten as(3.29)p(ɛ2∣r,δ,ϕ2)=∏t=1T12π(2ht2)exp⁡[-(ɛ¯t2-ζtα)22(2ht2)].
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
(3.30)α∣Y,X,Σ,ϕ2-α~N(μ̂α,Σ̂α)Iα,
where μ̂α=Σ̂α(Xα′Λ-1Yα+Σα-1μα), Σ̂α=(Xα′Λ-1Xα+Σα-1)-1, and Λ=diag(2h12,…,2hT2). Iα imposes the restriction that α>0 and α+β<1.
Hence a candidateα̃ is sampled from this proposal density and accepted with probability:
(3.31)min⁡{p(α̃,β,γ∣δ,r)q(α*∣α̃,β,γ,δ,r)p(α*,β,γ∣δ,r)q(α̃∣α*,β,γ,δ,r),1},
where α* is the previous draw.
Similar procedure is used for the generation ofβ and γ.Generation ofβ
Following Nakatsuma [29], we linearize wt by the first-order Taylor expansion
(3.32)wt(β)≈wt(β*)+ξt(β*)(β-β*),
where ξt is the first-order derivative of wt(β) evaluated at β* the previous draw of the M-H sampler.
Define as(3.33)rt=wt(β*)+gt(β*)β*,
where gt(β*)=-ξt(β*) which is computed by the recursion:
(3.34)gt=ɛt-12-wt-1+β*gt-1,ξt=0 for t≤0 [30]. Then,
(3.35)wt(β)≈rt-gt(β)β.
Let the following two vectors be(3.36)Yβ=[r1,…,rT]′,Xβ=[g1,…,gT]′.
Then, the likelihood function of the approximated model is rewritten as(3.37)p(ɛ2∣δ,r,ϕ2)=∏t=1T12π(2ht2)exp⁡[-{wt(β*)+ξt(β*)(β-β*)}22(2ht2)].
We have the following proposal distribution forβ (for a proof see Nakatsuma [29] or Ardia [30]).Proposition 3.5.
The proposal distribution forβ is
(3.38)β∣Y,X,σβ2,ϕ2-β~N(μ̂β,σ̂β2)Iβ,
where μ̂β=σ̂β2(Xβ′Λ-1Yβ+(μβ/σβ2)), σ̂β2=(Xβ′Λ-1Xβ+(1/σβ2))-1, and Λ=diag(2h14,…,2hT4). Iβ imposes the restriction that β>0 and α+β<1.
Hence, a candidateβ̃ is sampled from this proposal density and accepted with probability:
(3.39)min⁡{p(β̃,α,γ∣δ,r)q(β*∣β̃,α,γ,δ,r)p(β*,α,γ∣δ,r)q(β̃∣β*,α,γ,δ,r),1}Finally, we explain the generation ofγ.Generation ofγ
As withβ, we linearize wt by a first-order Taylor, expansion at a point γ* the previous draw in the M-H sampler. In this case,
(3.40)rt=wt(γ*)-gt(γ*)γ*,
where gt(γ*)=-ξt(γ*) which is computed by the recursion:
(3.41)gt=-2α(ɛt-1-γ*)+βgt-1,
and gt=0 for t≤0.
Then(3.42)wt(γ)≈rt-gtγ.
Let again(3.43)Yγ=[r1,…,rT]′,Xγ=[g1,…,gT]′
and the likelihood function of the approximated model is rewritten as
(3.44)p(ɛ2∣δ,r,ϕ2)=∏t=1T12π(2ht2)exp⁡[-{wt(γ*)-gt(γ*)γ*}22(2ht2)].
Thus, we have the following proposal distribution forγ (for proof see Nakatsuma [29] and Ardia [30]).Proposition 3.6.
The proposal distribution ofγ is
(3.45)γ∣Y,X,σγ2,ϕ2-γ~N(μ̂γ,σ̂γ2),
where μ̂γ=σ̂γ2(Xγ′Λ-1Yγ+(μγ/σγ2)), σ̂γ2=(Xγ′Λ-1Xγ+(1/σγ2))-1, and Λ=diag(2h14,…,2hT4). A candidate γ̃ is sampled from this proposal density and accepted with probability:
(3.46)min⁡{p(γ̃,α,β∣δ,r)q(γ*∣γ̃,α,β,δ,r)p(γ*,α,β∣δ,r)q(γ̃∣γ*,α,β,δ,r),1}.The algorithm described above is a special case of a MCMC algorithm, which converges as it iterates, to draws from the required densityp(ϕ,δ∣r). 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δ∣ϕ,r,ℱ0.


### 3.2.3. MCMC Simulation ofɛ∣ϕ,r,ℱ0

For a given set of parameter values and initial conditions, it is generally simpler to simulate{ɛt} for t=1,…,T and then compute {δt}t=1T than to simulate {δt}t=1T directly. For that matter, we concentrate on simulators of ɛt given r and ϕ. We set the mean and the variance of ɛ0 equal to their unconditional values, and, given that ht is a sufficient statistic for ℱt-1 and the unconditional variance is a deterministic function of ϕ, ℱ0 can be eliminated from the information set without any information loss.Now sampling fromp(ɛ∣r,ϕ)∝p(r∣ɛ,ϕ)p(ɛ∣ϕ) is feasible by using an M-H algorithm where we update each time only one ɛt leaving all the other unchanged [20]. In particular, let us write the nth  iteration of a Markov chain as ɛtn. Then we generate a potential new value of the Markov chain ɛtnew by proposing from some candidate density g(ɛt∣ɛ∖tn,r,ϕ) where ɛ∖tn={ɛ1n+1,…,ɛt-1n+1,ɛt+1n,…,ɛTn} which we accept with probability(3.47)min⁡[1,p(ɛtnew∣ɛ∖tn,r,ϕ)g(ɛtnew∣ɛ∖tn,r,ϕ)p(ɛtn∣ɛ∖tn,r,ϕ)g(ɛtn∣ɛ∖tn,r,ϕ)].If it is accepted then, we setɛtn+1=ɛtnew and otherwise we keep ɛtn+1=ɛtn. 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:
(3.48)p(ɛtnew∣ɛ∖tn,r,ϕ)p(ɛtn∣ɛ∖tn,r,ϕ)=p(rt∣ɛtnew,htnew,t,ϕ)p(ɛtnew∣htnew,t,ϕ)p(rt∣htn,t,ϕ)p(rt∣htnew,t,ϕ)p(rt∣ɛtn,htn,t,ϕ)p(ɛtnew∣htn,t,ϕ)*∏s=t+1Tp(rs∣ɛsr,hsnew,t,ϕ)p(ɛsr∣hsnew,t,ϕ)p(rs∣hsn,t,ϕ)p(rs∣hsnew,t,ϕ)p(rs∣ɛsn,hsn,t,ϕ)p(ɛsn∣hsn,t,ϕ)=p(rt∣ɛtnew,htnew,t,ϕ)p(ɛtnew∣htnew,t,ϕ)p(rt∣ɛtn,htn,t,ϕ)p(ɛtnew∣htn,t,ϕ)*∏s=t+1Tp(rs∣ɛsn,hsnew,t,ϕ)p(ɛsn∣hsnew,t,ϕ)p(rs∣ɛsn,hsn,t,ϕ)p(ɛsn∣hsn,t,ϕ),
where for s=t+1,…,T,(3.49)hsnew,t=V(ɛs∣ɛs-1n,ɛs-2n,…,ɛt+1n,ɛtnew,ɛt-1n+1,…,ɛ1n+1),hsn,t=V(ɛs∣ɛs-1n,ɛs-2n,…,ɛt+1n,ɛtn,ɛt-1n+1,…,ɛ1n+1)
while(3.50)htnew,t=htn,t.Nevertheless, each time we revise oneɛt, we have also to revise T-t conditional variances because of the recursive nature of the GARCH model which makes hsnew,t depend upon ɛtnew for s=t+1,…,T. And since t=1,…,T, it is obvious that we need to calculate T2 normal densities, and so this algorithm is O(T2). 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 O(T) calculations. The method is described in the following subsection.


## 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(δ,φu2,φ) is given by:(3.9)p(δ,φu2,φ)=p(δ∣φu2)p(φu2)p(φ),
which means that δ,φu2 is a priory independent of φ.Also the following holds for the prior distributions of the parameter subvectorϕ1:
(3.10)p(δ∣φu2)~N(δpr,φu2σδpr2),p(φu2)~IG(v02,d02),p(φ)~N(φ0,σφ02)Iφ,
where Iφ ensures that φ lies outside the unit circle, IG is the inverted gamma distribution, and the hyperparameters v0,d0,δpr,σδpr2,φ0,σφ02 have to be defined.Now, the joint posterior is proportional to(3.11)p(δ,φ,φu2∣r,δ)∝∏t=1T12πφu2exp⁡{-(δt-(1-φ)δ-φδt-1)22φu2}×N(δpr,φu2σδpr2)×IG(v02,d02)×N(φ0,σφ02)Iφ.
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
(3.12)δt*=δt-φδt-1,δt*∣Ft-1~N((1-φ)δ,φu2),
or, otherwise,
(3.13)δt*=(1-φ)δ+vt,vt~N(0,φu2).
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
(3.14)δ∣δ,ϕ,φu2~N(δ̃,φu2σ̃δ2),
where
(3.15)δ̃=σ̃δ2(δprσδpr2+(1-φ)∑t=1Tδt*),σ̃δ2=(1σδpr2+(1-φ)2)-1.
Hence, the generation ofδ is completed, and we may turn on the generation of the other parameters.Generation ofφu2
For the generation ofφu2 and using [27] notation, we have the following.
Proposition 3.2.
The proposal distribution ofφu2 is
(3.16)φu2∣δ,ϕ,δ~IG(T-v02,d0+Q+d2),
where
(3.17)Q=(δ-δpr)2σδpr-2,d=∑t=2T[δt*-δ(1-φ)]2.
Finally, we turn on the generation ofφ.Generation ofφ
For the generation ofφ, we follow again Chib [27] and write
(3.18)δt=(1-φ)δ-φδt-1+vt,vt~N(0,φu2).
We may now state the following proposition (see Chib [27] for a proof).
Proposition 3.3.
The proposal distribution ofφ is
(3.19)φ2∣δ,δ,φu2~N(φ̃,σ̃φ2),
where
(3.20)φ̃=σ̃φ2(σφ0-2φ0+φu-2∑t=1T(δt-1-δ)(δt-δ)),σ̃φ-2=σφ0-2+φu-2∑t=1T(δt-1-δ)2.
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:(3.21)p(α)~N(μα,Σα)Iα,p(β)~N(μβ,σβ2)Iβ
and similarly
(3.22)p(γ)~N(μγ,σγ2),
where α=(ω,α)′,Iα,Iβ 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:(3.23)p(α,β,γ∣r)∝∏t=1T12πhtexp⁡(-ɛt22ht)×N(μα,Σα)Iα×N(μβ,σβ2)Iβ×N(μγ,σγ2).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]:
(3.24)ɛt2=ω+α(ɛt-1-γ)2+βɛt-12+wt-βwt-1,
where wt=ɛt2-ht with wt~N(0,2ht2).Then the corresponding approximated likelihood is written as(3.25)p(ɛ2∣δ,r,ϕ2)=∏t=1T12htπexp⁡[-(ɛt2-ω-α(ɛt-1-γ)2-βɛt-12+βwt-1)24ht2],
and the generation of α,β,γ is based on the above likelihood where we update {ht} 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 wt in (3.32), below, can be written as a linear function of α:
(3.26)wt=ɛ¯t2-ζtα,
where ζt=[ιt̃,ɛ̂t2] with
(3.27)ɛ¯t2=ɛ̃t2-βɛ̃t-12,ɛ̃t2=ɛt2+βɛ̃t-12,ɛ̂t2=(ɛt-1-γ)2+βɛ̂t-12,ιt̃=1+βι̃t-1.
Now, let the two following vectors be(3.28)Yα=[ɛ¯12,…,ɛ¯T2]′,Xα=[ζ1′,…,ζT′]′.
Then the likelihood function of the approximated model is rewritten as(3.29)p(ɛ2∣r,δ,ϕ2)=∏t=1T12π(2ht2)exp⁡[-(ɛ¯t2-ζtα)22(2ht2)].
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
(3.30)α∣Y,X,Σ,ϕ2-α~N(μ̂α,Σ̂α)Iα,
where μ̂α=Σ̂α(Xα′Λ-1Yα+Σα-1μα), Σ̂α=(Xα′Λ-1Xα+Σα-1)-1, and Λ=diag(2h12,…,2hT2). Iα imposes the restriction that α>0 and α+β<1.
Hence a candidateα̃ is sampled from this proposal density and accepted with probability:
(3.31)min⁡{p(α̃,β,γ∣δ,r)q(α*∣α̃,β,γ,δ,r)p(α*,β,γ∣δ,r)q(α̃∣α*,β,γ,δ,r),1},
where α* is the previous draw.
Similar procedure is used for the generation ofβ and γ.Generation ofβ
Following Nakatsuma [29], we linearize wt by the first-order Taylor expansion
(3.32)wt(β)≈wt(β*)+ξt(β*)(β-β*),
where ξt is the first-order derivative of wt(β) evaluated at β* the previous draw of the M-H sampler.
Define as(3.33)rt=wt(β*)+gt(β*)β*,
where gt(β*)=-ξt(β*) which is computed by the recursion:
(3.34)gt=ɛt-12-wt-1+β*gt-1,ξt=0 for t≤0 [30]. Then,
(3.35)wt(β)≈rt-gt(β)β.
Let the following two vectors be(3.36)Yβ=[r1,…,rT]′,Xβ=[g1,…,gT]′.
Then, the likelihood function of the approximated model is rewritten as(3.37)p(ɛ2∣δ,r,ϕ2)=∏t=1T12π(2ht2)exp⁡[-{wt(β*)+ξt(β*)(β-β*)}22(2ht2)].
We have the following proposal distribution forβ (for a proof see Nakatsuma [29] or Ardia [30]).Proposition 3.5.
The proposal distribution forβ is
(3.38)β∣Y,X,σβ2,ϕ2-β~N(μ̂β,σ̂β2)Iβ,
where μ̂β=σ̂β2(Xβ′Λ-1Yβ+(μβ/σβ2)), σ̂β2=(Xβ′Λ-1Xβ+(1/σβ2))-1, and Λ=diag(2h14,…,2hT4). Iβ imposes the restriction that β>0 and α+β<1.
Hence, a candidateβ̃ is sampled from this proposal density and accepted with probability:
(3.39)min⁡{p(β̃,α,γ∣δ,r)q(β*∣β̃,α,γ,δ,r)p(β*,α,γ∣δ,r)q(β̃∣β*,α,γ,δ,r),1}Finally, we explain the generation ofγ.Generation ofγ
As withβ, we linearize wt by a first-order Taylor, expansion at a point γ* the previous draw in the M-H sampler. In this case,
(3.40)rt=wt(γ*)-gt(γ*)γ*,
where gt(γ*)=-ξt(γ*) which is computed by the recursion:
(3.41)gt=-2α(ɛt-1-γ*)+βgt-1,
and gt=0 for t≤0.
Then(3.42)wt(γ)≈rt-gtγ.
Let again(3.43)Yγ=[r1,…,rT]′,Xγ=[g1,…,gT]′
and the likelihood function of the approximated model is rewritten as
(3.44)p(ɛ2∣δ,r,ϕ2)=∏t=1T12π(2ht2)exp⁡[-{wt(γ*)-gt(γ*)γ*}22(2ht2)].
Thus, we have the following proposal distribution forγ (for proof see Nakatsuma [29] and Ardia [30]).Proposition 3.6.
The proposal distribution ofγ is
(3.45)γ∣Y,X,σγ2,ϕ2-γ~N(μ̂γ,σ̂γ2),
where μ̂γ=σ̂γ2(Xγ′Λ-1Yγ+(μγ/σγ2)), σ̂γ2=(Xγ′Λ-1Xγ+(1/σγ2))-1, and Λ=diag(2h14,…,2hT4). A candidate γ̃ is sampled from this proposal density and accepted with probability:
(3.46)min⁡{p(γ̃,α,β∣δ,r)q(γ*∣γ̃,α,β,δ,r)p(γ*,α,β∣δ,r)q(γ̃∣γ*,α,β,δ,r),1}.The algorithm described above is a special case of a MCMC algorithm, which converges as it iterates, to draws from the required densityp(ϕ,δ∣r). 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δ∣ϕ,r,ℱ0.


## 3.2.3. MCMC Simulation ofɛ∣ϕ,r,ℱ0

For a given set of parameter values and initial conditions, it is generally simpler to simulate{ɛt} for t=1,…,T and then compute {δt}t=1T than to simulate {δt}t=1T directly. For that matter, we concentrate on simulators of ɛt given r and ϕ. We set the mean and the variance of ɛ0 equal to their unconditional values, and, given that ht is a sufficient statistic for ℱt-1 and the unconditional variance is a deterministic function of ϕ, ℱ0 can be eliminated from the information set without any information loss.Now sampling fromp(ɛ∣r,ϕ)∝p(r∣ɛ,ϕ)p(ɛ∣ϕ) is feasible by using an M-H algorithm where we update each time only one ɛt leaving all the other unchanged [20]. In particular, let us write the nth  iteration of a Markov chain as ɛtn. Then we generate a potential new value of the Markov chain ɛtnew by proposing from some candidate density g(ɛt∣ɛ∖tn,r,ϕ) where ɛ∖tn={ɛ1n+1,…,ɛt-1n+1,ɛt+1n,…,ɛTn} which we accept with probability(3.47)min⁡[1,p(ɛtnew∣ɛ∖tn,r,ϕ)g(ɛtnew∣ɛ∖tn,r,ϕ)p(ɛtn∣ɛ∖tn,r,ϕ)g(ɛtn∣ɛ∖tn,r,ϕ)].If it is accepted then, we setɛtn+1=ɛtnew and otherwise we keep ɛtn+1=ɛtn. 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:
(3.48)p(ɛtnew∣ɛ∖tn,r,ϕ)p(ɛtn∣ɛ∖tn,r,ϕ)=p(rt∣ɛtnew,htnew,t,ϕ)p(ɛtnew∣htnew,t,ϕ)p(rt∣htn,t,ϕ)p(rt∣htnew,t,ϕ)p(rt∣ɛtn,htn,t,ϕ)p(ɛtnew∣htn,t,ϕ)*∏s=t+1Tp(rs∣ɛsr,hsnew,t,ϕ)p(ɛsr∣hsnew,t,ϕ)p(rs∣hsn,t,ϕ)p(rs∣hsnew,t,ϕ)p(rs∣ɛsn,hsn,t,ϕ)p(ɛsn∣hsn,t,ϕ)=p(rt∣ɛtnew,htnew,t,ϕ)p(ɛtnew∣htnew,t,ϕ)p(rt∣ɛtn,htn,t,ϕ)p(ɛtnew∣htn,t,ϕ)*∏s=t+1Tp(rs∣ɛsn,hsnew,t,ϕ)p(ɛsn∣hsnew,t,ϕ)p(rs∣ɛsn,hsn,t,ϕ)p(ɛsn∣hsn,t,ϕ),
where for s=t+1,…,T,(3.49)hsnew,t=V(ɛs∣ɛs-1n,ɛs-2n,…,ɛt+1n,ɛtnew,ɛt-1n+1,…,ɛ1n+1),hsn,t=V(ɛs∣ɛs-1n,ɛs-2n,…,ɛt+1n,ɛtn,ɛt-1n+1,…,ɛ1n+1)
while(3.50)htnew,t=htn,t.Nevertheless, each time we revise oneɛt, we have also to revise T-t conditional variances because of the recursive nature of the GARCH model which makes hsnew,t depend upon ɛtnew for s=t+1,…,T. And since t=1,…,T, it is obvious that we need to calculate T2 normal densities, and so this algorithm is O(T2). 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 O(T) 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 O(T) calculations. Following their transformation, we augment the state vector with the variables ht+1 and then sample the joint Markov process {ht+1,st}∣r,ϕ∈ℱt where
(3.51)st=sign⁡(ɛt-γ),
so that st=±1 with probability one. The mapping is one to one and has no singularities. More specifically, if we know {ht+1} and φ, then we know the value of(3.52)(ɛt-γ)2=ht+1-ω-βhtα∀t≥1.Hence the additional knowledge of the signs of(ɛt-γ) would reveal the entire path of {ɛt} so long as h0 (which equals the unconditional value in our case) is known, and, thus, we may now reveal also the unobserved random variable {δt}∣r,ϕ,{ht+1}.Now we have to sample from(3.53)p({st,ht+1}∣r,ϕ)∝∏t=1Tp(st∣ht+1,ht,ϕ)p(ht+1∣ht,ϕ)p(rt∣st,ht,ht+1,ϕ),
where the second and the third term come from the model, and the first comes from the fact that ɛt∣ℱt-1~N(0,ht) but ɛt∣{ht+1},ℱt-1 takes values(3.54)ɛt=γ±dt,
where
(3.55)dt=ht+1-ω-βhtα.From the above, it is seen that we should first simulate{ht+1}∣r,ϕ since we do not alter the volatility process when we flip from st=-1 to st=1 (implying that the signs do not cause the volatility process), but we do alter ɛt and then simulate {st}∣{ht+1},r,ϕ. The second step is a Gibbs sampling scheme whose acceptance rate is always one and also conditional on {ht+1},r,ϕ the elements of {st} 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{st}∣{ht+1},r, ϕ

First, we see how to sample from{st}∣{ht+1},r,ϕ. To obtain the required conditionally Bernoulli distribution, we establish first some notation. We have the following (see Appendix A):
(3.56)ct=1vt∣rt,ht[φ(γ+dt-ɛt∣rt,htvt∣rt,ht)+φ(γ-dt-ɛt∣rt,htvt∣rt,ht)],ɛt∣rt,ht=E(ɛt∣rt,ht)=(1-φ2)(rt-δht)φu2ht+1-φ2,vt∣rt,ht=Var⁡(ɛt∣rt,ht)=φu2ht2φu2ht+1-φ2.Using the above notation, we see that the probability of drawingst=1 conditional on {ht+1} is equal to the probability of drawing ɛt=γ+dt conditional on ht+1,ht,rt,ϕ, where dt is given by (3.70), which is given by(3.57)p(st=1∣{ht+1},r,ϕ)=p(ɛt=γ+dt∣ht+1,ht,rt,ϕ)=1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht).Similarly for the probability of drawingst=-1. Both of these quantities are easy to compute; for example,(3.58)φ(γ+dt-ɛt/rr,htvt/rt,ht)=12πexp⁡{-12(γ+dt-ɛt/rr,htvt/rt,ht)2}
and so we may simulate {st}∣{ht+1},r,ϕ using a Gibbs sampling scheme. Specifically, since conditional on {ht+1},r,ϕ the elements of {st} 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{st}∣{ht+1},r,ϕ may be described as below.(1)
Specify an initial values(0)=(s1(0),…,sT(0)).(2)
Repeat fork=1,…,M.
(a)
Repeat fort=0,…,T-1.
(i)
Draws(k)=1 with probability,
(3.59)1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht),
and s(k)=-1 with probability,
(3.60)1-1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht).(3)
Return the values{s(1),…,s(M)}.


### 3.3.2. Simulations of{ht+1}∣r, ϕ (Single Move Samplers)

On the other hand, the first step involves simulating from{ht+1}∣r,ϕ. To avoid large dependence in the chain, we use an M-H algorithm where we simulate one ht+1 at a time leaving the others unchanged [20, 31]. So if (ht+1)n is the current value of the nth iteration of a Markov chain, then we draw a candidate value of the Markov chain ht+1new by proposing it from a candidate density (proposal density) g(ht+1∣(h)/t+1n,r,ϕ) where (h)/t+1n={h1n+1,h2n+1,…,htn+1,ht+2n,…,hT+1n}. We set (ht+1)n+1=(ht+1)new with acceptance probability(3.61)min⁡[1,p(ht+1new∣(h)/t+1n,r,ϕ)g(ht+1n∣(h)/t+1n,r,ϕ)p(ht+1n∣(h)/t+1n,r,ϕ)g(ht+1new∣(h)/t+1n,r,ϕ)],
where we have used the fact that(3.62)p(h∣r,ϕ)=p((h)/t∣r,ϕ)p(ht∣(h)/t,r,ϕ).However, we may simplify further the acceptance rate. More specifically, we have that(3.63)p(ht+1∣(h)/t+1,r,ϕ)∝p(ht+2∣ht+1,ϕ)p(ht+1∣ht,ϕ)p(rt+1∣ht+2,ht+1,ϕ)p(rt∣ht+1,ht,ϕ).Now, since the following should hold:(3.64)ht+1≥ω+βht
and similarly
(3.65)ht+1≤β-1(ht+2-ω),
we have the support of the conditional distribution of ht+1 given that ht is bounded from below by ω+βht, and the same applies to the distribution of ht+2 given ht+1 (lower limit corresponds to dt=0 and the upper limit to dt+1=0). This means that the range of values of ht+1 compatible with ht and ht+2 in the GQARCH case is bounded from above and below; that is:(3.66)ht+1∈[ω+βht,β-1(ht+2-ω)].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(3.67)g(ht+1∣(h)/t+1,r,ϕ)=p(ht+1∣ht,ϕ)
appropriately truncated from above (since the truncation from below will automatically be satisfied). But the above proposal density ignores the information contained in rt+1, and so according to Fiorentini et al. [10] we can achieve a substantially higher acceptance rate if we propose from
(3.68)g(ht+1∣(h)/t+1,r,ϕ)=p(ht+1∣rt,ht,ϕ).A numerically efficient way to simulateht+1 from p(ht+1∣rt,ht,ϕ) is to sample an underlying Gaussian random variable doubly truncated by using an inverse transform method. More specifically, we may draw(3.69)ɛt∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
doubly truncated so that it remains within the following bounds:(3.70)ɛtnew∈[γ-lt,γ+lt],
where
(3.71)lt=ht+2-ω-βω-β2htβα,
using an inverse transform method and then compute(3.72)ht+1new=ω+α(ɛtnew-γ)2+βht,
which in turn implies a real value for dt+1new=(ht+2-ω-βht+1new)/α and so guarantees that ht+1new lies within the acceptance bounds.The inverse transform method to draw the doubly truncated Gaussian random variable first draws a uniform random number(3.73)u~U(0,1)
and then computes the following:(3.74)u¯=(1-u)Φ(γ-lt-((1-φ2)(rt-δht)/(φu2ht+1-φ2))φu2ht2/(φu2ht+1-φ2))+uΦ(γ+lt-((1-φ2)(rt-δht)/(φu2ht+1-φ2))φu2ht2/(φu2ht+1-φ2)).A draw is then given by(3.75)ɛtnew=Φ-1(u¯).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(3.76)ɛtnew∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
and accept the draw if γ-lt≤ɛtnew≤γ+lt, 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ɛtnew will be given according to the definition of a truncated normal distribution: (3.77)p(ɛtnew∣|ɛtnew-γ|≤lt,rt,ht,ϕ)=1vt/rt,htφ(ɛtnew-ɛt/rr,htvt/rt,ht)×[Φ(γ+lt-ɛt/rr,htvt/rt,ht)-Φ(γ-lt-ɛt/rr,htvt/rt,ht)]-1,
where Φ(·) is the cdf of the standard normal.By using the change of variable formula, we have that the density ofht+1new will be(3.78)p(ht+1new∣ht+1new∈[ω+βht,β-1(ht+2-ω)],rt,ht,ϕ)=ctnew|2αdtnew|×[Φ(γ+lt-ɛt/rr,htvt/rt,ht)-Φ(γ-lt-ɛt/rr,htvt/rt,ht)]-1.Using Bayes theorem we have that the acceptance probability will be(3.79)min⁡(1,p(ht+2∣ht+1new,rt+1,ϕ)p(rt+1∣ht+1new,ϕ)p(rt+1∣ht+1n,ϕ)p(ht+2∣ht+1n,rt+1,ϕ)).Since the degree of truncation is the same for old  and new, the acceptance probability will be(3.80)min⁡(1,p(rt+1∣ht+1new)p(rt+1∣ht+1n)ct+1newct+1ndt+1ndt+1new),
where p(rt+1∣ht+1) is a mixture of two univariate normal densities so(3.81)rt+1∣ht+1~N(δht+1,(φu21-φ2ht+1+1)ht+1).Hence,(3.82)p(rt+1∣ht+1n)=12π((φu2/(1-φ2))ht+1n+1)ht+1nexp⁡(-(rt+1-δht+1n)22((φu2/(1-φ2))ht+1n+1)ht+1n),
and the acceptance probability becomes(3.83)min⁡[1,(ht+1nht+1new)3/2ht+2n-ω-βht+1nht+2n-ω-βht+1newκ(ht+1new)κ(ht+1n)],
where
(3.84)κ(ht+1i)=exp⁡[-(rt+1-δht+1i)22(φu21-φ2ht+1i+1)ht+1i-φu21-φ2ht+1i+12φu21-φ2(ht+1i)2((γ-rt+1-δht+1iφu21-φ2ht+1i+1)2+ht+2r-ω-βht+1iα)]×[1+exp⁡[2φu21-φ2ht+1i+1φu21-φ2(ht+1i)2(γ-rt+1-δht+1iφu21-φ2ht+1i+1)ht+2n-ω-βht+1iα]exp⁡[φu21-φ2ht+1i+1φu21-φ2(ht+1i)2(γ-rt+1-δht+1iφu21-φ2ht+1i+1)ht+2n-ω-βht+1iα]].Overall the MCMC of{ht+1}∣r,ϕ includes the following steps.(1)
Specify an initial value{h(0)}.(2)
Repeat forn=1,…,M.
(a)
Repeat fort=0,…,T-1.
(i)
Use an inverse transform method to simulate(3.85)ɛtnew∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
doubly truncated.(ii)
Calculate(3.86)ht+1new=ω+α(ɛtc-γ)2+βht.
Steps (2)(a)(i) and (2)(a)(ii) are equivalent to draw
(3.87)(ht+1)new~p(ht+1new∣rt,ht,ϕ)
appropriately truncated.(iii)
Calculate(3.88)αr=min⁡[1,p(rt+1∣ht+1new)p(rt+1∣ht+1n)ct+1newct+1ndt+1ndt+1new].(iv)
Set(3.89)(ht+1)n+1={(ht+1)newifUnif(0,1)≤αr(ht+1)notherwiseRemark 3.7.
Every time we changeht+1, we calculate only one normal density since the transformation is Markovian, and, since t=0,…,T-1, we need O(T) calculations.Notice that if we retainht+1new, then ɛtnew is retained and we will not need to simulate st at a later stage. In fact we only need to simulate st at t=T since we need to know ɛT. The final step involves computing(3.90)δt+1(n)=rt+1-ɛt+1(n)ht+1(n),t=0,…,T-1,n=1,…,M.Using all the above simulated values, we may now take average of simulations and compute the quantities needed for the SEM algortihm. As for the Bayesian inference, having completed Step (2) we may now proceed to the Gibbs sampling and M-H steps to obtain draws from the required posterior density. Thus, the first-order Markov transformation of the model made feasible an MCMC algorithm which allows the calculation of a classical estimator via the simulated EM algorithm and a simulation-based Bayesian inference inO(T) computational operations.


## 3.3.1. Simulations of{st}∣{ht+1},r, ϕ

First, we see how to sample from{st}∣{ht+1},r,ϕ. To obtain the required conditionally Bernoulli distribution, we establish first some notation. We have the following (see Appendix A):
(3.56)ct=1vt∣rt,ht[φ(γ+dt-ɛt∣rt,htvt∣rt,ht)+φ(γ-dt-ɛt∣rt,htvt∣rt,ht)],ɛt∣rt,ht=E(ɛt∣rt,ht)=(1-φ2)(rt-δht)φu2ht+1-φ2,vt∣rt,ht=Var⁡(ɛt∣rt,ht)=φu2ht2φu2ht+1-φ2.Using the above notation, we see that the probability of drawingst=1 conditional on {ht+1} is equal to the probability of drawing ɛt=γ+dt conditional on ht+1,ht,rt,ϕ, where dt is given by (3.70), which is given by(3.57)p(st=1∣{ht+1},r,ϕ)=p(ɛt=γ+dt∣ht+1,ht,rt,ϕ)=1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht).Similarly for the probability of drawingst=-1. Both of these quantities are easy to compute; for example,(3.58)φ(γ+dt-ɛt/rr,htvt/rt,ht)=12πexp⁡{-12(γ+dt-ɛt/rr,htvt/rt,ht)2}
and so we may simulate {st}∣{ht+1},r,ϕ using a Gibbs sampling scheme. Specifically, since conditional on {ht+1},r,ϕ the elements of {st} 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{st}∣{ht+1},r,ϕ may be described as below.(1)
Specify an initial values(0)=(s1(0),…,sT(0)).(2)
Repeat fork=1,…,M.
(a)
Repeat fort=0,…,T-1.
(i)
Draws(k)=1 with probability,
(3.59)1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht),
and s(k)=-1 with probability,
(3.60)1-1ctvt/rt,htφ(γ+dt-ɛt/rr,htvt/rt,ht).(3)
Return the values{s(1),…,s(M)}.


## 3.3.2. Simulations of{ht+1}∣r, ϕ (Single Move Samplers)

On the other hand, the first step involves simulating from{ht+1}∣r,ϕ. To avoid large dependence in the chain, we use an M-H algorithm where we simulate one ht+1 at a time leaving the others unchanged [20, 31]. So if (ht+1)n is the current value of the nth iteration of a Markov chain, then we draw a candidate value of the Markov chain ht+1new by proposing it from a candidate density (proposal density) g(ht+1∣(h)/t+1n,r,ϕ) where (h)/t+1n={h1n+1,h2n+1,…,htn+1,ht+2n,…,hT+1n}. We set (ht+1)n+1=(ht+1)new with acceptance probability(3.61)min⁡[1,p(ht+1new∣(h)/t+1n,r,ϕ)g(ht+1n∣(h)/t+1n,r,ϕ)p(ht+1n∣(h)/t+1n,r,ϕ)g(ht+1new∣(h)/t+1n,r,ϕ)],
where we have used the fact that(3.62)p(h∣r,ϕ)=p((h)/t∣r,ϕ)p(ht∣(h)/t,r,ϕ).However, we may simplify further the acceptance rate. More specifically, we have that(3.63)p(ht+1∣(h)/t+1,r,ϕ)∝p(ht+2∣ht+1,ϕ)p(ht+1∣ht,ϕ)p(rt+1∣ht+2,ht+1,ϕ)p(rt∣ht+1,ht,ϕ).Now, since the following should hold:(3.64)ht+1≥ω+βht
and similarly
(3.65)ht+1≤β-1(ht+2-ω),
we have the support of the conditional distribution of ht+1 given that ht is bounded from below by ω+βht, and the same applies to the distribution of ht+2 given ht+1 (lower limit corresponds to dt=0 and the upper limit to dt+1=0). This means that the range of values of ht+1 compatible with ht and ht+2 in the GQARCH case is bounded from above and below; that is:(3.66)ht+1∈[ω+βht,β-1(ht+2-ω)].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(3.67)g(ht+1∣(h)/t+1,r,ϕ)=p(ht+1∣ht,ϕ)
appropriately truncated from above (since the truncation from below will automatically be satisfied). But the above proposal density ignores the information contained in rt+1, and so according to Fiorentini et al. [10] we can achieve a substantially higher acceptance rate if we propose from
(3.68)g(ht+1∣(h)/t+1,r,ϕ)=p(ht+1∣rt,ht,ϕ).A numerically efficient way to simulateht+1 from p(ht+1∣rt,ht,ϕ) is to sample an underlying Gaussian random variable doubly truncated by using an inverse transform method. More specifically, we may draw(3.69)ɛt∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
doubly truncated so that it remains within the following bounds:(3.70)ɛtnew∈[γ-lt,γ+lt],
where
(3.71)lt=ht+2-ω-βω-β2htβα,
using an inverse transform method and then compute(3.72)ht+1new=ω+α(ɛtnew-γ)2+βht,
which in turn implies a real value for dt+1new=(ht+2-ω-βht+1new)/α and so guarantees that ht+1new lies within the acceptance bounds.The inverse transform method to draw the doubly truncated Gaussian random variable first draws a uniform random number(3.73)u~U(0,1)
and then computes the following:(3.74)u¯=(1-u)Φ(γ-lt-((1-φ2)(rt-δht)/(φu2ht+1-φ2))φu2ht2/(φu2ht+1-φ2))+uΦ(γ+lt-((1-φ2)(rt-δht)/(φu2ht+1-φ2))φu2ht2/(φu2ht+1-φ2)).A draw is then given by(3.75)ɛtnew=Φ-1(u¯).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(3.76)ɛtnew∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
and accept the draw if γ-lt≤ɛtnew≤γ+lt, 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ɛtnew will be given according to the definition of a truncated normal distribution: (3.77)p(ɛtnew∣|ɛtnew-γ|≤lt,rt,ht,ϕ)=1vt/rt,htφ(ɛtnew-ɛt/rr,htvt/rt,ht)×[Φ(γ+lt-ɛt/rr,htvt/rt,ht)-Φ(γ-lt-ɛt/rr,htvt/rt,ht)]-1,
where Φ(·) is the cdf of the standard normal.By using the change of variable formula, we have that the density ofht+1new will be(3.78)p(ht+1new∣ht+1new∈[ω+βht,β-1(ht+2-ω)],rt,ht,ϕ)=ctnew|2αdtnew|×[Φ(γ+lt-ɛt/rr,htvt/rt,ht)-Φ(γ-lt-ɛt/rr,htvt/rt,ht)]-1.Using Bayes theorem we have that the acceptance probability will be(3.79)min⁡(1,p(ht+2∣ht+1new,rt+1,ϕ)p(rt+1∣ht+1new,ϕ)p(rt+1∣ht+1n,ϕ)p(ht+2∣ht+1n,rt+1,ϕ)).Since the degree of truncation is the same for old  and new, the acceptance probability will be(3.80)min⁡(1,p(rt+1∣ht+1new)p(rt+1∣ht+1n)ct+1newct+1ndt+1ndt+1new),
where p(rt+1∣ht+1) is a mixture of two univariate normal densities so(3.81)rt+1∣ht+1~N(δht+1,(φu21-φ2ht+1+1)ht+1).Hence,(3.82)p(rt+1∣ht+1n)=12π((φu2/(1-φ2))ht+1n+1)ht+1nexp⁡(-(rt+1-δht+1n)22((φu2/(1-φ2))ht+1n+1)ht+1n),
and the acceptance probability becomes(3.83)min⁡[1,(ht+1nht+1new)3/2ht+2n-ω-βht+1nht+2n-ω-βht+1newκ(ht+1new)κ(ht+1n)],
where
(3.84)κ(ht+1i)=exp⁡[-(rt+1-δht+1i)22(φu21-φ2ht+1i+1)ht+1i-φu21-φ2ht+1i+12φu21-φ2(ht+1i)2((γ-rt+1-δht+1iφu21-φ2ht+1i+1)2+ht+2r-ω-βht+1iα)]×[1+exp⁡[2φu21-φ2ht+1i+1φu21-φ2(ht+1i)2(γ-rt+1-δht+1iφu21-φ2ht+1i+1)ht+2n-ω-βht+1iα]exp⁡[φu21-φ2ht+1i+1φu21-φ2(ht+1i)2(γ-rt+1-δht+1iφu21-φ2ht+1i+1)ht+2n-ω-βht+1iα]].Overall the MCMC of{ht+1}∣r,ϕ includes the following steps.(1)
Specify an initial value{h(0)}.(2)
Repeat forn=1,…,M.
(a)
Repeat fort=0,…,T-1.
(i)
Use an inverse transform method to simulate(3.85)ɛtnew∣rt,ht,ϕ~N((1-φ2)(rt-δht)φu2ht+1-φ2,φu2ht2φu2ht+1-φ2)
doubly truncated.(ii)
Calculate(3.86)ht+1new=ω+α(ɛtc-γ)2+βht.
Steps (2)(a)(i) and (2)(a)(ii) are equivalent to draw
(3.87)(ht+1)new~p(ht+1new∣rt,ht,ϕ)
appropriately truncated.(iii)
Calculate(3.88)αr=min⁡[1,p(rt+1∣ht+1new)p(rt+1∣ht+1n)ct+1newct+1ndt+1ndt+1new].(iv)
Set(3.89)(ht+1)n+1={(ht+1)newifUnif(0,1)≤αr(ht+1)notherwiseRemark 3.7.
Every time we changeht+1, we calculate only one normal density since the transformation is Markovian, and, since t=0,…,T-1, we need O(T) calculations.Notice that if we retainht+1new, then ɛtnew is retained and we will not need to simulate st at a later stage. In fact we only need to simulate st at t=T since we need to know ɛT. The final step involves computing(3.90)δt+1(n)=rt+1-ɛt+1(n)ht+1(n),t=0,…,T-1,n=1,…,M.Using all the above simulated values, we may now take average of simulations and compute the quantities needed for the SEM algortihm. As for the Bayesian inference, having completed Step (2) we may now proceed to the Gibbs sampling and M-H steps to obtain draws from the required posterior density. Thus, the first-order Markov transformation of the model made feasible an MCMC algorithm which allows the calculation of a classical estimator via the simulated EM algorithm and a simulation-based Bayesian inference inO(T) computational operations.


## 3.4. A Comparison of the Simulators

In order to compare the performance of the inefficient and the efficient MCMC sampler introduced in the previous subsection, we have generated realizations of sizeT=240 for the simple GQARCH(1,1)-M-AR(1) model with parameters δ=0.1, φ=0.85, α=0.084, β=0.688, γ=0.314 (which are centered around typical values that we tend to see in the empirical literature). We first examine the increase in the variance of the sample mean of ɛt across 500,000 simulations due to the autocorrelation in the drawings relative to an ideal but infeasible independent sampler.We do so by recording the inefficient ratios for the observationst=80 and t=160 using standard spectral density techniques. In addition, we record the mean acceptance probabilities over all observations and the average CPU time needed to simulate one complete drawing. The behavior of the two simulators is summarized in Table 1 and is very much as one would expect. The computationally inefficient sampler shows high serial correlation for both t=80 and t=160 and a low acceptance rate for each individual t. Moreover, it is extremely time consuming to compute even though our sample size is fairly small. In fact, when we increase T from 240 to 2400, and 24000 the average CPU time increases by a factor of 100 and 10000, respectively, as opposed to 10 and 100 for the other one (the efficient), which makes it impossible to implement in most cases of practical interest. On the other hand, the single-move efficient sampler produces results much faster, with a reasonably high acceptance rate but more autocorrelation in the drawings for t=160.Table 1
Comparison between the efficient and inefficient MCMC simulator.
MAPIRTime (CPU effort)t=80t=160T=240T=2400T=24000Inefficient0.1633845.174847.36020.0312504.046875382.05532Single-move0.605858.74229110.2450.0010000.0156250.1083750Note: MAP denotes mean acceptance probability over the whole sample, while IR refers to inefficiency ratio of the MCMC drawings at observations 80 and 160. Time refers to the total CPU time (in seconds) taken to simulate one complete drawing.


## 4. Empirical Application: Bayesian Estimation of Weekly Excess Returns from Three Major Stock Markets: Dow-Jones, FTSE, and Nikkei

In this section we apply the procedures described above to weekly excess returns from three major stock markets: Dow-Jones, FTSE, and Nikkei for the period of the last week of 1979:8 to the second to the last week of 2008:5 (1,500 observations). To guarantee0≤β≤1-α≤1 and to ensure that ω>0 we also used some accept-reject method for the Bayesian inference. This means that, when drawing from the posterior (as well as from the prior), we had to ensure that α,β>0,α+β<1 and ω>0.In order to implement our proposed Bayesian approach, we first have to specify the hyperparameters that characterize the prior distributions of the parameters. In this respect, our aim is to employ informative priors that would be in accordance with the “received wisdom.” In particular, for all data sets, we set the prior mean forβ equal to 0.7, and for a, ω, and γ we decided to set their prior means equal to 0.15, 0.4, and 0.0, respectively. We had also to decide on the prior mean of δ. We set its prior mean equal to 0.05, for all markets. These prior means imply an annual excess return of around 4%, which is a typical value for annualized stock excess returns. Finally, we set the prior mean of φ equal to 0.75, of φu2 equal to 0.01, and the hyperparameters ν0 and do equal to 1550 and 3, respectively, for all three datasets, something which is consistent with the “common wisdom” of high autocorrelation of the price of risk. We employ a rather vague prior and set its prior variance equal to 10,000 for all datasets.We ran a chain for 200,000 simulations for the three datasets and decided to use every tenth point, instead of all points, in the sample path to avoid strong serial correlation. The posterior statistics for the Dow-Jones, FTSE, and Nikkey are reported in Table2. Inefficiency factors are calculated using a Parzen window equal to 0.1T (where, recall, T is the number of observations) and indicate that the M-H sampling algorithm has converged and well behaved (This is also justified by the ACFs of the draws. However, they are not presented for space considerations and are available upon request.) With the exception of the constants δ and the φu2's, there is low uncertainty with the estimation of the parameters. The estimated persistence, α+β, for all three markets is close to 0.8 with the highest being the one of FTSE (0.822), indicating that the half life of a shock is around 3.5. The estimated asymmetry parameters are round 0.4 with “t-statistics” higher than 3.2, indicating that the leverage effect is important in all three markets. In a nut shell, all estimated parameters have plausible values, which are in accordance with previous results in the literature.Table 2
Bayesian inference results.
Dow-JonesPMPSDϕ0.5ϕminϕmaxIFδ0.0520.0410.054−0.0480.1033.001φ0.8120.0820.854−0.9920.9992.115φu20.0100.0340.0130.0020.0181.509ω0.4050.0710.4610.0040.8162.367α0.1520.0400.2910.0010.8731.048β0.6510.1680.6290.0030.9842.994γ0.3920.1120.393−0.6810.4185.108FTSEPMPSDϕ0.5ϕminϕmaxIFδ0.0590.0360.059−0.0510.3953.111φ0.8110.0960.839−0.8090.9992.154φu20.0090.0290.0120.0050.0171.995ω0.2050.0870.3980.0030.9951.457α0.1400.0550.1870.0010.9313.458β0.6820.1530.7010.0010.9882.721γ0.3740.1020.381−0.6150.4011.254NikkeiPMPSDϕ0.5ϕminϕmaxIFδ0.0680.0510.068−0.0640.2112.998φ0.8090.0900.837−0.8310.9991.211φu20.0100.0310.0100.0040.0192.001ω0.1950.0790.2280.0040.5012.789α0.1490.0520.1970.0010.8933.974β0.6340.1190.6450.0060.9891.988γ0.4080.1230.409−0.5870.4874.007Note: PM denotes posterior mean, PSD posterior standard deviation,ϕ0.5 posterior median, ϕmin posterior minimum, ϕmax posterior maximum, and IF inefficiency factor.We have also performed a sensitivity analysis to our choice of priors. In particular, we have halved and doubled the dispersion of the prior distributions around their respective means. Figures1, 2, and 3 show the kernel density estimates for all parameters for all datasets for the posterior distributions for the three cases: when the variances are 10,000 (baseline posterior), when the variances are halved (small variance posterior), and when the variances are doubled (large variance posterior). We used a canonical Epanechnikov kernel, and the optimal bandwidth was determined automatically by the data. The results which are reported in Figures 1, 2, and 3 indicate that the choice of priors does not unduly influence our conclusions.Dow-Jones: posterior density estimates and sensitivity analysis.
(a)(b)(c)(d)(e)(f)(g)FTSE: posterior density estimates and sensitivity analysis.
(a)(b)(c)(d)(e)(f)(g)Nikkei: posterior density estimates and sensitivity analysis.
(a)(b)(c)(d)(e)(f)(g)Finally, treating the estimated posterior means as the “true” parameters, we can employ the formulae of Section2.1 and compare the moments implied by the estimates and the sample ones. One fact is immediately obvious. All order autocorrelations of excess returns implied by the estimates are positive but small, with the 1st one being around 0.04, which is in accordance with the (ii) stylized fact (see Section 1). However, for all the three markets, the sample skewness coefficients are negative, ranging from −0.89 (FTSE) to −0.12 (Nikkey), whereas the implied ones are all positive, ranging from 0.036 (FTSE) to 0.042 (Dow-Jones). Nevertheless, the model is matching all the other stylized facts satisfactorily, that is, the estimated parameter values accommodate high coefficient of variation, leptokurtosis as well the volatility clustering and leverage effect.


## 5. Conclusions

In this paper, we derive exact likelihood-based estimators for our time-varying GQARCH(1,1)-M model. Since in general the expression for the likelihood function is unknown, we resort to simulation methods. In this context, we show that MCMC likelihood-based estimation of such a model can in fact be handled by means of feasibleO(T) algorithms. Our samplers involve two main steps. First we augment the state vector to achieve a first-order Markovian process in an analogous manner to the way in which GARCH models are simulated in practice. Then, we discuss how to simulate first the conditional variance and then the sign given these simulated series so that the unobserved in mean process is revealed as a residual term. We also develop simulation-based Bayesian inference procedures by combining within a Gibbs sampler the MCMC simulators. Furthermore, we derive the theoretical properties of this model, as far as moments and dynamic moments are concerned.In order to investigate the practical performance of the proposed procedure, we estimate within a Bayesian context our time-varying GQARCH(1,1)-M-AR(1) model for weekly excess stock returns from the Dow-Jones, Nikkei, and FTSE index. With the exception of the returns' skewness, the suggested likelihood-based estimation method and our model is producing satisfactory results, as far as a comparison between sample and theoretical moments is concerned.Although we have developed the method within the context of an AR(1) price of risk, it applies much more widely. For example, we could assume that the market price of risk is a Bernoulli process or a Markov’s switching process. A Bernoulli’s distributed price of risk would allow a negative third moment by appropriately choosing the two values of the in-mean process. However, this would make all computations much more complicated.  In an earlier version of the paper, we assumed that the market price of risk follows a normal distribution, and we applied both the classical and the Bayesian procedure to three stock markets (where we decided to set the posterior means as initial values for the simulated EM algorithm). The results suggested that the Bayesian and the classical procedures are quite in agreement (see Anyfantaki and Demos [24]).Finally, it is known that (e.g., [32, pages 84 and 85]) the EM algorithm slows down significantly in the neighborhood of the optimum. As a result, after some initial EM iterations, it is tempting to switch to a derivative-based optimization routine, which is more likely to quickly converge to the maximum. EM-type arguments can be used to facilitate this switch by allowing the computation of the score. In particular, it is easy to see that
(5.1)E(∂ln⁡p(δ∣r,ϕ,F0)∂ϕ∣r,ϕ(n),F0)=0,
so it is clear that the score can be obtained as the expected value given r,ϕ,ℱ0 of the sum of the unobservable scores corresponding to lnp(r∣δ,ϕ,ℱ0) and lnp(δ∣ϕ,ℱ0). This could be very useful for the classical estimation procedure, not presented here, as even though our algorithm is an O(T) one, it is still rather slow. We leave these issues for further research.



---

*Source: 718647-2011-09-22.xml*