HITOTSUBASHI

UNIVERSITY

Center for Financial Engineering Education (CFEE)CFEE Discussion Paper Series No. 2013-1

Matrix exponential stochastic volatility with cross leverage

Tsunehiro Ishihara, Yasuhiro Omori and Manabu Asai

Center for Financial Engineering Education (CFEE)
Department of Economics, Hitotsubashi University
2-1 Naka, Kunitachi Tokyo, 186-8601 Japan
http://www.econ.hit-u.ac.jp/~finmodelMatrix exponential stochastic volatility with cross leverage
Tsunehiro Ishihara∗

Yasuhiro Omori † andManabu Asai ‡October 14, 2013

Abstract
A multivariate stochastic volatility model with the dynamic correlation and thecross leverage effect is described and its efficient estimation method using Markov chain
Monte Carlo is proposed. The time-varying covariance matrices are guaranteed to bepositive definite by using a matrix exponential transformation. Of particular interest is our approach for sampling a set of latent matrix logarithm variables from theirconditional posterior distribution, where we construct the proposal density based on
an approximating linear Gaussian state space model. The proposed model and its extended models with fat-tailed error distribution are applied to trivariate returns data
(daily stocks, bonds, and exchange rates) of Japan. Further, a model comparison isconducted including constant correlation multivariate stochastic volatility models withleverage.

Key words: Dynamic correlation, Leverage effect, Matrix exponential, Markov chain
Monte Carlo, Multi-move sampler, Multivariate stochastic volatility∗

Department of Economics, Hitotsubashi University, 2-1, Naka, Kunitachi, Tokyo 186-8601, Japan.
Email: T.Ishihara@r.hit-u.ac.jp†
Faculty of Economics, University of Tokyo. 7-3-1 Hongo, Bunkyo-Ku, Tokyo 113-0033, Japan. Tel:
+81-3-5841-5516. Email: omori@e.u-tokyo.ac.jp‡
Faculty of Economics, Soka University, Hachioji-shi, Tokyo 192-8577, Japan. E-mail: m-asai@soka.ac.jp1

1

Introduction

Over the last several decades, there has been a great deal of interest in modeling volatilities ofmultivariate stock market returns. The examples are multivariate generalized autoregressive
conditional heteroskedasticity (GARCH) models (see the review of Bauwens, Laurent, andRombouts (2006)), multivariate stochastic volatility (SV) models (see the review of Asai,
McAleer, and Yu (2006), Chib, Omori, and Asai (2009) ) and realized covariance models (e.g.
Andersen, Bollerslev, Diebold, and Ebens (2001), Barndorff-Nielsen and Shephard (2004),
Chiriac and Voev (2010), Audrino and Corsi (2010), Voev (2007), Bauer and Vorkink (2010),
Jin and Maheu (2010) and Sheppard (2007)). The realized covariance model uses the highfrequency data to estimate covariance matrices and regard them as observed covariancematrices, while they are latent variables in GARCH and SV models.
Various multivariate volatility models have been proposed in the literature to describe
the dynamic properties of the covariance matrices such as the volatility clustering, thedynamic correlations, and the leverage effects. The DCC models, (Engle (2002)) and BEKK
model (Engle and Kroner (1995)) are such multivariate GARCH models, and autoregressive
Wishart models (Philipov and Glickman (2006) and Gourieroux, Jasiak, and Sufana (2009))
are examples in multivariate SV models. The common difficulty in these models is to keep
the covariance matrices positive definite. To overcome this difficulty, reparameterization
methods are considered in Yu and Meyer (2006), Tsay (2005), and Jungbacker and Koopman(2006). The Choleski decomposition of the covariance matrix is also considered in Lopes,
McCulloch, and Tsay (2012).
However, there have been still few previous works on the multivariate volatility models
with both dynamic correlations and cross leverage effects. Cross leverage refers to thecorrelation between the i-th asset return at time t and the function of j-th asset volatility
at time t+1 (when i = j, we simply call it a leverage effect). Thus, to model these propertiesof covariance matrices, this paper considers the matrix logarithm transformation which is
known useful to model positive definite matrices in a flexible way. Since the seminal workof Chiu, Leonard, and Tsui (1996), the matrix exponential model for the covariance matrix
has been applied to the spatial model to simplify the calculation of log-likelihood functions(LeSage and Pace (2007)), and is extended to the GARCH model (Kawakatsu (2006)), theSV model (Asai, McAleer, and Yu (2006)) and the realized covariation model (Bauer andVorkink (2010) and Sheppard (2007)) for multivariate financial time series.
We consider the general multivariate volatility model using the matrix exponential SV
2model with cross leverage effects and propose an efficient computational algorithm. Thisis a generalization of Ishihara and Omori (2012) who propose the following multivariate
stochastic volatility (MSV) model with cross-asset leverage effect of the form
yt = diag (exp(α1t /2), . . . , exp(αpt /2)) εt ,αt+1 = Φαt + ηt ,¡ 0 0 ¢0
εt , ηt∼ N2p (0, Σ),

(1)
(2)
(3)

where yt = (y1t , . . . , ypt )0 , αt = (α1t , . . . , αpt )0 , Φ = diag(φ1 , . . . , φp ) and Np (µ, Σ) denotes
the p-dimensional normal distribution with mean µ and variance Σ. This is fairly generalin the sense that there is no restriction imposed on the covariance matrix Σ, while, inthe previous literature, various parameter restrictions are imposed (e.g. Asai and McAleer(2006), Asai and McAleer (2009), Chan, Kohn, and Kirby (2006), and Danı́elsson (1998))
to estimate parameters based on the Monte Carlo likelihood. We, further, model the dynamic covariance matrices (dynamic variances and correlations) using a matrix logarithm
transformation. Since it is difficult to implement a maximum likelihood estimation for our
proposed model without imposing restrictions on parameters, we take Bayesian approachand estimate posterior distributions of model parameters using Markov chain Monte Carlo(MCMC) method. The simple sampling algorithm for the latent covariance matrices isknown to be inefficient as discussed in Ishihara and Omori (2012). They showed that thesingle-move sampler which samples one volatility variable given others is highly inefficientand proposed the efficient multi-move sampler (block sampler) which divides the vector of
all latent variables into blocks and samples one block given other blocks based on Omori
and Watanabe (2008). Thus we construct the multi-move sampler for our matrix exponential model and show that it is efficient in comparison with the alternative simple samplingalgorithm.
The rest of the paper is organized as follows. In Section 2, we introduce an matrix
exponential stochastic volatility model with cross leverage effects. Bayesian estimation
method and the associated particle filter for calculating likelihood functions are describedin Section 3. Section 4 shows the efficiency of our proposed algorithm using the simulateddata, and, in Section 5, the empirical studies are given using the trivariate asset returns
data (stock indices, bond indices and foreign exchange rates). We conduct a model selection
among the proposed model, extended models with fat-tailed error distribution and someconstant correlation multivariate SV models. Section 6 concludes the paper.
3

2Matrix exponential stochastic volatility

This section proposes the matrix exponential stochastic volatility (MESV) model with crossleverage effects. The MESV model is based on the matrix exponential transformation as
below. A matrix exponential is widely studied in the context of multidimensional differential equations and Lie algebra. The statistical applications of the matrix exponential
transformation are given, for example, in Chiu, Leonard, and Tsui (1996), and Kawakatsu(2006). For any p × p matrix A, the matrix exponential is defined by the following power
series expansion
exp(A) ≡

∞
X
1
s=0

s!

As ,

where the series converges absolutely if all eigenvalues of A are finite. ( see e.g. Abadir andMagnus (2005) for various properties of the matrix exponential transformation). For any
real symmetric positive definite matrix C, there exists a real symmetric p × p matrix A suchthat C = exp(A), and the matrix A is obtained by the matrix logarithm transformation.
Conversely, for any real symmetric matrix A, C = exp(A) is a symmetric positive definite
matrix (Chiu, Leonard, and Tsui (1996)). If A is a p × p real symmetric matrix, there exists
a p × p orthogonal matrix U and a diagonal matrix Λ such that A = UΛU0 andexp(A) = U

Ã∞
X1
s=0

s!

!
Λs

U0 = U exp(Λ)U0 .

Now let yt = (y1t , . . . , ypt )0 denote the p-dimensional asset return vector at time t, and let
Ht denote the matrix logarithm of the variance-covariance matrix of yt . The MESV modelwith leverage effects is given by
yt = exp(Ht /2)εt ,εt ∼ i.i.d. N (0, Ip ),t = 1, . . . , n.

Ht+1 = M + Φ̃ ¯ (Ht − M) + Et ,εt
Ip Σεη ∼ i.i.d. Np+q (0, Σ), Σ = ,ηtΣηε Σηη

(5)



h1 ∼

Nq (µ, Σ0 ) ,(4)

t = 1, . . . , n − 1,

(6)
(7)

where ηt = vech(Et ), q = p(p + 1)/2, M = {µij }, and Φ̃ = {φij } are p × p symmetricmatrices of parameters, and ¯ denotes the Hadamard product. For the identifiability, we

4set the covariance matrix of εt equal to Ip .
If we let ht = vech(Ht ) = (h11,t , h21,t , . . . , hp1,t , h22,t . . . , hpp,t )0 denote the stacked column vector of the lower triangle elements of the Ht , then the “vech form” of (5) is givenbyht+1 = µ + Φ(ht − µ) + ηt ,(8)with µ = vech(M) = (µ11 , µ21 , . . . µp1 , µ22 , . . . , µpp )0 , Φ = diag(φ) (a diagonal matrix whose
diagonal elements are equal to φ) and φ = vech(Φ̃) = (φ11 , φ21 , . . . φp1 , φ22 , . . . , φpp )0 . The
number of parameters in the MESV model is q(q + 2p + 3)/2. The covariance matrix of theinitial latent variable, Σ0 , is assumed to satisfy a stationary condition such thatvec(Σ0 ) = (Ip2 − Φ ⊗ Φ)−1 vec(Σηη ).
We let Σηη = {ρij,ηη σi,ηη σj,ηη }, and Σεη = {ρij,εη σj,ηη } where σi,ηη is the standard deviationof ηit and ρij,xy is the correlation coefficient between xit and yjt . Further, for convenience, we
use the notation E(i, j) = k based on the relationship ηt = vech(Et ) such that the (i, j)-thelement of Et , Et (i, j), is equal to the k-th element of ηt , ηkt (i.e., E(1, 1) = 1, E(2, 1) = 2,
. . . , E(p, 1) = p, E(2, 2) = p + 1,. . . , E(p, p) = p(p + 1)/2). Thus, Cov(εlt , ηkt ) = ρlk,εη σk,ηη
is equal to Cov(εlt , Et (i, j)) = ρlE(i,j),εη σE(i,j),ηη .
Remark. Due to the nonlinearity of the matrix exponential transformation, the interpretation of the (untransformed) parameters will depend on the dimension of yt . Thus, weconsider estimates of transformed parameters to investigate the properties of interest, suchas volatilities, correlations, principal components, and the news impact curve.

