Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
90,334 characters · 21 sections · 31 citation commands
Multivariate Stochastic Volatility Model with Realized Volatilities and Pairwise Realized Correlations
Modelling the time-varying volatility and the correlations of multivariate time series is one of the most important problems in financial risk management, and there are numerous studies that model the time-varying volatility of univariate time series using the GARCH or stochastic volatility (SV) models. However, the extension of their models to multivariate model with dynamic correlations has not been straightforward due to the following several major problems.
First, there are too many parameters to estimate if the only available data are daily returns, which results in unstable estimates. An intuitive solution to reduce the number of parameters is to introduce the factor structure assuming that a small number of common factors describe the dynamics of time-varying covariance matrices as discussed in the factor stochastic volatility models (e.g. pitt1999time, chib2006analysis and lopes2007factor). However, factor modelling requires the a priori selection of the number of factors and we need to restrict the structure of the factors in order to identify the parameters (e.g. lopes2004bayesian). Furthermore, the estimation results and the predictive performance of the model are usually subject to the ordering of the asset returns.
An alternative effective approach is to incorporate additional observations based on the intraday asset returns, such as the realized covariances, which have recently become available in financial markets. In univariate SV models, the realized SV (RSV) models that estimate the time-varying volatilities using the daily returns and realized volatility simultaneously have been proposed to achieve more accurate parameter estimates than those of the SV models using only daily returns and the RSV models outperform SV models in forecasting volatilities (e.g. TakahashiOmoriWatanabe(09), DobrevSzerszen(10), KoopmanScharth(13), ZhengSong(14), TakahashiWatanabeOmori(16)). Although the realized volatilities are subject to microstructure noises and nontrading hours and hence are biased estimates of the integrated volatilities, such biases are automatically adjusted within the proposed model. Similarly, the univariate GARCH model is extended to the realized GARCH models which incorporates the realized volatilities into the variance equations and it is shown to lead to substantial improvements in the empirical fit and quantile forecasts over the standard GARCH model that only uses daily returns (HansenHuangShek(12)).
The extension to the multivariate RSV model is also considered in the Cholesky RSV model (ShirotaOmoriLopesPiao(17)). In this model, the Cholesky decompositions of the realized covariance matrices are used as additional sources for measurement equations, and it models the dynamics of the logarithm of the diagonal elements and the off-diagonal elements of Cholesky decomposed covariance matrices respectively. It is shown that the portfolio performances of the proposed model outperformed other SV models without realized measures in the empirical studies, but it should also be noted that the performance of the Cholesky RSV models may depend on the ordering of the asset returns in the vector of the response.
Second, high-frequency data are not always observed at the same time points, which causes difficulties in the extension of the univariate RSV model to the multivariate RSV model. For example, in the Cholesky RSV model, it is implicitly assumed that the all multivariate assets are traded every few minutes when computing the realized covariance matrices. If the multivariate assets are not traded synchronously, we have to use the largest time intervals so that all asset returns are observed when computing the realized covariance matrices. This nonsynchronous trading leads us to ignore some of the frequently traded asset return data, and hence we would fail to make full use of the available intraday informations when there are less frequently traded assets.
Third, it is not straightforward to guarantee that the estimated (and the realized) covariance matrices are positive definite. The model parameters may be difficult to estimate in practice under the constraints that satisfy the positive definiteness. Using the Cholesky decomposition of the time-varying covariance matrices is one way to guarantee the positive definiteness (ShirotaOmoriLopesPiao(17)), but it also requires that the multivariate assets are traded synchronously in order to compute the realized covariance matrices as mentioned above. Additionally, the interpretation of each latent variable of the decomposition is not straightforward since it does not correspond to each pair of asset returns and it is subject to the ordering of the asset returns. If we use each element of the realized covariances for each pair of asset returns, it may result in the nonpositive definite covariance matrices.
To overcome these difficulties, we propose a multivariate realized SV (MRSV) model with pairwise realized correlations, in which we incorporate the dynamic latent correlation variables in addition to latent volatility variables with realized measures for each pairwise correlation and the volatilities in the framework of multivariate SV models with realized volatilities. The model parameters are estimated using Markov chain Monte Carlo simulations, and we sample the latent correlation variables one at a time given the others so that we keep the covariance matrices positive definite. The realized Beta GARCH model proposed by HansenLundeVoev(14) is a promising multivariate GARCH model with realized measures for volatilities and co-volatilities in which they used measurement equations for the pairwise realized correlations with market returns and modelled dynamics of the Fisher transformed conditional correlation coefficients. However, they focused on the pairwise correlations between the market return and an individual asset return, assuming that the individual asset returns are conditionally independent given the market return. Another useful approach for the joint modelling of returns and realized covariances are based on Wishart processes (e.g. JinMaheu(13), WindleCarvalho(14), JinMaheu(16), SoLiAsaiJiang(16)). The covariance matrix is assumed to follow a Wishart distribution whose scale matrix depends on the past realized covariance matrices which are computed using larger time intervals than necessary in order for it to be positive definite.
Our approach, on the other hand, is based on simultaneously modeling the individual volatilities and pairwise covariances, rather than the covariance matrix, and we are able to make full use of the available intraday information, even when there are less frequently traded assets. This finding implies that our model still can be constructed even though some of the realized measures are missing. Additionally, our model is far more flexible in the sense that it is possible to restrict any correlation coefficients to be zero for very high dimensional asset returns data, which reduces the number of parameters and may improve the forecasting performances. Among the multivariate SV models in the literature, our model is a natural extension of a univariate RSV model to a multivariate model, and it gives us a straightforward interpretation of the estimated parameters.
Furthermore, we extend our model to incorporate the leverage effect, which is well-known to exist in stock markets. The leverage effect refers to the negative correlation between an asset return and its volatility. In other words, a decrease in the stock return is followed by an increase in its volatility. In forecasting the means and covariances of asset returns for portfolio optimization, it is expected that incorporating the leverage effect in econometric models improves the predictive accuracy. However, it may increase the number of parameters that are to be estimated and the realized measures are not available for such an effect. Thus we also consider the parsimonious parameterization for the leverage effect.
Our contributions are as follows: (1) we obtain the stable parameter estimates for the dynamic correlation models using the realized measures, (2) we make full use of the intraday informations by using pairwise realized correlations, (3) the estimated covariance matrices are guaranteed to be positive definite, (4) we avoid the arbitrariness of the ordering of asset returns, (5) we propose the flexible correlation structure model and set some correlations to be zero if necessary, and (6) we introduce the parsimonious specification for the leverage effect.
The structure of this paper is as follows. Section 2 introduces the multivariate realized SV model with daily returns, realized volatilities, and pairwise realized correlations. Section 3 describes the estimation algorithms using the Markov Chain Monte Carlo simulation. Section 4 extends it to incorporate the leverage effect. Finally, in Section 5, the proposed model is applied to nine U.S. stock return data and the model with the leverage effect is shown to outperform other competing models with regard to the portfolio performances.
This section introduces the multivariate realized stochastic volatility (MRSV) model, which uses realized measures for the volatility and pairwise correlations of asset returns. By using the additional information of the realized measure for asset returns, we can overcome the curse of dimensionality when estimating the dynamic covariance matrices. The Cholesky RSV model proposed by ShirotaOmoriLopesPiao(17) also uses the realized measure of variances and covariances (which we call the realized covariance matrix) in order to estimate the latent covariance matrix of asset returns. However, the realized covariance matrix is less informative when there are less frequent asset returns. This finding is observed because we require the synchronous observations of all asset return series in order to compute the realized covariance matrix. In order to utilize the full information of the realized measure for the correlations, we propose using the realized measures for the latent pairwise correlations. It should be noted that the pairwise correlation can be computed if that pair of series is synchronously observed. We call the realized measure for the correlation coefficient the pairwise realized correlation. Using the pairwise realized correlations, in order to guarantee the positive definiteness of the latent covariance matrices, we propose the MCMC algorithm in which we sample latent correlation coefficients from the conditional posterior distribution so that the matrices are positive definite.
First, we define the multivariate SV (MSV) model without the realized measures. Let $\operatorname{\mbox{\boldmath $y$}}_t= (y_{1t},\ldots,y_{pt})'$ and $\operatorname{\mbox{\boldmath $h$}}_t = (h_{1t},\ldots,h_{pt})'$ denote a $p \times 1$ stock return vector and its corresponding log volatility latent vector at time $t$. The basic MSV model is given by
where $\operatorname{\mathbf{R}}_t=\{\rho_{ij,t}\}$ is a correlation matrix, $\operatorname{\mbox{\boldmath $\epsilon$}}_t = (\epsilon_{1t},\ldots,\epsilon_{pt})'$, $\operatorname{\mbox{\boldmath $\eta$}}_t = (\eta_{1t},\ldots,\eta_{pt})'$, and
We assume that $h_{it}$ follows a stationary autoregressive process (with its coefficient $\lvert\phi_j\rvert < 1)$ and that the mean process $\bm{m}_t= (m_{1t},\ldots,m_{pt})'$ follows a random walk process. We denote a diagonal matrix $\mathbf{A}$ with diagonal elements $\bm{a}=(a_{11},\ldots,a_{mm})'$ as $\mathbf{A}=\mbox{diag}(\bm{a})$. For the initial distributions of $\bm{m}_1$ and $\bm{h}_1$, we set $\kappa$ to some large constant for $\bm{m}_1$ for simplicity and set $\operatorname{\mathbf{\Omega}}_0$ to satisfy the stationary condition $\operatorname{\mathbf{\Omega}}_0 = \operatorname{\mathbf{\Phi}} \operatorname{\mathbf{\Omega}}_0 \operatorname{\mathbf{\Phi}} + \operatorname{\mathbf{\Omega}}$ for $\bm{h}_1$ such that
where $\mathbf{I}_{p^2}$ denotes a $p^2\times p^2$ unit matrix. In order to model the dynamics of the correlation matrix, we consider the following Fisher transformation $g_{ij,t+1}$ of the correlation coefficient $\rho_{ij,t}$, and assume that it follows a random walk process for simplicity:
for $i,j=1,\ldots, p$ $(j<i)$ and we denote $\bm{\rho}_t = (\rho_{21,t},\ldots,\rho_{p\hspace{0.1mm}p-1, t})'$, $\bm{g}_t = (g_{21,t},\ldots,g_{p\hspace{0.1mm}p-1, t})'$, $\bm{\zeta}_t = (\zeta_{21,t},\ldots,\zeta_{p\hspace{0.1mm}p-1, t})'$, and $\bm{\sigma}^2_{\zeta} = (\sigma^2_{\zeta,21},\ldots,\sigma^2_{\zeta, p\hspace{0.1mm}p-1})'$. \\ {\it Non-arbitrary ordering of asset returns and the flexible correlation structure.} We note that above specifications ((ref)) -- ((ref)) are independent of the ordering of the asset returns in $\bm{y}_t$, while the conventional factor SV models or the Cholesky SV models (ShirotaOmoriLopesPiao(17)) may be affected by the ordering. Further, it allows us to model the structure of the correlations in a flexible way. For example, we can easily restrict some correlation coefficients to be zero when the dimension of $\bm{y}_t$ is very high.\\
{\it Remark} 1. It is easy to assume that $\bm{m}_t$ and $g_{ij,t}$ follow stationary autoregressive processes. However, since it imposes the mean reversion properties on these processes, we would rather consider random walk processes without such properties for simplicity. For the long term prediction, we may need such a stationarity condition.
{\it Realized measures as an additional source of information}. In the above MSV models, there are too many parameters to estimate using only daily asset returns, and the parameter estimates are often unstable. Recently, high frequency data in the financial markets have become available, and they play a more important role in the finance-related empirical studies, since the realized measures of the variances and covariances, are more informative estimators of the true variances and covariances (see e.g. AndersenBollerslevDieboldLabys(01), AndersenBollerslevDieboldEbens(01), BarndorffShephard(02), BarndorffShephard(04)).
Let $x_{it}= \log RV_{it}$ and $w_{ij,t} = \log \{(1 + RCOR_{ij,t})/(1 - RCOR_{ij,t})\}$ where $RV_{it}$ and $RCOR_{ij,t}$ are the realized measures of the volatility of the $i$-th asset return and the correlation between $i$-th and $j$-th asset returns at time $t$. Thus we introduce the following additional measurement equations based on the realized measures:
for $i,j=1,\ldots, p$ $(i>j)$. The terms $\xi_{j}$ and $\delta_{ij}$ are included in order to adjust the biases due to the microstructure noise, nontrading hours, nonsynchronous trading and so forth. The multivariate realized stochastic volatility model with pairwise realized correlations is defined by ((ref)) -- ((ref)). We denote $\bm{x}_t = (x_{1t},\ldots,x_{pt})', \quad \bm{w}_t = (w_{21,t},\ldots,w_{p\hspace{0.1mm}p-1,t})'$, $\operatorname{\mbox{\boldmath $\xi$}} = (\xi_1,\ldots,\xi_p)'$, $\operatorname{\mbox{\boldmath $\delta$}} = (\delta_{21},\ldots,\delta_{p\hspace{0.1mm}p-1})'$, $\bm{u}_t = (u_{1t},\ldots,u_{pt})'$, $\bm{v}_t = (v_{21,t},\ldots,v_{p\hspace{0.1mm}p-1,t})'$, $\bm{\sigma}^2_u = (\sigma^2_{u,1},\ldots,\sigma^2_{u,p})'$, and $\bm{\sigma}^2_v = (\sigma^2_{v,21},\ldots,\sigma^2_{v,p\hspace{0.1mm}p-1})'$. \\ {\it Use of pairwise realized correlations.} Given the realized correlation $RCOR_{ij,t}$, we will use the pairwise realized correlations. If there is less frequent series of asset returns, the realized covariance matrix may lose a large part of the information since it is calculated only when all the series are synchronously observed. On the other hand, the pairwise realized correlation coefficients can be respectively calculated for each pair of series of returns; therefore, we can use the full information of the realized measures for the correlations. Moreover, we can estimate the parameters even if we cannot obtain the realized measures for some pairs. \\ {\it Bias corrections of the realized measures}. The realized volatilities and pairwise realized correlations have more information about the true volatilities and correlations, but there may be biases due to the market microstructure noise, nontrading hours, nonsynchronous trading and so forth. In order to correct these biases in the realized measures, we model the observation equations of the realized volatilities and pairwise realized correlations with bias adjustment terms, $\xi_j$ and $\delta_{ij}$. Although daily returns have relatively less information about the true volatilities and correlations, they are less subject to the biases that are caused by the high frequency data. Therefore, we can estimate the biases in the realized measures using the information of daily returns and also get additional information with regard to the true volatilities and correlations using the realized measures.
Since there are many latent variables in our proposed model and hence it is difficult to evaluate the likelihood, we take the Bayesian approach and estimate the model parameters using the Markov chain Monte Carlo simulation. First we assume the prior distribution of $\operatorname{\mbox{\boldmath $\theta$}}\equiv (\bm{\phi},\operatorname{\mbox{\boldmath $\mu$}},\operatorname{\mbox{\boldmath $\xi$}},\operatorname{\mbox{\boldmath $\delta$}}, \operatorname{\mbox{\boldmath $\sigma$}}_{u}^2,\operatorname{\mbox{\boldmath $\sigma$}}_{v}^2,\operatorname{\mbox{\boldmath $\sigma$}}_{\zeta}^2,\mathbf{\Sigma}_m,\operatorname{\mathbf{\Omega}})$ as follows. For the prior distributions of $\mu_i, \xi_i$ and $\delta_{ij}$, we assume multivariate independent normal distributions. The prior distributions of $\sigma_{u,i}^2,\sigma_{v,ij}^2,\sigma_{\zeta,ij}^2$ and $\sigma_{m,i}^2$ are assumed to be independent inverse gamma distributions. For $\phi_i$ and $\operatorname{\mathbf{\Omega}}$, we assume $(1+\phi_i)/2 \sim \operatorname{\operatorname{Beta}}(a,b)$ and an inverse Wishart distribution respectively. In summary, we assume the following prior distributions:
for $i, j=1,\ldots,p$ $(j<i)$, and $a, b,m_\mu,s_{\mu}, m_{\xi},s_{\xi}, m_{\delta},s_{\delta}, n_u, d_u, n_v, d_v,n_{\zeta}, d_{\zeta}, n_m, d_m, \nu,\mbox{\boldmath $S$}$ are hyperparameters. \\ {\it Remark} 2. The particle MCMC may be a possible alternative estimation method to the MCMC below for the univariate models, but it may not be appropriate for the multivariate models since the discrete approximation to the high dimensional state distribution often results in the degeneracy of the particles.
Let $\bm{g}=(\bm{g}_1',\ldots,\bm{g}_T')'$, $\bm{h}=(\bm{h}_1',\ldots,\bm{h}_T')'$ and $\bm{m}=(\bm{m}_1',\ldots,\bm{m}_T')'$. Further, let $\bm{w}=(\bm{w}_1',\ldots,\bm{w}_T')'$, $\bm{x}=(\bm{x}_1',\ldots,\bm{x}_T')'$ and $\bm{y}=(\bm{y}_1',\ldots,\bm{y}_T')'$. In order to conduct the statistical analysis of the parameters, we implement the Markov chain Monte Carlo simulation in nine blocks. The MCMC sampling algorithm is described in more details in the following subsections. Let $\operatorname{\mbox{\boldmath $\theta$}}_{\backslash \operatorname{\mbox{\boldmath $\beta$}}}$ denote the parameter $\operatorname{\mbox{\boldmath $\theta$}}$ excluding $\operatorname{\mbox{\boldmath $\beta$}}$. Then,
The conditional posterior probability density function of $g_{ij,t}$ given other parameters and latent variables is
where
and
{\it Positive definiteness of $\mathbf{R}_t$}. We use an identity matrix for the initial value of $\mathbf{R}_t$ when implementing the MCMC. Thus, given the current correlation matrix $\mathbf{R}_t$, we generate each correlation coefficient $\rho_{ij,t}$ (or equivalently $g_{ij,t})$ so that we guarantee that the proposed $\mathbf{R}_t^*$ is the correlation matrix. We first state the condition for $\rho_{ij,t}$ to guarantee that the proposed $\mathbf{R}_t^{*}$ is positive definite given the other elements of $\operatorname{\mathbf{R}}_t$ and other $\rho_{ij,s}$ $(s\neq t)$ .\\
{\bf Proposition 1}. Suppose that $\operatorname{\mathbf{R}}_t=\{\rho_{ij,t}\}$ is a correlation matrix and let $\operatorname{\mbox{\boldmath $\rho$}}_{it}$ denote the transpose of the $i$-th row vector of $\operatorname{\mathbf{R}}_t$ excluding 1, $\operatorname{\mbox{\boldmath $\rho$}}_{it} = (\rho_{i1,t},\ldots,\rho_{i\hspace{0.1mm}i-1,t},\rho_{i\hspace{0.1mm}i+1,t},\ldots,\rho_{ip,t})'$, and $\operatorname{\mathbf{R}}_{it}$ denotes the submatrix excluding the $i$-th row and the $i$-th column from $\operatorname{\mathbf{R}}_t$. The condition for $\rho_{ij,t}$ to guarantee that $\mathbf{R}_t$ is positive definite is $\rho_{ij,t} \in (L_{ijt}, U_{ijt})$ where bounds $L_{ijt}$ and $U_{ijt}$ are given by
and $\operatorname{\mbox{\boldmath $\rho$}}_{i,-j,t}$ is the vector excluding the $j$-th element of $\operatorname{\mbox{\boldmath $\rho$}}_{it}$, $a_j$ is the $(j,j)$-th element of $\operatorname{\mathbf{R}}_{it}^{-1}$, $\mbox{\boldmath $b$}_j$ is the vector excluding $a_j$ from the $j$-th column of $\operatorname{\mathbf{R}}_{it}^{-1}$, and $\mathbf{C}_j$ is the matrix excluding the $j$-th row and $j$-th column from $\mathbf{R}_{it}^{-1}$. \\ {\bf Proof:} See Appendix (ref) \\
Thus we propose a candidate $g_{ij,t}^{\dagger}$ from normal distribution truncated on the interval $(a_{ijt}, b_{ijt}) $, $TN_{(a_{ijt},b_{ijt})}(m_{t*},\sigma_{t*}^2)$, and accept it with probability $\min\{1, \exp(r(g_{ij,t}^{\dagger})-r(g_{ij,t}))\}$, where
We use a single-move sampler for $\bm{h}_t$ in which we sample $\bm{h}_t$ given the other parameters and latent variables. Such a sampler is efficient when the realized measures are available as the additional information source for $\bm{h}_t$. The conditional posterior probability density function of $\operatorname{\mbox{\boldmath $h$}}_t$ is given by
where
where $\bm{1}_p$ denotes a $p\times 1$ vector with all elements equal to one. Therefore, we generate a candidate $\bm{h}_t^{\dagger}$ from $\operatorname{\operatorname{N}}(\bm{m}_{t*}, \operatorname{\mathbf{\Omega}}_{t*})$, and accept it with probability $\min\{1,\exp(l(\bm{h}_t^{\dagger}) - l(\bm{h}_t))\}$. See Appendix (ref) for the generations of $\mbox{\boldmath $\theta$}$ and $\operatorname{\mbox{\boldmath $m$}}_t$.
This section extends our model in order to incorporate the leverage effect. The leverage effect, which corresponds to the well-known negative correlation between asset returns and their volatilities in the stock market, is expected to improve the performance of the forecast of the mean processes and volatility processes of asset returns.
We first define the matrix variate normal distribution and show its probability density function, which will be used in modelling the leverage effect. \\
{\bf Definition 1}. The random matrix $\operatorname{\mathbf{X}}$ ($p \times n$) is said to have a matrix variate normal distribution with mean matrix $\operatorname{\mathbf{M}}$ ($p \times n$) and covariance matrix $\operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Sigma}}$ where $\operatorname{\mathbf{\Psi}}$ $(p \times p)$ and $\operatorname{\mathbf{\Sigma}}$ $(n \times n)$ are positive definite matrices if $\operatorname{vec}{(\operatorname{\mathbf{X}}')} \sim \operatorname{\operatorname{N}}(\operatorname{vec}{(\operatorname{\mathbf{M}}')}, \operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Sigma}})$ and we denote $\operatorname{\mathbf{X}} \sim \operatorname{\operatorname{N}}_{p,n}(\operatorname{\mathbf{M}}, \operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Sigma}})$.
We extend our proposed model to incorporate the leverage effect as follows. The joint distribution of $(\operatorname{\mbox{\boldmath $y$}}_t, \operatorname{\mbox{\boldmath $h$}}_{t+1})$ is given by
The marginal distributions of $\bm{y}_t$ and $\bm{h}_{t+1}$ given $\bm{h}_t$ are the same as before with $\mathbf{\Omega}=\mathbf{\Psi}+\operatorname{\mathbf{\Lambda}}\operatorname{\mathbf{\Lambda}}'$, but we note that
If $\mathbf{\Lambda}=\mathbf{O}$, it reduces to the model without leverage effect. The matrix $\mathbf{\Lambda}$ is the coefficient of the leverage for $\bm{z}_t = \mathbf{R}_t^{-1/2}\mathbf{V}_t^{-1/2}(\bm{y}_t-\bm{m}_t)$. We assume that the prior distribution of $\mathbf{\Lambda}$ given $\mathbf{\Psi}$ is $\operatorname{\operatorname{N}}_{p,p}(\operatorname{\mathbf{M}}_0, \operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Gamma}}_0)$. That is, $\mathbf{\Lambda} |\mathbf{\Psi} \sim \operatorname{\operatorname{N}}_{p,p}(\operatorname{\mathbf{M}}_0, \operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Gamma}}_0).$ \\
{\it Remark} 3. There are several ways to choose $\operatorname{\mathbf{R}}_t^{1/2}$. For example, we can use the spectral decomposition of the correlation matrix $\operatorname{\mathbf{R}}_t=\operatorname{\mathbf{P}}_t\operatorname{\mathbf{Q}}_t\operatorname{\mathbf{P}}_t'$ and $\operatorname{\mathbf{R}}_t^{1/2} = \operatorname{\mathbf{P}}_t \operatorname{\mathbf{Q}}_t^{1/2}$ where the $i$-th diagonal element of the diagonal matrix $\operatorname{\mathbf{Q}}_t$ is the $i$-th largest eigenvalue of $\operatorname{\mathbf{R}}_t$ and the $i$-th column of $\operatorname{\mathbf{P}}_t$ is the corresponding $i$-th eigenvector (and we set the first elements of the eigenvectors to be positive for the identification purpose). Thus the $i$-th element of $\bm{z}_t$ can be interpreted as the $i$-th market factor among $p$ asset returns. Alternatively, Cholesky decomposition, $\operatorname{\mathbf{R}}_t = \operatorname{\mathbf{R}}_t^{1/2}\operatorname{\mathbf{R}}_t^{1/2'}$, can be used so that $\operatorname{\mathbf{R}}_t^{1/2}$ is a lower triangular matrix where all the diagonal elements are equal to one, but we note that it is affected by the ordering of the asset returns.
The conditional posterior distribution of $\mathbf{\Lambda}$ is derived in the following Proposition and we generate $\operatorname{vec}(\operatorname{\mathbf{\Lambda}}) \vert \cdot \sim \operatorname{\operatorname{N}}(\operatorname{vec}(\operatorname{\mathbf{M}}_1'), \operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Gamma}}_1)$.\\
{\bf Proposition 2}. Suppose that the prior distribution of $\mathbf{\Lambda}$ given $\mathbf{\Psi}$ is $\operatorname{\operatorname{N}}_{p,p}(\operatorname{\mathbf{M}}_0, \operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Gamma}}_0)$. Then the conditional posterior distribution of $\mathbf{\Lambda}$ given other parameters and latent variables is $\operatorname{\operatorname{N}}_{p,p}(\operatorname{\mathbf{M}}_1, \operatorname{\mathbf{\Psi}} \otimes \operatorname{\mathbf{\Gamma}}_1)$ where
and $\bm{z}_t = \mathbf{R}_t^{-1/2}\mathbf{V}_t^{-1/2}(\bm{y}_t-\bm{m}_t)$ and $\bm{\eta}_t=\bm{h}_{t+1}-\bm{\mu}-\mathbf{\Phi}(\bm{h}_t-\bm{\mu})$. \\ {\bf Proof:} See Appendix (ref). \\
For the generations of other parameters and latent variables, see Appendix (ref) .
This subsection proposes the parsimonious specification for $\operatorname{\mathbf{\Lambda}}= [\operatorname{\mbox{\boldmath $\lambda$}}_1, \cdots, \bm{\lambda}_p]$, in order to reduce the number of leverage parameters from $p^2$ to $pq$ $(q\ll p)$ by setting $\operatorname{\mathbf{\Lambda}} = [\operatorname{\mbox{\boldmath $\lambda$}}_1, \cdots, \bm{\lambda}_q,\operatorname{\mbox{\boldmath $0$}},\ldots,\operatorname{\mbox{\boldmath $0$}}]$ since we do not have additional measurement equations for the leverage effect. Using the spectral decomposition to compute $\mathbf{R}_t^{1/2}$, we can interpret that the $i$-th column corresponds to the $i$-th market factor among asset returns $(i=1,\ldots,q)$. The number of factors, $q$, is expected to be small, e.g., $q=1$ or $q=2$.
The following proposition and the corollary shows the conditional posterior distribution of the parameters for the leverage effect under parsimonious specifications. \\
{\bf Proposition 3}. Let $\operatorname{\mathbf{\Lambda}} = [\operatorname{\mbox{\boldmath $\lambda$}}_1, \cdots, \bm{\lambda}_q,\operatorname{\mbox{\boldmath $0$}},\ldots,\operatorname{\mbox{\boldmath $0$}}]$ and $\bm{\lambda}=(\bm{\lambda}_1',\ldots,\bm{\lambda}_q')'$. If the prior distribution of $\bm{\lambda}$ is assumed to be normal, $\operatorname{\mbox{\boldmath $\lambda$}} \sim \operatorname{\operatorname{N}}(\operatorname{\mbox{\boldmath $m$}}_{0}, \operatorname{\mathbf{\Gamma}}_{0})$, then the conditional posterior distribution of $\operatorname{\mbox{\boldmath $\lambda$}}$ is $\operatorname{\mbox{\boldmath $\lambda$}} \vert \cdot \sim \operatorname{\operatorname{N}}(\operatorname{\mbox{\boldmath $m$}}_{1},\operatorname{\mathbf{\Gamma}}_{1})$ where
$\mathbf{A},\mathbf{B}$ are defined in ((ref)), $\operatorname{\mathbf{A}}_{1:q,1:q}$ denotes the first $q$ rows and the $q$ columns of $\operatorname{\mathbf{A}}$, $\text{vec}(\mathbf{X})\equiv (\bm{x}_1', \ldots, \bm{x}_m')'$ denotes a vectorization of the matrix $\mathbf{X}=\{\bm{x}_1, \ldots, \bm{x}_m\}$, and $\otimes$ denotes Kronecker product. \\ {\bf Proof:} See Appendix (ref). \\
{\bf Corollary 1}. Let $q=1$ and $\operatorname{\mathbf{\Lambda}} = [\operatorname{\mbox{\boldmath $\lambda$}}, \operatorname{\mbox{\boldmath $0$}},\ldots,\operatorname{\mbox{\boldmath $0$}}]$. If the prior distribution of $\bm{\lambda}$ is assumed to be normal, $\operatorname{\mbox{\boldmath $\lambda$}} \sim \operatorname{\operatorname{N}}(\operatorname{\mbox{\boldmath $m$}}_{0}, \operatorname{\mathbf{\Gamma}}_{0})$, then the conditional posterior distribution of $\operatorname{\mbox{\boldmath $\lambda$}}$ is $\operatorname{\mbox{\boldmath $\lambda$}} \vert \cdot \sim \operatorname{\operatorname{N}}(\operatorname{\mbox{\boldmath $m$}}_{1},\operatorname{\mathbf{\Gamma}}_{1})$ where
and $z_{1t}$ is the first element of $\bm{z}_t = \mathbf{R}_t^{-1/2}\mathbf{V}_t^{-1/2}(\bm{y}_t-\bm{m}_t)$.
See Appendix (ref).
This section applies our proposed model to the daily returns of nine U.S. stocks ($p=9$) with realized volatilities and pairwise realized correlations. The nine series of stock returns are JP Morgan (JPM), International Business Machine (IBM), Microsoft (MSFT), Exxon Mobil (XOM), Alcoa (AA), American Express (AXP), Du Pont (DD), General Electric (GE), and Coca Cola (KO). The sample period is from February 1, 2001 to December 31, 2009, and the number of observation is $T=2242$. The daily returns for the $i$-th stocks are defined as $y_{it} = 100 \times (\log p_{it} - \log p_{i, t-1})$, where $p_{it}$ is the closing price of the $i$-th asset at time $t$.
Time series plots of $y_{it}$ are shown in Figure (ref), and show that there is a very high volatility period in 2008 (the financial crisis when Lehman Brothers filed for Chapter 11 bankruptcy protection). Additionally, there are other relatively high volatility periods in 2001 (the dot-com bubble and the September 11 attacks) and in 2002 (the market turmoil during which Worldcom filed for Chapter 11 bankruptcy protection). The realized volatilities and pairwise realized correlations are computed from the realized covariance matrices for these assets which can be downloaded from the Oxford Man Institute website (see, Section 5 of noureldin2012multivariate for details). The prior distributions are assumed to be vague and flat in order to reflect the fact that we have little information with regard to the parameters:
for $i=1,\ldots,p$, $j= 1,\ldots,i-1$. The proposed model is estimated and we use the parsimonious specification of the leverage effect with the number of factors $q=1$.
We run 12,000 MCMC iterations and the first 2,000 iterations are discarded as the burn-in period. Table (ref) shows the posterior means, 95% credible intervals and inefficiency factors\footnote{The inefficiency factor is defined as $1+2\sum_{g=1}^\infty \rho(g)$, where $\rho(g)$ is the sample autocorrelation at lag $g$. This is interpreted as the ratio of the numerical variance of the posterior mean from the chain to the variance of the posterior mean from hypothetical uncorrelated draws. The smaller the inefficiency factor becomes, the closer the MCMC sampling is to the uncorrelated sampling.} for $\operatorname{\mbox{\boldmath $\mu$}}, \operatorname{\mbox{\boldmath $\xi$}}, \operatorname{\mbox{\boldmath $\phi$}},\operatorname{\mbox{\boldmath $\sigma$}}_u,\operatorname{\mbox{\boldmath $\sigma$}}_m$ and $\bm{\lambda}$. The inefficiency factors are relatively small (less than 130 ) in the multivariate stochastic volatility models and our algorithm works well. The posterior means and posterior standard deviations for $\bm{\delta}$, $\bm{\sigma}_v$ and $\bm{\sigma}_{\zeta}$ are shown in Tables (ref), (ref) and (ref), respectively. \\ {\it Mean processes and volatilities}. The posterior means of $\sigma_{m,i}$ are around $0.067\sim 0.098$, which reflects that the magnitude of the mean process, $\bm{m}_t$, is much smaller than that of the stochastic volatility component, $\mathbf{V}_t^{1/2}\bm{\epsilon}_t$, as we expected. The unconditional means of the log volatilities, $\mu_i$, are estimated to be from $0.203$ to $1.582$ and the posterior mean of $\mu_5$ (corresponding to Alcoa) is much larger than those of others. The stock returns of Alcoa are found to be the most volatile among others, while those of Coca Cola are the least volatile.
Since all posterior means of the autoregressive coefficients, $\phi_i$, are approximately 0.9, the log volatilities are found to have high persistence. The elements of $\mathbf{\Psi}$ (the conditional covariance matrix of $\bm{h}_{t+1}$ given $\bm{y}_t$) are all approximately 0.1 and the probability that $\psi_{ij}$ is positive is greater than 0.975 for all $i$s and $j$s. The log volatilities, $h_{i,t+1}$, are positively correlated with each other given $\bm{y}_t$. Figure (ref) shows the 95% credible intervals for $h_{1t}$ with $x_{1t}-\xi_1$ where $\xi_1$ is the estimated posterior mean of the first bias correction term. The figures for $h_{it}$ $(i=2,\ldots,9)$ are similar and hence are omitted. The estimated 95 % credible intervals have smaller fluctuation than those of the bias-adjusted realized measures. These estimates succeeded at automatically extracting the mean trends of the volatilities and adjusting the measurement errors. Overall, the 95% credible intervals captures the traceplot of the (bias-corrected) realized volatilities, suggesting that our proposed model is successful at describing the dynamics of the latent log volatilities.
{\it Biases in realized volatilities and correlations}. The bias correction terms, $\xi_i$, of the realized volatilities are estimated to be negative, thereby indicating that the realized volatilities have downward biases and underestimate the volatilities by ignoring the overnight nontrading hours. Since the realized volatilities tend to overestimate the volatilities due to the microstructure noises, the effect of nontrading hours seems to dominate in the direction of the biases. We also note that the magnitudes of the biases depend on the series of stock returns. Table (ref) shows the estimation result of the bias term $\operatorname{\mbox{\boldmath $\delta$}}$ of the correlation coefficients. All $\delta_{ij}$ are estimated to be negative, and the posterior probability that $\delta_{ij}$ is negative is greater than 0.975. This implies that the realized correlations underestimate the latent correlations, thereby suggesting the existence of the Epps effect. \\ {\it Dynamic correlations}. The posterior means of the standard deviations of the disturbance terms in the state equations corresponding to the dynamic correlations, $\sigma_{\zeta,ij}$, are shown in Table (ref). They are $0.039\sim 0.071$ and are much smaller than the posterior means of the standard deviations for the measurement errors of the realized measures, $\sigma_{u,i}$ and $\sigma_{v,ij}$ (as shown in Tables (ref) and (ref)) which are found to be similar for all $i$s and $j$s at approximately 0.30. Figure (ref) shows the time series plots of the 95% credible intervals of the selected dynamic correlations, $\rho_{21,t}$ with $\{\exp(w_{21,t}-\delta_{21})-1\}/\{\exp(w_{21,t}-\delta_{21})+1\}$ where $\delta_{21}$ is the estimated posterior mean. The figures for the other $\rho_{ij,t}$ are similar and hence are omitted. Again, the estimated 95% credible intervals of $\rho_{ij,t}$ have much smaller fluctuation those of bias-adjusted realized measures, $\{\exp(x_{ij,t}-\delta_{ij})-1\}/\{\exp(x_{ij,t}-\delta_{ij})+1\}$. These intervals seem to extract the mean trends of the bias-adjusted realized measures that have relatively large noises in the measurement equation. The correlations between the asset returns are found to be time-varying in the sample period, and they seem to increase after the financial crisis in 2008. This result corresponds to our intuition that each asset return has a larger positive correlation with others when the market faces stress, rather than when it is in a usual period. \\ {\it Leverage effect and the selection of the number of factors $q$}. The parameters for the leverage effect, $\lambda_i$, are estimated to be negative in Table (ref) and the posterior probability that $\lambda_i$ is negative is greater than 0.975 for all $i$s. This implies the existence of the leverage effect. Table (ref) also shows the estimation results for the correlation between the first element of $\bm{z}_t=\mathbf{V}_t^{-1/2}\mathbf{R}_t^{-1/2}\bm{y}_t$ and $h_{i,t+1}$, {\it i.e.}, $\rho_i^*=Corr(z_{1t}, h_{i,t+1} )=\lambda_{ii}/\sqrt{\lambda_{ii}^2+\psi_{ii}}$ for $i=1,\ldots,9$. The posterior means of $\rho_i^*$ are estimated to be negative ranging from $-0.22$ to $-0.15$. If we regard $z_{1t}$ as the market factor, a decrease in the market return ($z_{1t}$) is followed by an increase in the log volatility ($h_{i,t+1}$), which implies the existence of the leverage effect. The estimation results for $\mathbf{\Psi}$ are omitted in order to save space where the posterior probability that $\psi_{ij}>0$ is found to be greater than 0.975 for all $i$s and $j$s.
To investigate whether the number of factors is $q=1$, we also fit the proposed model using $q=2$ and we set $\mathbf{\Lambda}=[\bm{\lambda}_1,\bm{\lambda}_2,\bm{0},\ldots,\bm{0}]$ for the leverage effect. Table (ref) shows the posterior means, the 95% credible intervals and inefficiency factors for $\rho_{1i}^*=Corr(z_{1t}, h_{i,t+1})$ and $\rho_{2i}^*=Corr(z_{2t}, h_{i,t+1})$ for $i=1,\ldots,9$ where $\bm{z}_t=\mathbf{R}_t^{-1/2}\mathbf{V}_t^{-1/2}(\bm{y}_t-\bm{m}_t)$. The estimation results for $\rho_{1i}^*$ are almost the same as those for $\rho_i^*$ in Table (ref). Conversely, the posterior means of $\rho_{2i}^*$ are close to zeros, and the 95% credible intervals include zero. This suggests that one factor ($q=1$) is enough to describe the leverage effect for our dataset.
{\it Cholesky and spectral decompositions for computing $\mathbf{R}_t^{-1/2}$}. We also estimated our proposed models with $q=1$ and 2 using the Cholesky decomposition instead of the spectral decomposition. The estimation results using the Cholesky decomposition are very similar to those using the spectral decomposition (and hence are omitted) except for the parameters of the leverage effect. Table (ref) shows the estimation results for the correlation, $\rho_i^*$, with $q=1$. All posterior means are estimated to be negative and the posterior probability that $\rho_i^*$ is negative is greater than 0.975 for all $i$s. However, we note that the absolute values of $\rho_i^*$ are smaller than those in the model using the spectral decomposition.
Table (ref) shows the estimation results for the correlations, $\rho_{1i}^*$ and $\rho_{2i}^*$, with $q=2$. The estimation results for $\rho_{1i}^*$ are similar to those for $\rho_i^*$ in Table (ref), but all posterior means of $\rho_{2i}^*$ are estimated to be negative, and the posterior probability that $\rho_{2i}^*$ is negative is greater than 0.975 for $i=2,3$ and $7$. This implies that we need to include more factors when we use the Cholesky decomposition. Since the number of factors $q$ depends on the order of the asset return, we have to find the order that minimizes $q$ for the parsimonious specification. Conversely, the spectral decomposition does not depend on the order of the asset returns and it is much faster at finding a parsimonious specification. We will compare these models using different decompositions with regard to their portfolio performances.
In order to compare the forecasting performance of our proposed models and other existing models, we consider the minimum-variance portfolio strategy (see Han2006). We denote the conditional mean and the conditional covariance matrix of the stock return $\bm{y}_{t+1}$ given the information set $\operatorname{\mathcal{F}}_t$ at time $t$ as
Let $r_{p,t+1}$ denote the portfolio return at time $t+1$. Further, we denote the conditional mean and conditional variance of $r_{p,t+1}$ given the information set $\operatorname{\mathcal{F}}_t$ at time $t$ by
where $r_f$ is the risk free asset return, and $\operatorname{\mbox{\boldmath $w$}}_t$ is a portfolio weight vector for the stock return $\bm{y}_{t+1}$. In the minimum-variance strategy, we minimize the conditional variance $\sigma^2_{p, t+1}$ for the target level $\mu_{p}^*$ of the conditional expected return $\mu_{p,t+1}$. Then the optimal weight $\operatorname{\mbox{\boldmath $w$}}_t$ is given by
The portfolio performances are compared based on the rolling forecast:
To compute the optimal weight $\hat{\bm{w}}_t$, we also need the estimates of $\operatorname{\mbox{\boldmath $m$}}_{t+1 \vert t}$ and $\operatorname{\mathbf{\Sigma}}_{t+1 \vert t}$. We let $N$ denote the number of MCMC iterations, and $(\theta^{(i)}, \{\operatorname{\mbox{\boldmath $h$}}_t^{(i)} \}_{t=1}^{T}, \{\operatorname{\mathbf{R}}_t^{(i)} \}_{t=1}^{T}, \{\operatorname{\mbox{\boldmath $m$}}_t^{(i)} \}_{t=1}^{T} )$ denote the $i$-th MCMC sample ($i=1,\ldots,N$). Using $\operatorname{\mbox{\boldmath $m$}}_{t+1 \vert t}^{(i)}, \operatorname{\mathbf{V}}_{t+1\vert t}^{(i)}, \operatorname{\mathbf{R}}_{t+1 \vert t}^{(i)}, \operatorname{\mathbf{\Sigma}}_m^{(i)}$, we estimate $\operatorname{\mbox{\boldmath $m$}}_{t+1 \vert t}$ and $\operatorname{\mathbf{\Sigma}}_{t+1 \vert t}$ by
In our empirical study, we set $N=1500$ and we discard $500$ samples as the burn-in period for each MCMC rolling estimation (Steps 2 and 3)\footnote{The number of samples being discarded as the burn-in period is sufficient after we obtain the MCMC posterior samples from the previous sample period since we use the posterior means of the parameters and latent variables for the initial values of the next MCMC runs.}. We compare the following multivariate stochastic volatility models as follows.
{\it Cumulative realized objective functions.} Table (ref) shows the cumulative values of the realized objective functions. The MRSV-L1-S model outperforms the other models. Among the MRSV models, the models with leverage outperform the models without leverage, thereby indicating the existence and the importance of the leverage effect. If we assume a constant mean for $\bm{y}_t$, the performance becomes poor in this prediction period, which implies that the random walk process is more flexible for describing the dynamics of the mean of the return vector. The MRSV-L2-C model outperforms the MRSV-L1-C model, but its performance is not as good as that of the MRSV-L1-S model. We could improve the performance of the MRSV models using the Cholesky decomposition by changing the order of the assets or increasing the number of nonzero columns of $\mathbf{\Lambda}$, but it is more efficient to use the spectral decomposition to compute $\mathbf{R}_t^{-1/2}$. Finally, in comparison with the DCC-GARCH model, the HEAVY model, the HAR model and the equally weighted model, we found that the classes of the MSV and MRSV models perform much better.
{\it Time series plots of the portfolio weights}. Figure (ref) shows the time series plots of the portfolio weights in the MRSV-L1-S. The weights for Exxon Mobil are large among all stocks, but towards the end of the period, the weights for IBM and Microsoft tend to become large. However, the weights for the FF rate ($1-\sum_{i=1}^9w_{it}$) are the largest throughout the forecasting period.\\
{\it Comparison of performances before and after the financial crisis}. In order to illustrate the portfolio performances of the models (excluding the CRSV model) in more detail before and after the financial crisis, we divide the forecasting periods into two subperiods: (1) Jan 9, 2008 - July 31, 2008 and (2) Aug 1, 2008 - Dec 31, 2009. As shown in Table (ref), in both subperiods (1) and (2), the portfolio performances are similar to those that we found in the whole period.
{\it Remark} 4. As suggested by an anonymous referee and the Editor, we conducted a predictive ability test based on GiacominiWhite(06) in order to investigate whether the realized objective function of each model is significantly different from that of the MRSV-L1-S in Tables 9 and 10. We found that the differences are all significant except for the MRSV-L1-C models ($\mu_p^*=0.004,0.01,0.1$) and MRSV-L2-C models ($\mu_p^*=0.004,0.1$) in the subperiod (1) (which is before the financial crisis).
The multivariate SV model with flexible dynamic correlation structures that uses the Markov chain Monte Carlo estimation method is proposed. By making full use of the realized variances and realized pairwise correlations, we obtain stable parameter estimates where the covariance matrices are guaranteed to be positive definite. The spectral decomposition is used for the correlation matrices in order to avoid the arbitrariness of the ordering of asset returns. The parsimonious specification for the leverage effect is also proposed. Our models are applied to the daily returns of nine U.S. stocks with their realized volatilities and pairwise realized correlations and are shown to outperform the existing models with regard to portfolio optimizations under a minimum-variance strategy.\\ {\bf Acknowledgements}\\ We thank anonymous referees, the Editor, John Maheu, Hideo Kozumi, Masahiko Sagae and Shuji Tanaka for providing useful comments and discussions. The computational results were obtained by using Ox version 7 (Doornik(06)). This work was supported by JSPS KAKENHI Grant Numbers 25245035, 26245028.