3

Bayesian estimation and associated particle filter

In this section, we describe an efficient Bayesian estimation method and an associatedparticle filter to compute the likelihood for the MESV model. Let θ = (φ, µ, Σ) andh = (h01 , . . . , h0n )0 and Yn = (y1 , . . . , yn ). Then the joint probability density function of Yn

5and h given θ for (4) and (8) is given byf (Yn , h|θ) = f (h1 |θ)n−1Yf (yt , ht+1 |ht , θ)f (yn |hn , θ)

t=1
− 21

∝ |Σ0 |

− n−1
2

|Σ|

"

1
× exp −2− 21|Σεε |

exp

( n
X
t=1

n−1
X

(9)

)
1
0 −1
lt − (h1 − µ) Σ0 (h1 − µ)2#

{ht+1 − µ − Φ(ht − µ)}0 Σ−1ηη {ht+1 − µ − Φ(ht − µ)} ,

t=1

©
ª
where lt = − 12 tr(Ht ) + (yt − µt )0 Σ−1t (yt − µt ) and3.1µt = exp(Ht /2)mt ,(10)

Σt = exp(Ht /2)St exp(Ht /2),

(11)

mt = Σεη Σ−1ηη (ht+1 − µ − Φ(ht − µ))I(t < n),(12)

St = Ip − Σεη Σ−1ηη Σηε I(t < n).

(13)

Prior distributions

For prior distributions of (φ, µ), we assume
φij + 1
∼ B(aij , bij ), i = 1, . . . , p,2
µ ∼ Nq (m∗0 , V0∗ ),j = 1, . . . , i,

(14)
(15)

where B(a, b) denotes a beta distribution with parameters a and b. To define a priordistribution of Σ, we first denoteΣ−1 = 


Σ11

Σ12

Σ21

Σ22

,

where Σ11 , Σ12 and Σ22 are p × p, p × q and q × q matrices. Noting that Σ11 = Ip +
Σ12 Σ22−1 Σ21 , we assume the prior distributions such thatvec(Σ21 )|Σ22 ∼ Npq (vec(Σ22 ∆0 ), Ω0 ⊗ Σ22 ),Σ22 ∼ W(n0 , R0 ),where W(n, R) denotes Wishart distribution with parameters n and R.

6

(16)

3.2

MCMC algorithm

Using Equations (9), (14), (15) and (16), we obtain the joint posterior density function of(θ, h) given byπ(θ, h|Yn )∝ f (Yn , h|θ) ×

p Y
i
Y

(1 + φij )aij −1 (1 − φij )bij −1 × fN (µ|m∗0 V0∗ )i=1 j=1¾½¡¢
1 ¡ −1 22 ¢21
22
22
22 n−q−1
×fN vec(Σ )|vec(Σ ∆0 ), Ω0 ⊗ Σ
× |Σ | 2 exp − tr R0 Σ,(17)2
where fN (·|µ, Σ) denotes a normal density with mean µ and covariance matrix Σ. To
obtain the posterior quantities of the parameters θ and volatility variables {ht }nt=1 fromthe posterior distribution, we implement the MCMC algorithm in six blocks:
1. Initialize h, φ, µ, Σ.
2. Generate h|φ, µ, Σ, Yn .
3. Generate µ|φ, Σ, h, Yn .
4. Generate Σ|φ, µ, h, Yn .
5. Generate φ|µ, Σ, h, Yn .
6. Go to Step 2.
3.2.1

Generation of h

As is often pointed out in the literature, it is important to sample the latent volatility
variables {ht }nt=1 in an efficient way. The simple sample method, which samples one ht
at a time given the other hs ’s and parameters, is known to be inefficient, often producing
highly autocorrelated MCMC samples. This is because the the estimates of autoregressive
parameters φi are often found to be very close to one in empirical studies. Thus, we propose
the sampling method based on a multi-move sampler which samples a set of ht ’s as one
block at a time (see e.g.Shephard and Pitt (1997), Watanabe and Omori (2004), Omori
and Watanabe (2008), Ishihara and Omori (2012)). We first describe a simple algorithm
which we call a single-move sampler as we use it as a benchmark to evaluate the estimation
efficiencies of the multi-move algorithm.

7

Single-move sampler. Let Σh = Σηη −Σηε Σεη and µh,t+1 = µ+Φ(ht −µ)+Σηε exp(−Ht /2)yt .
Then, the conditional posterior density of ht given {hs }s6=t , Φ and Σ is
½
¾
1
0 −1
π(ht |{hs }s6=t , Φ, Σ, Yn ) ∝ exp − (ht − γt ) Γt (ht − γt ) + g(ht ) ,
2
where
1
1
g(ht ) = − tr(Ht ) − yt0 exp(−Ht /2)S−1
t exp(−Ht /2)yt22+(ht+1 − µ − Φ(ht − µ))0 Σ−1h Σηε exp(−Ht /2)yt I(t < n),
and

Γt

γt


−1
−1 −1


 (ΦΣh Φ + Σ0 ) , t = 1,
−1
−1
=
(Σ−1
h + ΦΣh Φ) , 1 < t < n,



Σh ,
t = n,
−1

t = 1,
 Γ1 ΦΣh (h2 − (Ip − Φ)µ),−1−1=Γt (ΦΣh (ht+1 − (Ip − Φ)µ) + Σh µh,t ), 1 < t < n,µh,t ,
t = n.

We generate a candidate h†t from h†t ∼ Nq (γt , Γt ) and accept it with probabilitynomin exp{g(h†t ) − g(ht )}, 1 ,t = 1, . . . , n.

Multi-move sampler. In this algorithm, we first divide h into several blocks, and samples
one block at a time from its conditional posterior distribution given other blocks. Using the
Taylor expansion of the logarithm of the conditional posterior density around the conditional
posterior mode, we derive a candidate distribution as a posterior distribution for some linear
Gaussian state space model to exploit various smoothing and simulation algorithms as in
Omori and Watanabe (2008).
To generate {ht }s+mt=s+1 given other ht ’s, for example, we sample the normalized disturs+m−1
bances {xt }t=sinstead of {ht }s+mt=s+1 since such a sampling method is known to reduce

the MCMC sample autocorrelations wherext = Σ−1/2ηη ηt ,

t = 1, . . . , n − 1,

8

−1/2

x0 = Σ0

η0 ,

1/2

1/2

and Σηη , Σ0

1/2

1/20

denote Choleski decompositions such that Σηη = Σηη Σηηand Σ0 =1/2 1/20
Σ0 Σ0 .

Consider the following linear Gaussian state space model:ŷt = Zt αt + Gt ut ,t = s + 1, . . . , s + m,

(18)

αt+1 = Φαt + Kt ut ,t = s + 1, . . . , s + m − 1,

(19)

ut = (ξt0 , x0t )0 ∼ N2p (0, I) ,
h
i
h
i
1/2
1/2
where Kt = O, Σηη for t = 1, . . . , n−1 and K0 = O, Σ0 , and ŷt , Zt , Gt are computed
as follows:
1. First we set x̂t = xt (t = s, . . . , s + m − 1) where xt is a current sample.
2. Define dt , At and Bt as in Appendix A. Let dˆt , Ât and B̂t denote the vector and1/2

matrices evaluated at xt = x̂t (t = s, . . . , s+m−1). Further, let α̂t+1 = Φα̂t +Σηη x̂t(t = s, . . . , s + m − 1) with α̂s = hs − µ.
3. Set bs = 0 and B̂s+m+1 = O. Compute0
Dt = Ât − B̂t D−1
t−1 B̂t ,bt = dˆt − B̂t D−1
t−1 bt−1 ,t = s + 1, . . . , s + m.

4. Let0ŷt = α̂t + D−1
t (bt + B̂t α̂t+1 ),−1/20

Gt = [Dt1/2

where Dt0
, D−1
t B̂t Σηη ],0
Zt = Ip + D−1t B̂t Φ,t = s + 2, . . . , s + m,
1/2

1/20

denotes a Choleski decomposition such that Dt = Dt Dt.

5. Implement the disturbance smoother (Koopman (1993)) to obtain {x̂t }s+m−1
, the
t=s
mode of the conditional posterior density of {xt }s+m−1
for the model (18) and (19).
t=s
If the the mode converges (however, usually several iterations will be sufficient toconstruct a proposal distribution), save ŷt , Zt and Gt . Otherwise, go to Step 2.
Then we apply the simulation smoother (e.g. de Jong and Shephard (1995), Durbin andKoopman (2002)) to generate a candidate {x†t }s+m−1
from this state space model fort=s

9Metropolis-Hastings algorithm. We accept a candidate with probability(
min 1,f ({x†t }s+m−1)f ∗ ({xt }s+m−1
)
t=s
t=s
f ({xt }s+m−1)f ∗ ({x†t }s+m−1
)
t=s
t=s

)
,

where f ({xt }s+m−1) and f ∗ ({xt }s+m−1) denote conditional posterior densities of {xt }s+m−1
t=s
t=s
t=s
for MESV model and the linear Gaussian state space model (18)–(19).
3.2.2

Generation of Σ

Since the p × p leading principal submatrix of Σ is an identity matrix, we first generate Σ22
and then sample vec(Σ21 ) conditional on Σ22 , to conduct MH algorithm using the propertyof Wishart distribution (see, e.g., Theorem 3.3.9 of Gupta and Nagar (2000)).
Let
R−1 = 



R21

 P

n−1
ε ε0
 =  P t=1 t t
n−1
0
R22
t=1 ηt εt

R11 R12

Pn−1



Pn−1

,

0
t=1 εt ηt0
t=1 ηt ηtwhere R11 , R12 = R210 , R22 are p × p, p × q, q × q matrices, εt = exp(−Ht /2)yt and ηt =ht+1 − µ − Φ(ht − µ). Using tr(AB) = vec(A0 )0 vec(B) and vec(AXB) = (B0 ⊗ A)vec(X),for X(n × n), A(m × n) and B(n × m), the joint conditional posterior probability densityof Σ12 and Σ22 is obtained as follows.
π(Σ12 , Σ22 |Φ, h, Yn ) = π(Σ22 |Φ, h, Yn )π(vec(Σ21 )|Σ22 , Φ, h, Yn )
½
¾
½
¾
¢
¢
n0 −q−1
n−1
1 ¡
1 ¡
22
∝ h(Σ) × |Σ|− 2 exp − tr R−1 Σ × |Σ22 | 2 exp − tr R−1
Σ0
2
2
×fN (vec(Σ21 )|vec(Σ22 ∆0 ), Ω0 ⊗ Σ22 )·¸ª
n +n−1−q−1
1 © 22
−1
22 0
21
11
−1
12
22
2
∝ h(Σ) × |Σ |exp − tr (R − R R
R + R0 )Σ
2¾½
1
× exp − vec(Σ21 + Σ22 R21 R11 −1 )0 (R11 ⊗ Σ22 −1 )vec(Σ21 + Σ22 R21 R11 −1 )2
×fN (vec(Σ21 )|vec(Σ22 ∆0 ), Ω0 ⊗ Σ22 )½¾
n −q−11 ¡ −1 22 ¢22 1 2∝ h(Σ) × |Σ |exp − tr R1 Σ× fN (vec(Σ21 )|vec(Σ22 ∆1 ), Ω1 ⊗ Σ22 )
2

10

−1 0
−1
−1 0
22
where n1 = n0 + n − 1, R−1
1 = R − ∆1 Ω1 ∆1 + R0 + ∆0 Ω0 ∆0 and

¡ 11
¢−1
,
R + Ω−10¡¢Ω1∆1 = −R21 + ∆0 Ω−1
0
½
¾
1
1
h(Σ) = |Σ0 |− 2 exp − (h1 − µ)0 Σ−1
(h
−
µ)
.
1
0
2
Ω1 =

Thus, we generate a candidate Σ† in three steps.
1. Draw Σ22† ∼ W(n1 , R1 ).
2. Draw vec(Σ21† )|Σ22† ∼ Npq (vec(Σ22† ∆1 ), Ω1 ⊗ Σ22† ).
3. Compute Σ†εη = −Σ12† Σ22†−1 and Σ†ηη = Σ22†−1 + Σ†ηε Σ†εη .
and accept it with probability
½
min
3.2.3

¾
h(Σ† ),1 .
h(Σ)

Generation of (µ, φ)Generation of µ. The conditional posterior distribution of µ isµ|Σ, Φ, h, Yn ∼ Nq (m∗1 , V1∗ ),
where
Ã
V1∗ =

Σ0−1 +

!−1

n−1
X

∗ −1
(Ip − Φ)0 Σ−1h (Ip − Φ) + V0

,

t=1

m∗1

=

V1∗

"n−1
#
X
−1
∗ −1 ∗
0 −1(Ip − Φ) Σh {ht+1 − Φht − Σεη exp(−Ht /2)yt } + Σ0 h1 + V0 m0 .
t=1

Generation of φ. Let A =Pn−10
t=1 (ht − µ)(ht − µ) , B =

Pn−1

12
0
t=1 {(ht − µ)yt exp(−Ht /2)Σ +(ht − µ)(ht+1 − µ)0 Σ22 } and b denote a vector whose i-th element is equal to the (i, i)-thelement of B. Then the conditional posterior probability density function of φ is½¾122π(φ|Σ, α, Yn ) ∝ exp − tr(ΦΣ ΦA) − 2tr(ΦB) × k(φ)
2∝ fN (φ|µffi , Σffi ) × k(φ),½¾
1 0 −1− 12
k(φ) = |Σ0 | π(φ) exp − α1 Σ0 α1 ,
2
11

22
where µffi = Σffi b, Σ−1
ffi = Σ ¯ A and ¯ denotes the Hadamard product. To sample φfrom its conditional posterior distribution using MH algorithm, we generate a candidatefrom a truncated normal distribution over the region R, φ† ∼ T N R (µffi , Σffi ), R = {φ :|φj | < 1, j = 1, . . . , p} and accept it with probability min{k(φ† )/k(φ), 1}.

3.3

Associated particle filter

We describe the associated auxiliary particle filter introduced by Pitt and Shephard (1999)
for the MESV model to compute the log likelihood function. In Section 5, we use thisalgorithm to calculate DIC (deviance information criterion proposed by Spiegelhalter, Best,
Carlin, and van der Linde (2002)). Let½¾
1 0
f (yt |ht ) = (2π) | exp (Ht ) | exp − yt exp(−Ht )y ,
2
½
¾
q
1
1
−2
−2
0 −1
f (ht+1 |yt , ht , θ) = (2π) |Σh | exp − (ht+1 − µh,t+1 ) Σh (ht+1 − µh,t+1 ) ,
2
− p2

− 12

and f (ht |Yt , θ) denote a conditional density of ht given (Yt , θ). Then the conditional joint
density function of ht+1 , ht , given (Yt+1 , θ) is
f (ht+1 , ht |Yt+1 , θ) ∝ f (yt+1 |ht+1 )f (ht+1 |yt , ht , θ)f (ht |Yt , θ).
We first construct an importance function to sample from the conditional joint distribution.
Let fˆ(ht |Yt , θ) denote a discrete probability mass function approximating f (ht |Yt , θ) andg(ht+1 , hit |Yt+1 , θ) ∝ f (yt+1 |µih,t+1 )f (ht+1 |yt , hit , θ)fˆ(hit |Yt , θ)
∝ f (ht+1 |yt , hit , θ)g(hit |Yt+1 , θ),wheref (yt+1 |µih,t+1 )fˆ(hit |Yt , θ)
,
PI
j
ˆ j
j=1 f (yt+1 |µh,t+1 )f (ht |Yt , θ)
½
¾
1 0
− p2
i
− 12
i
i
f (yt+1 |µh,t+1 ) = (2π) | exp(Mh,t+1 )| exp − yt+1 exp(−Mh,t+1 )yt+1 ,2g(hit |Yt+1 , θ) =µih,t+1 = vech(Mih,t+1 ) = µ + Φ(hit − µ) + Σηε exp(−Hit /2)yt .
Using this importance function, we implement the auxiliary particle filter as follows.
1. Set t = 1.
12

(a) Generate hi1 ∼ N (µ, Σ0 ) (i = 1, . . . , I).
P(b) Compute wi = f (y1 |hi1 , θ) and save w̄1 = I1 Ii=1 wi .
P(c) Let fˆ(hi1 |Y1 , θ) = π1i = wi / Ij=1 wj (i = 1, . . . , I).
¡¢
2. Generate hit+1 , hit ∼ g(ht+1 , hit |Yt+1 , θ) (i = 1, . . . , I):¡¢(a) Compute µih,t+1 = µ + Φ(hit − µ) + Σηε exp −Hit /2 yt .
(b) Generate hit ∼ g(hit |Yt+1 , θ).
(c) Generate hit+1 ∼ f (ht+1 |yt , hit , θ).
Then computef (yt+1 |hit+1 )f (hit+1 |yt , hit , θ)fˆ(hit |Yt , θ)¡¢,g hit+1 , hit |y t+1 , θf (yt+1 |hit+1 )fˆ(hit |Yt , θ)
, i = 1, . . . , I,g(hit |Yt+1 , θ)wi ==and save

I

wt =

1X
wi .
I
i=1

P
i
Further let fˆ(hit+1 |Yt+1 , θ) = πt+1
= wi / Ij=1 wj (i = 1, . . . , I).
3. Set t ← t + 1. Go to Step 2.
Then, as I → ∞, we obtain

Pn

p Pn
t=1 log w t →t=1 log f (yt |Yt−1 , θ).

13

4

Illustrative example with simulated data

This section shows the efficiency of our proposed method using a simulated data. To
replicate the dynamics of the stock return series, we use the MESV model with
 0.5, i = j,µij =,
 0.2, i 6= j, 0.97, i = j,φij =
 0.85, i 6= j, −0.3, i = j = k,
ρiE(j,k),εη ≡ Corr(εit , Et (j, k)) = −0.1, otherwise, 0.3, i = j 6= k = l
ρE(i,j)E(k,l),ηη ≡ Corr(Et (i, j), Et (k, l)) = 0.1, otherwise, 0.2, i = jp
σE(i,j),ηη ≡ V ar(Et (i, j)) = 0.15, otherwise,
which are based on typical values in our empirical studies where E(1, 1) = 1, E(2, 1) = 2,
E(3, 1) = 3, E(2, 2) = 4, E(3, 2) = 5 and E(3, 3) = 6.
We generate n = 4, 000 observations with p = 3 and q = 6. Prior distributions are
assumed to be as follows.
µ ∼ Nq (0, 5Iq ), B(20, 3/2), i = j,
φij + 1∼
 B(1, 1),
2
i=
6 j,
Σ22 ∼ W(6, (6Σ22∗ )−1 ),vec(Σ21 )|Σ22 ∼ Npq (0, (5Ip ) ⊗ Σ22 ),where Σ22∗ is a true covariance matrix satisfying E(Σ22 ) = Σ22∗ . The mean and the standard deviation of the prior distribution of φii , j = 1, 2, 3 are set 0.86 and 0.11 respectively.
Using the multi-move (single-move) sampler, we draw 110,000 (550,000) posterior samples
and discard the first 10,000 (50,000) samples as burn-in periods.
Parameter estimation. Tables 1, 2, and 3 show the estimation summaries for all parameters
via the multi-move sampler. The posterior means and 95% credible intervals suggest that theestimates are sufficiently close to true values, which indicates that our proposed estimation
algorithm works well. The inefficiency factors are also shown for the multi-move sampler

14and the single-move sampler.
The inefficiency factor is the ratio of the numerical variance of the estimate from the
MCMC samples relative to that from hypothetical uncorrelated samples, and is definedP
as 1 + 2 ∞s=1 ρs where ρs is the sample autocorrelation at lag s. It suggests the relative
number of correlated draws necessary to attain the same variance of the posterior sample
mean from the uncorrelated draws (see e.g. Chib (2001)).
Table 1: Posterior means, 95% credible intervals, and inefficiency factors.
Param.

True
0.97

φij
0.85

0.5
µij
0.2

0.2
σE(i,j),ηη
0.15

ij

Mean 95% interval

11
22
33
21
31
32
11
22
33
21
31
32
11
22
33
21
31
32

0.967 [0.955, 0.978]0.964 [0.951, 0.975]
0.975 [0.965, 0.984]
0.802 [0.701, 0.877]
0.841 [0.734, 0.910]
0.837 [0.751, 0.900]
0.572 [0.390, 0.755]
0.456 [0.266, 0.649]
0.485 [0.225, 0.744]
0.201 [0.159, 0.243]
0.203 [0.158, 0.247]
0.204 [0.156, 0.252]
0.189 [0.158, 0.223]
0.216 [0.185, 0.251]
0.211 [0.181, 0.242]
0.158 [0.118, 0.204]
0.145 [0.106, 0.194]
0.175 [0.129, 0.230]p
σE(i,j),ηη = Var(Et (i, j))Inefficiency
multi [single]
149
[334]
119
[115]
105
[437]
288
[795]
511
[1622]
332
[435]
6
[38]
6
[75]
3
[42]
20
[465]
17
[355]
18
[304]
285
[572]
219
[306]
180
[657]
460
[1182]
646[1648]
438[761]

The inefficiency factors for the single-move sampler are about three times larger thanthose for the multi-move sampler. This implies that our proposed multi-move sampler ismore efficient than the single-move sampler as we expected.

15

Table 2: Posterior means, 95% credible intervals and inefficiency factors.
ρiE(j,k),εη = Corr(εit , Et (j, k))
True
−0.3

−0.1

i jk

Mean

95% interval

1 11
2 22
3 33
1 21
1 31
1 22
1 32
1 33
2 11
2 21
2 31
2 32
2 33
3 11
3 21
3 31
3 22
3 32

-0.277-0.255-0.3740.046-0.041-0.113-0.0350.005-0.018-0.017-0.172-0.072-0.067-0.092-0.194-0.118-0.080-0.027[-0.388,-0.159][-0.372,-0.134][-0.479,-0.264][-0.100, 0.188][-0.190, 0.107][-0.226, 0.002][-0.163, 0.092][-0.114, 0.123][-0.147, 0.111][-0.172, 0.131][-0.328,-0.022][-0.207, 0.061][-0.192, 0.057][-0.217, 0.030][-0.341,-0.047][-0.263, 0.028][-0.192, 0.035][-0.155, 0.104]

Inefficiency
multi [single]
61
[133]
84
[126]
81
[220]
73
[164]
96
[218]
34
[143]
63
[206]
66
[155]
60
[231]
94
[164]
86
[289]
65
[189]
71
[184]
61
[159]
87
[113]
91
[293]
51
[174]
72
[149]

Table 3: Posterior means, 95% credible intervals and inefficiency factors.
ρE(i,j)E(k,l),ηη = Corr(Et (i, j)), Et (k, l))
True
0.3

0.1

ij kl

Mean

95% interval

11 22
11 33
22 33
11 21
11 31
11 32
21 31
21 22
21 32
21 33
31 22
31 32
31 33
22 32
32 33

0.248
0.332
0.266
-0.055
0.150
0.064
0.092
-0.034-0.127
0.203
0.055
0.041
0.213
-0.0060.053[ 0.080, 0.411][ 0.170, 0.481][ 0.104, 0.420][-0.266, 0.177][-0.066, 0.355][-0.135, 0.267][-0.199, 0.377][-0.239, 0.188][-0.386, 0.142][-0.028, 0.424][-0.173, 0.269][-0.192, 0.275][ 0.010, 0.413][-0.203, 0.189][-0.134, 0.239]

16

Inefficiency
multi [single]
102
[340]
135
[382]
68
[186]
184
[971]
215
[662]
138
[109]
301
[110]
205
[133]
239
[100]
177
[106]
234
[111]
213
[161]
224
[931]
137
[536]
138
[830]

5

Application to trivariate asset returns data

5.1

Data

This section applies our proposed MESV model to returns of three assets: (1) Tokyo stock
price index (TOPIX), (2) the Japan government bond clean price index (JGB) provided
by Thomson Reuters Datastream), and (3) the currency exchange rate of Japanese Yen to
U.S. Dollar (Yen/USD) announced by the Federal Reserve Bank at noon in New York. We
excluded those days when at least one of three observations is not reported. The sample
period is from January 4, 1995 to July 30, 2010 for a total of 3710 observations. Figure
1 shows the time series plot of three returns which are 100 times the differences of thelogarithm of the asset values.

5

TOPIX

0
−5
Jan. 1995
1

Feb.1997

Apr.1999

May.2001

Jul.2003

Sep.2005

Oct.2007

Dec.2009

Feb.1997

Apr.1999

May.2001

Jul.2003

Sep.2005

Oct.2007

Dec.2009

Apr.1999

May.2001

Jul.2003

Sep.2005

Oct.2007

Dec.2009

JGB

0−1
Jan. 1995
2.5

Yen/Doller

0.0
−2.5
Jan. 1995

Feb.1997

Figure 1: TOPIX, JGB and Yen/USD

5.2
5.2.1

Estimation results
Univariate SV and MSV models

We first estimate the (univariate) SV model with leverage effect defined by (1)–(3) with
p = 1 for each series, and the constant correlation multivariate SV (MSV) model definedby (1)–(3) with p = 3. For the SV model, the prior distributions for φi and Σi are assumed17

to be
φi + 1
2∼ B(20, 1.5),Σ∗ = 

Σi ∼ IW(5, (5Σ∗ )−1 ),
10

0 0.04,i = 1, 2, 3.

For the MSV model, we assume prior distributions such thatφi + 12

Σ ∼ IW(6, (6Σ∗ )−1 ),∼ B(20, 1.5), i = 1, . . . , 5,

where


Σ∗ = 


1.52

0

0

0

0.5I2

0



.


0.22 (0.3I3 + 0.713 103 )Tables 4 and 5 show summary statistics of posterior distributions of the parameters for
the SV and the MSV models. In both models, the posterior means of the autoregressive
parameters (φj ’s) are very high (between 0.938 and 0.972 in the SV model and between
0.948 and 0.966 in the MSV model) showing that volatilities are highly persistent.
Table 4: Univariate SV model.
Posterior means and 95% credible intervals.
Param.
φi
σiε

σiη

ρi

i
1
2
3
1
2
3
1
2
3
1
2
3

Mean
0.971
0.972
0.938
1.177
0.267
0.652
0.177
0.184
0.239
-0.511-0.176-0.201

95% interval[0.960, 0.981][0.959, 0.983][0.897, 0.966][1.078, 1.285][0.240, 0.298][0.609, 0.699][0.148, 0.207][0.150, 0.220][0.175, 0.321][-0.604,-0.408][-0.293,-0.052][-0.316,-0.079]

1: TOPIX, 2: JGB, 3: Yen/USD.

Also, the leverage effects are estimated to be negative as in the previous literature wherethat of the stock return (−0.511 in the SV model and −0.445 in the MSV model) is muchstronger than those of the bond return and the foreign exchange return (−0.176 and −0.201
in the SV model, −0.177 and −0.167 in the MSV model). In the MSV model, the cross18

leverage effect, ρ12,εη , (from the stock return to the bond return volatility) is estimated to
be positive (0.189), while the opposite effect, ρ21,εη , (from the bond return to the stock
return volatility) is not credible in the sense that its 95% credible interval doesn’t include
zero. This implies the increase in the stock return at time t causes the high volatility in thebond return at time t + 1, but the fall of the bond return seems to have a limited impact
on the stock return volatility.
Table 5: MSV model.
Posterior means and 95% credible intervals.
Param.
φi
σi,εε
σi,ηη

ρii,εη

i
1
2
3
1
2
3
1
2
3
1
2
3

Mean
0.965
0.966
0.948
1.212
0.275
0.662
0.191
0.210
0.225
-0.445-0.177-0.167

Param.

95% interval[0.953, 0.975][0.952, 0.978][0.918, 0.969][1.115, 1.322][0.249, 0.305][0.615, 0.716][0.163, 0.224][0.174, 0.249][0.173, 0.293][-0.543,-0.339][-0.281,-0.067][-0.282,-0.045]

ρij,εερij,ηη

ρij,εη

ij
12
13
23
12
13
23
12
13
21
23
31
32

Mean
-0.3120.062-0.014
0.196
0.559
0.390
0.189
-0.0030.091
0.035-0.036-0.05095% interval[-0.343,-0.281][ 0.028, 0.095][-0.048, 0.020][ 0.024, 0.363][ 0.402, 0.694][ 0.221, 0.552][ 0.073, 0.301][-0.127, 0.117][-0.021, 0.200][-0.077, 0.145][-0.152, 0.079][-0.167, 0.067]

1: TOPIX, 2: JGB, 3: Yen/USD.
ppσi,εε = Var(εit ), σi,ηη = Var(ηit ),ρij,εε = Corr(εit , εjt ), ρij,εη = Corr(εit , ηjt ), ρij,ηη = Corr(ηit , ηjt ).

5.2.2

MESV model

Next, we estimate the MESV model assuming that prior distributions are
φij + 12∼


 B(20, 3/2),

if i = j, B(33/4, 11/4), otherwise,µ ∼ Nq (0, 5Iq ),Σ22 ∼ W(6, (6Σ22∗ )−1 ),vec(Σ21 )|Σ22 ∼ Npq (0, (5Ip ) ⊗ Σ22 ).

∗∗The (i, j)-th element of Σ22∗−1 is ρ∗E(k,l)E(m,n),ηη σE(k,l),ηη
σE(m,n),ηη
such thatρ∗E(k,l)E(m,n),ηη 0.6, if k = l 6= m = n,=
 0.2, otherwise,∗
σE(k,l),ηη= 0.2,

19

1 ≤ l ≤ k ≤ 3.

The mean and the standard deviation of the prior distribution of φij , i 6= j are set 0.5 and
0.5 respectively.
We draw 220,000 samples for the multi-move sampler discarding the first 20,000 samples
as a burn-in-period. The number of blocks is set to 185 based on several trials.
Table 6: MESV model.
Posterior means, standard deviations, 95% credibleintervals, and inefficiency factors.

φij

µij

σE(i,j),ηη

ij
11
22
33

Mean
0.968
0.968
0.953

Stdev
0.005
0.006
0.012

95% interval[0.958, 0.977][0.955, 0.979][0.927, 0.972]

IF
141
94
648

21
31
32
11
22
33

0.860
0.845
0.590
0.286
-2.923-0.929

0.050
0.060
0.199
0.087
0.113
0.079

[0.738, 0.923][0.699, 0.924][0.051, 0.826][ 0.114, 0.457][-3.143,-2.699][-1.084,-0.773]

693
560
741
6
11
16

21
31
32
11
22
33

-0.222
0.054
0.007
0.176
0.211
0.204

0.018
0.023
0.017
0.014
0.018
0.027

[-0.258,-0.186][ 0.009, 0.099][-0.027, 0.041][0.150, 0.204][0.176, 0.249][0.156, 0.263]

24
25
75
253
277
645

21
31
32

0.106
0.133
0.144

0.016
0.022
0.022

[0.080, 0.144][0.097, 0.185][0.107, 0.193]

742
828
801

1:TOPIX, 2:JGB, 3:Yen/USD
p
σE(i,j),ηη = V ar(Et (i, j))

The estimation results are summarized in the Tables 6, 7 and 8. We notice that the parameters of the diagonal elements of Ht (the 1st, 4th, and 6th elements of ht ) are similar to
those of MSV models. The autoregressive parameters of log volatilities, (φ11 , φ22 , φ33 ), are(0.968, 0.968, 0.953) while (φ1 , φ2 , φ3 ) for MSV models are (0.965, 0.966, 0.948). The poste2 , log σ 2 , log σ 2 )rior means of (µ11 , µ22 , µ33 ) are (0.286, −2.923, −0.929), while (log σ1,εε
2,εε
3,εε

evaluated at the posterior means in Table 5 are (0.385, −2.582, −0.825). Further, the estimates of standard deviations of the ηit , (σE(1,1),ηη , σE(2,2),ηη , σE(3,3),ηη ), are (0.176, 0.211,
0.204), while those of (σ1,ηη , σ2,ηη , σ3,ηη ) are (0.191, 0.210, 0.225). Regarding the leverage effects, the estimates of (ρ1E(1,1),εη , ρ2E(2,2),εη , ρ3E(3,3),εη ) are (−0.467, −0.109, −0.166),
while (ρ11,εη , ρ22,εη , ρ33,εη ) for the MSV models are (−0.445, −0.177, −0.167).

20

Table 7: MESV model. Posterior means, standard deviations,
95% credible intervals and inefficiency factors.
ρiE(j,k),εη = Corr(εit , Et (j, k))
i jk
1 11
2 22
3 33
1 21
1 31
1 22
1 32
1 33
2 11
2 21
2 31
2 32
2 33
3 11
3 21
3 31
3 22
3 32

Mean
-0.467-0.109-0.166-0.100-0.050
0.218
0.058
0.026
-0.0050.023
0.012-0.0460.055-0.0240.078-0.128-0.061
0.020

Stdev
0.054
0.059
0.063
0.065
0.073
0.060
0.085
0.065
0.061
0.068
0.072
0.080
0.063
0.061
0.065
0.074
0.062
0.085

95% interval[-0.567,-0.358][-0.224, 0.008][-0.287,-0.041][-0.230, 0.027][-0.191, 0.094][ 0.097, 0.332][-0.108, 0.225][-0.107, 0.150][-0.125, 0.114][-0.108, 0.160][-0.130, 0.152][-0.204, 0.109][-0.069, 0.177][-0.143, 0.095][-0.051, 0.206][-0.273, 0.016][-0.183, 0.061][-0.148, 0.188]

IF
90
64
87
71
89
65
162
99
70
74
74
120
76
101
69
109
67
215

1:TOPIX, 2:JGB, 3:Yen/USD

Table 8: MESV model. Posterior means, standard deviations,
95% credible intervals and inefficiency factors.
ρE(i,j)E(k,l),ηη = Corr(Et (i, j), Et (k, l))
ij kl
11 22
11 33
22 33
11 21
11 31
11 32
21 31
21 22
21 32
21 33
31 22
31 32
31 33
22 32
32 33

Mean
0.009
0.466
0.263
0.151
0.155
-0.2330.020
0.081-0.023-0.0100.045-0.151
0.134
0.062-0.145

Stdev
0.091
0.083
0.097
0.088
0.099
0.131
0.108
0.096
0.137
0.104
0.109
0.140
0.109
0.144
0.136

95% interval[-0.168, 0.187][ 0.293, 0.616][ 0.064, 0.443][-0.025, 0.320][-0.041, 0.350][-0.459, 0.056][-0.185, 0.233][-0.108, 0.267][-0.296, 0.240][-0.211, 0.197][-0.175, 0.252][-0.419, 0.135][-0.079, 0.348][-0.217, 0.346][-0.399, 0.139]

1:TOPIX, 2:JGB, 3:Yen/USD

21

IF
155
221
237
151
181
362
230
154
369
202
193
395
225
386
472

Estimated volatilities and correlations. However, as mentioned in Section 2, the parameterssuch as Φ, M and Σεη in the MESV model do not always correspond to those of the stock,
the bond, and the exchange rate as they are in the MSV model. Thus, to interpret theestimation results of the MESV model in more intuitive way, we consider the posteriormeans of time-varying volatilities of each series, dynamic correlations among three returns
and news impact curves using MCMC simulation technique.
First, we consider the posterior means of the square root of the time-varying variances asshown in Figure 2. The estimated volatility series of the TOPIX returns sharply increasedin September 2008, corresponding to the financial crisis during which Lehman Brothers filed
for Chapter 11 bankruptcy protection. The volatilities of the JGB index returns increasedin December 1998 and in September 2003 when the index dropped (or equivalently, the JGB
interest rate ran up) in both periods. In December 1998, the JGB is supplied excessively
because the Japanese government issued a large amount of JGB for the economic-stimulus
measure, and the Ministry of Finance Japan announced to stop buying the new bonds inthis month. Moreover, the Moody’s downgraded the JGB rating in November 1998. In
this month, the Nippon Credit Bank was brought under government control because of alarge amount of the bad debt. In the mid 2003s, following the increase of the US bondinterest rate and the economic boom, the deflationary concerns of Japan toned down. The
expectation for the lifting of the zero-interest-rate policy of the Bank of Japan, the JGB
interest rate shot up in June and September of 2003. The volatilities of the Yen/USD
increase after the August 1998 when the ruble devaluation and the Long Term CapitalManagement report a large loss. Especially, the USD fell from 135.6 yen to 117 yen in fivedays of early October 1998. This is just after the decision of the monetary relaxation policyin the USA on September 29th and issuing the G-7 communique which urged the injectionof taxpayers’ money to financial institutions in Japan on October 5th.
Next, we investigate the posterior means of dynamic correlations among three returns
as shown in Figure 3. These correlations are computed using the MCMC samples of thecovariance and the variances that are elements of exp(Ht ) which is the matrix exponential
transformation of the log volatility matrix Ht . The correlations between the stock and the
JGB returns largely fluctuate taking negative values where they drop to less than −0.86
in January 2008 during the downturn of the stock market. The correlations between thestock market and the exchange rate returns fluctuates around zero. It is noted that thistakes negative values during Asian crisis period from July 1997 to August 1998. In thisperiod, the yen kept weakening and Japanese stock prices dropped. The JGB returns and22the exchange rate returns seem to have no correlation throughout the sample period.

stdev(y1)42Jan. 1995
1.00

Feb. 1997

Apr. 1999

May. 2001

Jul. 2003

Sep. 2005

Oct. 2007

Dec. 2009

Apr. 1999

May. 2001

Jul. 2003

Sep. 2005

Oct. 2007

Dec. 2009

Apr. 1999

May. 2001

Jul. 2003

Sep. 2005

Oct. 2007

Dec. 2009

stdev(y2)

0.75
0.50
0.25
Jan. 1995
2

Feb. 1997

stdev(y3)1Jan. 1995

Feb. 1997

Figure 2: Posterior means of the square root of time-varying variances. Top: TOPIX,
middle: JGB, bottom: Yen/USD.

1.0

Corr(y1,y2)

0.5
0.0
−0.5
Jan. 1995
1.0

Feb. 1997

Apr. 1999

May. 2001

Jul. 2003

Sep. 2005

Oct. 2007

Dec. 2009

Apr. 1999

May. 2001

Jul. 2003

Sep. 2005

Oct. 2007

Dec. 2009

Apr. 1999

May. 2001

Jul. 2003

Sep. 2005

Oct. 2007

Dec. 2009

Corr(y1,y3)

0.5
0.0
−0.5
Jan. 1995
1.0

Feb. 1997

Corr(y2,y3)

0.5
0.0
−0.5
Jan. 1995

Feb. 1997

Figure 3: Posterior means of time-varying correlations. Top: (TOPIX, JGB), middle:(TOPIX, Yen/USD), bottom: (JGB, Yen/USD).

23

Dynamic principal components implied by the MESV model. Furthermore, we can interpret our estimation results for MESV model in terms of the principal component analysis.
Using the diagonalization to calculate the matrix exponential, we calculate the principal
components of yt . Since Ht = Ut Λt U0t where Λt = diag(λ1t , λ2t , λ3t ) and Ut U0t = I3 , weobtainU0t yt = diag(exp(λ1t /2), exp(λ2t /2), exp(λ3t /2))ε∗t ,(20)

where ε∗t := U0t εt ∼ N (0, I3 ), and we assume λ1t ≥ λ2t ≥ λ3t . The i-th element of U0t yt in(20) can be interpreted as the i-th principal component of the returns at time t. Moreover,
the row vector of the loading matrix U0t can be a interpreted as a weight vector for the
portfolio which represents the stock market principal component. Top and middle panelsof Figure 4 show time series plots of the principal components, U0t yt , t = 1, . . . , n and theirstandard deviations.

10

U’y_1
U’y_3

U’y_2

0
−10
Jan.1995

Feb.1997

Apr.1999

May.2001

Jul.2003

Sep.2005

Oct.2007

Dec.2009

Feb.1997

Apr.1999

May.2001

Jul.2003

Sep.2005

Oct.2007

Dec.2009

Feb.1997

Apr.1999

May.2001

Jul.2003

Sep.2005

Oct.2007

Dec.2009

5.0

2.5

Jan.1995
1.0
0.8
0.6
0.4
0.2
Jan.1995

Figure 4: Posterior means of principal components series (top) and their standard deviations(middle), and contribution ratio series (bottom).

24

The first principal component (blue lines) fluctuates sharply and its variance shot up whenat least one of the volatilities of the original series increase around 1998, 2003 and 2008. The
bottom panel of Figure 4 shows the contribution ratio of each principal component, defined
1t )
2t )
3t )
by P3exp(λ, P3exp(λ, P3exp(λfor t = 1, . . . , n. The ratio of the first principal
exp(λ )
exp(λ )
exp(λ )
i=1

it

i=1

it

it

i=1

component (blue part) increases when one of three return series became more volatile. On
the other hand, the contribution ratio of the second series (green part) decreases in thosevolatile periods. The third principal component (red part) has the ratio less than 10 % overalmost all the sample period. This implies that most of the variation of the original threeseries can be explained by two principal components.
News impact curves. Finally, to show how the shocks in the returns at time t affect thevolatilities at time t + 1, we describe the news impact curve following Engle and Ng (1993).
Similar ideas for stochastic volatility models are discussed by Yu (2005) and Asai and
McAleer (2009). Let Ht = M and Et = O and consider the case where ht+1 = µ +Σηε exp(−M/2)yt .
y_1
y_1U
y_2L
y_3
y_3U

stdev(y1)
1.50

y_1L
y_2
y_2U
y_3L −0.2

1.25

−0.3

1.00

−0.4
−5.0
−2.5
stdev(y2)

0.0

2.5

5.0

0.30

Corr(y1,y2)

−5.0
−2.5
Corr(y1,y3)

0.0

2.5

5.0

−5.0

0.0

2.5

5.0

0.0

2.5

5.0

0.2

0.25
0.0
0.20
−5.0
−2.5
stdev(y3)

0.0

2.5

5.0

0.7

0.1

0.6

−0.1
−5.0

−2.5

−2.5

Corr(y2,y3)

0.0

2.5

5.0

−5.0

−2.5

Figure 5: Posterior mean (solid line) and 95% interval (dotted lines) of the news impact
curve∗ for the one-step-ahead conditional covariances of yt+1 when Ht = M and Et = O.
(∗ The domain of the estimated curve is restricted to the range within ±4 sample standard deviationof actual returns for each yit .)25

Figure 5 shows the posterior news impact curves on standard deviations and correlationsof yt+1 from yt , obtained by calculating the posterior means and the 95% intervals of thenormalized exp(Ht+1 ) under the various shocks of yt . The horizontal and vertical axes show
the values of yt and the standard deviations and correlations calculated from exp(Ht+1 )respectively.
The left three panels in Figure 5 show the news impacts on the standard deviationsof yt+1 caused by the elements, y1t (red lines), y2t (blue lines) and y3t (black lines). The
negative return on the i-th asset increases its own (i-th) future volatility, indicating theexistence of leverage effects. The red lines in the middle and bottom left panels imply thatthe positive TOPIX return increases the future volatility of the JGB and the exchange rate
returns. However, the impact of the TOPIX return on the exchange rate return volatility
is smaller than that on the volatility of the JGB return. The positive shocks of the JGB
return increase the future volatility of the exchange rate, while the negative shocks of theexchange rate return cause the higher future volatility of the JGB return. These results
are generally consistent with those obtained with MSV model, taking account of the 95%
credible intervals.
The right three panels in Figure 5 show the news impacts on the correlations among
y1,t+1 , y2,t+1 and y3,t+1 . The top right panel shows the impact on the correlation between
the TOPIX and JGB returns. It is noted that the correlations are strongly affected by theexchange rate return but hardly affected by the TOPIX and the JGB returns. The blackline in the middle right panel shows the impact on the correlation between the TOPIX andexchange rate returns caused by the shock of the exchange rate return. Interestingly, thesign of the correlation between the TOPIX and exchange rate returns strongly depends on
the impact of the shock. More precisely, the large positive shock (greater than one) on
the exchange rate return tends to produce the negative correlation between the TOPIX
and exchange rate return, while the small or negative shock tends to produce the positive
correlation. The sign of the correlation between the JGB and exchange rate returns also
depends on the impacts of the shocks. However, we note that the impacts by shocks of the
TOPIX and exchange rate returns are very small and the impacts by the JGB return shock
have wide 95% intervals.
It is noted that the 95% intervals of the impact on the standard deviations and correlations looks very large and that some of them include the horizontal line (i.e. no leveragecase). Thus, we calculate the posterior probability that the standard deviations and thecorrelations calculated from exp(Ht+1 ) are larger than the posterior means of the standard
26

deviations and the correlations from exp(M). It means that the news impacts are larger
than the posterior means of them at yt = 0. Figure 6 shows the posterior probabilities ofsuch events. The left three panels show the results for standard deviations. The posterior
probability for the standard deviation of the TOPIX (exchange rate) return approaches one
as the TOPIX (exchange rate) return decreases in the top (middle) left panel. The probability for the standard deviation of the JGB return also approaches one as the TOPIX return
increases in the middle left panel. The right three panels show the results for correlation
coefficients. The posterior probability for the correlation between the TOPIX and exchange
rate returns rapidly goes to one as the exchange rate return decreases in the middle rightpanel.
1.00

stdev(y1)
y_1
y_3

0.75

y_2

1.00
0.75

0.50

0.50

0.25

0.25

1.00

−5.0
−2.5
stdev(y2)

0.0

2.5

5.0
1.00

0.75

0.75

0.50

0.50

0.25

0.25

1.00

−5.0
−2.5
stdev(y3)

0.0

2.5

5.0
1.00

0.75

0.75

0.50

0.50

0.25

0.25
−5.0

−2.5

0.0

2.5

Corr(y1,y2)

5.0

−5.0
−2.5
Corr(y1,y3)

0.0

2.5

5.0

−5.0
−2.5
Corr(y2,y3)

0.0

2.5

5.0

−5.0

0.0

2.5

5.0

−2.5

Figure 6: Posterior probability that news impact curves are larger than the posterior meanof exp(M).

5.3

Extensions to fat-tailed error distributions and model comparison

Finally, we conduct a model comparison of the proposed MESV model with MSV models.
In addition to the MESV and MSV models with normal errors, we consider extended models
with fat-tailed error distribution given by
1/2

εt = ξt et ,et ∼ Np (0, Ip ),27

where ξt is a random variable which takes positive values and independent of et . We consider¡¢
two cases: (i) multivariate Student-t error with ξt−1 ∼ G ν2 , ν2 and (ii) log-Gaussian scale³ 2´
mixture error with log ξt ∼ N − δ2 , δ 2 . The extensions to these models are straightforward
and hence we omit details of MCMC algorithms (similar specifications and MCMC sampling
are also discussed in the Omori, Chib, Shephard, and Nakajima (2007) for univariate SVmodels, and Ishihara and Omori (2012) for MSV models). Thus, we consider the following
six models:• MSV-n model: MSV model with normal error distribution.
• MSV-t model: MSV model with multivariate Student-t error distribution.
• MSV-g model: MSV model with log-Gaussian scale mixture error distribution.
• MESV-n model: MESV model with normal error distribution.
• MESV-t model: MESV model with multivariate Student-t error distribution.
• MESV-g model: MESV model with log-Gaussian scale mixture error distribution.
We assume the prior distributions ν ∼ G(0.001, 0.001) and δ 2 ∼ G(1, 1) respectively.
The estimation results for ν and δ 2 are summarized in Table 9. The posterior means of ν
for the MSV and MESV models are small, suggesting fat-tailed error distributions. The
estimate of ν (δ 2 ) for the MSV model is smaller (larger) than that of the MESV model,
probably because the MSV model fails to capture the dynamics of time-varying correlations.
Other parameter estimates are similar to those of models with normal error and hence areomitted.
Table 9: The estimation results of ν and δ 2
ModelParam.

Mean

Stdev

95% intervalParam.

Mean

Stdev

95% interval

MESV
MSV

ν

15.8
11.8

2.4
1.2

[11.8, 21.5][ 9.8, 14.5]

δ2

0.17
0.34

0.02
0.02

[0.13, 0.22][0.30, 0.38]

We compute the DIC (deviance information criterion) for the model comparison defined by
DIC = E„|Yn [D(θ)] + pD ,pD = E„|Yn [D(θ)] − D(E„|Yn [θ]),28

D(θ) = −2 log f (Yn |θ) + Cy ,

where Cy is a constant term which depends only on the dataset Yn . Since it cancels out inall calculations that compare different models, we set Cy = 0 for convenience. To estimate
1 PM(m) ), where we set M = 100, andEθ|Yn [D(θ)], we use a sample analogue M
m=1 D(θθ (m) s are resampled from the posterior samples generated by the MCMC method. To
calculate D(E„|Yn [θ]), which equals to D(θ) evaluated at the posterior mean, we implement
an auxiliary particle filter to compute the log-likelihood ordinate log f (Yn |θ), where we set
the number of particles I = 10, 000 for the MSV and the MESV models. We repeat thisprocedure ten times to obtain the numerical standard error.
Table 10: The averages of DIC estimates, their standard errors,
the maximum and the minimum of DIC values.
Model
MESV-tMESV-g
MESV-nranking1
2
3

DIC
20066.7
20087.4
20090.1(s.e.)
(1.7)
(1.8)
(1.5)

DICmax20074.2
20099.6
20098.4

DICmin20058.3
20079.6
20084.1

MSV-t
MSV-g
MSV-n

4
5
6

20118.4
20197.3
20201.2

(1.3)
(3.2)
(1.0)

20124.320214.3
20206.3

20112.4
20180.5
20196.0

Table 10 shows the averages of DIC, their standard errors, and the maximum and theminimum of DIC values computed for six competing models. The DIC values for the MESV
models are much smaller than those for the MSV models, and hence the MESV models outperform the MSV models. Among MESV models, models with fat-tailed error outperformthe model with normal error, and the model with multivariate-t error distribution has thesmallest DIC. This empirical study shows that our proposed model with multivariate-t error
distribution performs quite well to describe the multivariate asset returns data.

6

Conclusion

In this paper, we extend the MSV model to allow the time-varying correlations and propose an efficient MCMC algorithm using a multi-move sampler. To sample a block of statevectors, we construct a proposal density using the normal approximation via a Taylor expansion of the logarithm of the target posterior density for the MH algorithm where theexpectations of Hessian matrices are derived analytically. Moreover, to calculate the loglikelihood, we describe an auxiliary particle filter. An empirical analysis is presented using29

three returns of the TOPIX, the Japanese bond price index and the Yen/USD exchange
rate. The correlation between returns of the TOPIX and the Japanese bond index is foundto be time-varying. In contrast, the correlation between returns of the Japanese bond priceindex and the Yen/USD exchange rate is shown to be stable and less volatile. The positive cross leverage effects from the TOPIX on the Japanese bond price index is also found.
The news impact curves for the MESV model are presented and investigated in detail. A
model comparison between the MESV model with constant correlation MSV models including heavy-tailed error models is conducted. The MESV model with multivariate Student-t
distributed error is found to outperform other models based on DIC.

Acknowledgment
The authors are grateful to Siddhartha Chib, Toshiaki Watanabe, Teruo Nakatsuma, andYoshi Baba for helpful comments and discussions. This work is supported by the Research
Fellowship (DC1) from the Japan Society for the Promotion of Science and the Grants-inAid for Scientific Research (A) 21243018 and (A) 24243031 from the Japanese Ministry of
Education, Science, Sports, Culture and Technology, and Global COE Program ”Research
Unit for Statistical and Empirical Analysis in Social Sciences” from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. The computational
results are generated using Ox (Doornik (2006)).

AppendixA

Multi-move sampler

The logarithm of the full conditional joint density of {xt }s+m−1
excluding constant terms
t=s
is given bylog f ({xt }s+m−1|hs , hs+m+1 , ys , . . . , ys+m ) = −
t=s

30

s+m−1
1 X 0
xt xt + L,
2 t=s

(21)

where
L =

s+m
X

1 0
lt − ηs+m
Σ−1ηη ηs+m I(s + m < n),
2
t=s

(22)

ηs+m = hs+m+1 − µ − Φ(hs+m − µ),(23)xt = Σ−1/2ηη (ht+1 − µ − Φ(ht − µ)).
Let αt = ht − µ and α̂t = ĥt − µ where ĥt is ht evaluated at xt = x̂t . Consider theapproximation via the second order Taylor expansion around the xt = x̂t (t = s, . . . , s +m − 1), and replace the Hessian matrices by the negative of information matrices to obtainlog f ({xt }s+m−1|hs , hs+m+1 , ys , . . . , ys+m )
t=s
s+m−1
1 X 0
≈ const. −xt xt + L̂
2 t=s
¸
s+m
X ·
1
0
0
0
ˆ
+
dt − (αt − α̂t ) Ât + (αt−1 − α̂t−1 ) B̂t (αt − α̂t )
2
t=s+1
∗

= const. + log f (x|αs , αs+m+1 , ys , . . . , ys+m ),(24)

where d0t = ∂L/∂α0t ,¸∂2L, t = s + 1, . . . , s + m,
= −E
∂αt ∂α0t
·
¸
∂2L
= −E
, t = s + 2, . . . , s + m,
∂αt ∂α0t−1
·

At
Bt

(25)
Bs+1 = O,(26)

and dˆt , Ât , B̂t are those evaluated at α̂t . The expectations are taken with respect to yt
given parameters and other latent variables.

A.1

Derivation of dt , At and Bt

Summary of matrix differentiation
We first summarize definitions for the first and second derivatives of a matrix and someresults (Magnus and Neudecker (1999), and Magnus and Abadir (2007)). Let F be a twice
differentiable m×p matrix function of an n×q matrix X. Then the first derivative (Jacobianmatrix) of F at X is defined by the mp × nq matrix
DF (X) =∂F (X)∂vec(F (X))=,∂X
∂vec(X)031

where vec(·) is a vectorizing operator, and the second derivative (Hessian matrix) of F at
X is defined by the mnpq × nq matrix¡¢
HF (X) = D (DF (X))0 =∂vec∂(vec(X))0µµ∂vec(F (X))∂(vec(X))0¶0 ¶.

Chain rule: Let S a subset of Rn×q , and assume that F : S → Rm×p is differentiable at aninterior point C of S. Let T be a subset of Rm×p such that F (X) ∈ T for all X ∈ S, andassume that G : T → Rr×s is differentiable at an interior point B = F (C) of T . Then thecomposite function H : S → Rr×s defined by H(X) = G(F (X)) is differentiable at C, andDH(X) = (DG(F (X)))(DF (X)) =∂vec(G(F (X))) ∂vec(F (X)).
∂(vec(F (X)))0 ∂(vec(X))0(27)

When q = 1, x ∈ Rn×1 , f : Rn×1 → Rm×p , g : Rm×p → Rr×s ,∂g(f (x))∂x0=∂vec(g(f (x))) ∂vec(f (x)).
∂vec(f (x))0 ∂vec(x)0(28)

Product rule: Let S a subset of Rn×q , and assume that F : S → Rm×p and G : S → Rp×r
are differentiable at an interior point C of S. Then∂vec(F )∂vec(G)∂vec(F G)= (G0 ⊗ Im )+ (Ir ⊗ F ).
∂(vec(X))0∂(vec(X))0∂(vec(X))0(29)

Derivation of dt
Let Ft = − 12 Ht and zt = exp(Ft )yt . The logarithm of the conditional probability densityof yt given ht excluding the constant term is
11
lt = − tr(Ht ) − (zt − mt )0 S−1
t (zt − mt ),2200
where E[zt ] = mt = Σεη Σ−1ηη (αt+1 − Φαt )I(t < n) and E[zt zt ] = St + mt mt with St =2
Ip −Σεη Σ−1ηη Σηε I(t < n). Further, let Dp denote a p ×q duplication matrix (whose elements

are 0 or 1) such that vec(A) = Dp vech(A) for a symmetric matrix A. Then∂tr(Ht )
∂α0t
∂mt
∂α0t
∂mt−1
∂α0t

=

∂tr(Ht ) ∂vec(Ht )= vec(Ip )0 Dp = vech(Ip )0 ,∂vec(Ht )0 ∂h0t(30)

= −Σεη Σ−1ηη ΦI(t < n),(31)= Σεη Σ−1ηη I(t > 1),

(32)

32

where we used the chain rule and D0p vec(A) = vech(A + A0 − (A ¯ Ip )) for a p × p matrix
A in (30) (e.g. Magnus and Neudecker (1999), Magnus (1988)). Further, define

Vt ≡

∞
i
i
∂vec(exp(Ft )) X 1 X h i−j
j−1
=
F
⊗
F
,
t
t
∂vec(Ft )0
i!
i=1

j=1

using the product rule and let Qt = {exp(−Ft )⊗Ip }Vt . Noting that zt = (yt0 ⊗Ip )vec(exp(Ft )),
∂zt
∂α0t

1
∂vec(Ht )1= − (yt0 ⊗ Ip )Vt= − (zt0 ⊗ Ip )Qt Dp .
0
2
∂ht
2

(33)

Using (30) – (33) and ∂x0 Ax/∂x0 = 2x0 A for a p × p symmetric matrix A and a p × 1
vector x, we obtain
·

dt

¸ ·
¸
∂lt 0
∂lt−1 0=++ ΦΣ−1ηη (αt+1 − Φαt )I(t = s + m < n)
∂α0t
∂α0t
½
¾
1 0 0
1
Dp Qt (zt ⊗ Ip ) − ΦΣ−1
Σ
I(t
<
n)
S−1
= − vech(Ip ) +ηη ηεt (zt − mt )
2
2
−1
−1
+Σ−1
ηη Σηε St−1 (zt−1 − mt−1 )I(t > 1) + ΦΣηη (αt+1 − Φαt )I(t = s + m < n). (34)

Although Qt involves an infinite series of matrices, its computation is easy as shown inAppendix A.2.
Derivation of At
By (31)–(34) and Q0t (zt ⊗ Ip ) = Vt (yt ⊗ Ip ),At·¸∂Vt (yt ⊗ Ip )S−1
1 0
1t (zt − mt )= − Dp E− ΦN0t02∂αt2+ΦMt Φ + Mt−1 + ΦΣ−1ηη ΦI(t = s + m < n),

(35)

where
−1
−1
Mt = Σ−1ηη Σηε St Σεη Σηη I(1 ≤ t < n),

(36)

−1
Nt = D0p Q0t (mt ⊗ Ip )S−1t Σεη Σηη .

(37)

33

Using a product rule and ∂vec(Ft )/∂α0t = − 12 Dp , we obtain∂Vt (yt ⊗ Ip )S−1
t (zt − mt )
∂α0t
ª© 0
ª
1©
∂(zt − mt )= − (zt − mt )0 S−1(yt ⊗ Ip ) ⊗ Ip2 Wt Dp + Vt (yt ⊗ Ip )S−1,
t ⊗ Ip2t2∂α0twhere Wt = ∂vec(Vt )/∂vec(Ft )0 . Since£¤0
E {(zt − mt )0 S−1t ⊗ Ip2 }{(yt ⊗ Ip ) ⊗ Ip2 }£¤0= E (zt − mt )0 S−1
t (yt ⊗ Ip ) ⊗ Ip2
£
¤
0 0
0
= E vec(S−1t (zt − mt )yt ) ⊗ Ip2 = vec(exp(−Ft )) ⊗ Ip2 ,and·E∂(zt − mt )Vt (yt ⊗ Ip )S−1
t
∂α0t

¸

£
¤
1
−1
0
0
−1
= − Q0t E (zt ⊗ Ip )S−1
t (zt ⊗ Ip ) Qt Dp + Qt (mt ⊗ Ip )St Σεη Σηη Φ
2
£
¤
1
−1
Qt Dp + Q0t (mt ⊗ Ip )S−1
= − Q0t E (zt zt0 ) ⊗ S−1tt Σεη Σηη Φ
2
ª
1 ©
−1
Qt Dp + Q0t (mt ⊗ Ip )S−1
= − Q0t (St + mt m0t ) ⊗ S−1tt Σεη Σηη Φ,
2Equation (35) reduces to
At =

©
ª ¤
1
1 0£
Dp Pt + Q0t (St + mt m0t ) ⊗ S−1
Qt Dp − (Nt Φ + ΦN0t )t
42+ΦMt Φ + Mt−1 + ΦΣ−1
ΦI(t
=
s
+
m
<
n),
ηη

(38)

©
ª
where D0p Pt = D0p vec(exp(−Ft ))0 ⊗ Ip2 Wt . The computation of Pt as well as Qt isdiscussed in Appendix A.2.
Derivation of Bt
By (31), (33) and (34),
·
Bt =

−1
−Σ−1
ηη Σηε St−1 E¸∂(zt−1 − mt−1 )1= N0t−1 − Mt−1 Φ,
0
∂αt−1
2

34

t = 2, . . . , n. (39)A.2

Computation of Pt and Qt(1) QtSince Ft is a symmetric matrix, there exists a p × p orthogonal matrix Ut such that
U0t Ft Ut = Λt ,where Λt = diag(λ1t , λ2t , . . . , λpt ) and λ1t ≥ λ2t ≥ . . . ≥ λpt are the ordered eigenvalues ofFt . Then
Qt

i
∞
i
X
1 X h i−j
= {exp(−Ft ) ⊗ Ip }Ft ⊗ Ftj−1 ,
i!
i=1

=

∞
X
i=1

1
i!

j=1

i h
X

i
j−1
exp(−Ft )Fi−j
⊗
F
,
t
t

j=1



i n
∞
o
X
X
1
exp(−Λt )Λi−j⊗ Λtj−1  (U0t ⊗ U0t ).
= (Ut ⊗ Ut ) 
t
i!
i=1

(40)

j=1

The second factor in (40) is a diagonal matrix with its (k, k)-th element given by
i
∞
X
1X
i=1

i!

j−1
exp(−λat )λi−j
at λbt =

j=1


 1,

if λat = λbt ,if λat 6= λbt ,exp(λbt −λat )−1,
λbt −λatwhere a = b(k − 1)/pc + 1 and b = k − pb(k − 1)/pc, and bxc denotes the integer part of x.
Note that λat = λbt for k = (i − 1)p + i (i = 1, . . . , p).
(2) PtLet Kmn denote a mn × mn 0-1 matrix called a commutation matrix such that vec(A0 ) =Kmn vec(A) holds for a m × n matrix A (see Chapter 3 of Magnus and Neudecker (1999)).
Using∂vec(A ⊗ B)
∂vec(A)0∂vec(A ⊗ B)
∂vec(B)0= (In ⊗ Kqm ⊗ Ip )(Imn ⊗ vec(B)),= (In ⊗ Kqm ⊗ Ip )(vec(A) ⊗ Ipq ),35

where A, B, X are m × n, p × q, n × p matrices respectively, we compute∞i∂vec(Vt ) X 1 X ∂vec(Fi−j⊗ Ftj−1 )t=∂vec(Ft )0i!
∂vec(Ft )0
i=2
j=1
( i−1
∞
X
X
1
∂vec(Fti−j )(Ip2 ⊗ vec(Ftj−1 ))=(Ip ⊗ Kpp ⊗ Ip )i!
∂vec(Ft )0

Wt =

i=2

j=1

)
i
j−1
X
∂vec(F)t+(vec(Fti−j ) ⊗ Ip2 ),∂vec(Ft )0
j=2

and, using the product rule,kX∂vec(Fkt )=Fk−m⊗ Fm−1.
tt∂vec(Ft )0
m=1

Note thatno ∂vec(Fi−j )
©
ª
j−1
0
t
vec(exp(−Ft )) ⊗ Ip2 (Ip ⊗ Kpp ⊗ Ip ) Ip2 ⊗ vec(Ft )∂vec(Ft )0
i−j= Kpp (Fj−1
exp(−Ft ) ⊗ Ip )t

X j−1
∂vec(Fi−jt )Ft exp(−Ft )Fti−j−h ⊗ Fh−1
= Kpp
,
t
0
∂vec(Ft )
h=1

and
n
o ∂vec(Fj−1 )©
ªtvec(exp(−Ft ))0 ⊗ Ip2 (Ip ⊗ Kpp ⊗ Ip ) vec(Fti−j ) ⊗ Ip2∂vec(Ft )0j−1
no ∂vec(Fj−1 ) Xt=Fi−jexp(−F)⊗I=
Fti−j exp(−Ft )Fi−1−h⊗ Fh−1,tpttt∂vec(Ft )0
h=1

by using
Knn (BA ⊗ In ) = (vec(A)0 ⊗ In2 )(In ⊗ Knn ⊗ In )(In2 ⊗ vec(B)),(41)(BA0 ⊗ In ) = (vec(A)0 ⊗ In2 )(In ⊗ Knn ⊗ In )(vec(B) ⊗ In2 ),

(42)

36

where A, B are m × n, p × q matrices (e.g. Theorem 4.4 of Rogers (1980) on p. 23). Thus,
using D0p Kpp = D0p ,©
ª
D0p Pt = D0p vec(exp(−Ft ))0 ⊗ Ip2 Wt
Ã i−1 i−j
!
j−1
∞
i X
X
XX
X
1
= D0p
exp(−Ft )Fi−h−1⊗ Fh−1+exp(−Ft )Fti−k−1 ⊗ Fk−1
t
t
t
i!
i=2

= 2D0pi−1 i−j
∞
X
1 XX
i=2

= 2D0pj=1 h=1i!

j=2 k=1

exp(−Ft )Fi−h−1⊗ Fh−1tt

j=1 h=1∞i−1
X
1X(i − j) exp(−Ft )Fi−j−1⊗ Ftj−1 ,
t
i!
i=2

j=1

and we obtain∞i−1 i−jX
1 XX(i − j) exp(−Ft )Fi−j−1⊗ Ftj−1ti!
i=2j=1 h=1


∞
i−1
n
o
X
X
1
= (Ut ⊗ Ut ) 2(i − j) exp(−Λt )Λti−j−1 ⊗ Λtj−1  (U0t ⊗ U0t ). (43)i!

Pt = 2i=2

j=1

The second factor in (43) is a diagonal matrix with its (k, k)-th element given by

2

i−1
∞
X
1X
i=2

i!

j−1
(i − j) exp(−λat )λi−j
at λbt =

j=1


 1,

if λat = λbt ,bt −λat )} 2{exp(λbt −λat )−1−(λ, if λat 6= λbt ,(λbt −λat )2

where a = b(k − 1)/pc + 1 and b = k − pb(k − 1)/pc.

References
Abadir, K. M. and J. R. Magnus (2005). Matrix algebra. New York: Cambridge University Press.
Andersen, T. G., T. Bollerslev, F. X. Diebold, and H. Ebens (2001). Thedistribution of stock
return volatility. Journal of Financial Economics 61, 43–76.
Asai, M. and M. McAleer (2006). Asymmetric multivariate stochastic volatility. Econometric
Reviews 25, 453–473.
Asai, M. and M. McAleer (2009). Multivariate stochastic volatility, leverage and news impact
surfaces. Econometrics Journal 12, 292–309.
Asai, M., M. McAleer, and J. Yu (2006). Multivariate stochastic volatility: A review. Econometric
Reviews 25, 145–175.

37

Audrino, F. and F. Corsi (2010). Modelling tick-by-tick realized correlations. Computational
Statistics and Data Analysis 54, 2372–2382.
Barndorff-Nielsen, O. E. and N. Shephard (2004). Econometric analysis of realized covariation:high frequency covariance, regression and correlation in financial economics. Econometrica 72,
885–925.
Bauer, G. H. and K. Vorkink (2010). Forecasting multivariate realized stock market volatility.
JoE . doi:10.1016/j.jeconom.2010.03.021.
Bauwens, L., S. Laurent, and J. V. K. Rombouts (2006). Multivariate GARCH: a survey. Journalof Applied Econometrics 21, 79–109.
Chan, D., R. Kohn, and C. Kirby (2006). Multivariate stochastic volatility models with correlated
errors. Econometric Reviews 25, 245–274.
Chib, S. (2001). Markov chain Monte Carlo methods: computation and inference. In J. J. Heckman
and E. Leamer (Eds.), Handbook of Econometrics, Volume 5, pp. 3569–3649. Amsterdam:
North-Holland.
Chib, S., Y. Omori, and M. Asai (2009). Multivariate stochastic volatility. In T. G. Andersen,
R. A. Davis, J. P. Kreiss, and T. Mikosch (Eds.), Handbook of Financial Time Series, pp.
365–400. New York: Springer-Verlag.
Chiriac, R. and V. Voev (2010). Modelling and forecasting multivariate realized volatility. Journalof Applied Econometrics. DOI: 10.1002/jae.1152.
Chiu, T., T. Leonard, and K. Tsui (1996). The matrix-logarithmic covariance model. Journal of
the American Statistical Association 91, 198–210.
Danı́elsson, J. (1998). Multivariate stochastic volatility models: Estimation and a comparison
with VGARCH models. Journal of Empirical Finance 5, 155–173.
de Jong, P. and N. Shephard (1995). The simulation smoother for time series models.
Biometrika 82, 339–350.
Doornik, J. (2006). Ox: Object Oriented Matrix Programming. London: Timberlake Consultants
Press.
Durbin, J. and S. J. Koopman (2002). A simple and efficient simulation smoother for state spacetime series analysis. Biometrika 89, 603–616.
Engle, R. F. (2002). Dynamic conditional correlation: A simple class of multivariate generalized autoregressive conditional heteroskedasticity models. Journal of Business and Economic
Statistics 20, 339–350.
Engle, R. F. and K. F. Kroner (1995). Multivariate simultaneous generalized ARCH. Econometric
Theory 11, 122–150.
Engle, R. F. and V. Ng (1993). Measuring and testing the impact of news on volatility. Journalof Finance 48, 1749–1778.

38

Gourieroux, C., J. Jasiak, and R. Sufana (2009). The Wishart autoregressive process of multivariate stochastic volatility. Journal of Econometrics.
Gupta, A. K. and D. K. Nagar (2000). Matrix variate Distributions. Boca Raton: Chapman &Hall/CRC.
Ishihara, T. and Y. Omori (2012). Efficient Bayesian estimation of a multivariate stochastic volatility model with cross leverage and heavy-tailed errors. Computational Statistics and Data
Analysis 56, 3674–3689.
Jin, X. and J. M. Maheu (2010). Modelling realized covariances and returns. Working Paper 408,
University of Toronto.
Jungbacker, B. and S. J. Koopman (2006). Monte Carlo likelihood estimation for three multivariate stochastic volatility models. Econometric Reviews 25, 385–408.
Kawakatsu, H. (2006). Matrix exponential GARCH. Journal of Econometrics 134, 95–128.
Koopman, S. J. (1993). Disturbance smoother for state space models. Biometrika 80, 117–126.
LeSage, J. P. and R. K. Pace (2007). A matrix exponential spatial specification. JoE 140, 190–214.
Lopes, H. F., R. E. McCulloch, and R. S. Tsay (2012). Cholesky stochastic volatility models forhigh-dimensional time series. Discussion papers.
Magnus, J. R. (1988). Linear Structures. Oxford: Oxford University press.
Magnus, J. R. and K. M. Abadir (2007). On some definitions in matrix algebra. Discussion paper,
CIRJE-F-476, Faculty of Economics, University of Tokyo.
Magnus, J. R. and H. Neudecker (1999). Matrix differential culculus with applications in statistics
and econometrics, revised edition. Chichester: John Wiley.
Omori, Y., S. Chib, N. Shephard, and J. Nakajima (2007). Stochastic volatility with leverage:
fast likelihood inference. Journal of Econometrics 140, 425–449.
Omori, Y. and T. Watanabe (2008). Block sampler and posterior mode estimation for asymmetricstochastic volatility models. Computational Statistics and Data Analysis 52-6, 2892–2910.
Philipov, A. and M. E. Glickman (2006). Multivariate stochastic volatility via Wishart processes.
Journal of Business and Economic Statistics 24, 313–328.
Pitt, M. K. and N. Shephard (1999). Filtering via simulation: Auxiliary particle filters. Journalof the American statistical association 94 (446), 590–599.
Rogers, G. S. (1980). Matrix Derivarives: Lecture Notes in Statistics, Vol. 2. New York: Marcel
Dekker.
Shephard, N. and M. K. Pitt (1997). Likelihood analysis of non-Gaussian measurement time series.
Biometrika 84, 653–667.
Sheppard, K. (2007). Positive semi-definite matrix multiplicative error models. University of Oxford.

39

Spiegelhalter, D. J., N. G. Best, B. P. Carlin, and A. van der Linde (2002). Bayesian measuresof model complexity and fit (with discussion). Journal of the Royal Statistical Society, Series
B 64, 583–639.
Tsay, R. S. (2005). Analysis of Financial Time Series: Financial Econometrics (2 ed.). New York:
Wiley.
Voev, V. (2007). Dynamic modelling of large-dimensional covariance matrices. In L. Bauwens,
W. Pohlmeier, and D. Veredas (Eds.), Recent Developments in High Frequency Financial
Econometrics, pp. 365–402. Berlin: Springer-Verlag.
Watanabe, T. and Y. Omori (2004). A multi-move sampler for estimating non-Gaussian times
series models: Comments on Shephard and Pitt (1997). Biometrika 91, 246–248.
Yu, J. (2005). On leverage in a stochastic volatility model. Journal of Econometrics 127, 165–178.
Yu, J. and R. Meyer (2006). Multivariate stochastic volatility models: Bayesian estimation andmodel comparison. Econometric Reviews 25, 361–384.

40