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.
106,409 characters · 23 sections · 94 citation commands
A Multivariate Realized GARCH Model
{\smallKeywords:}{ Multivariate GARCH, Realized Measures, Correlation Factor Model, Block Correlation Matrix}
{\smallJEL Classification:}{ G11, G17, C32, C58 }
Univariate GARCH models have enjoyed considerable empirical success since the ARCH model was introduced by Engle:1982. Subsequently, many univariate GARCH-type models have been proposed in the literature, whereas the research on multivariate GARCH models is less extensive. Generalizing a univariate GARCH model to higher dimensions involves several choices, and it is not always clear how to do this in the most natural way. One challenge to modeling the conditional covariance matrix, $H_{t}=\mathrm{var}(r_{t}|\mathcal{F}_{t-1})$, is the need for $H_{t}$ to be positive (semi) definite. This requirement amounts to nonlinear restrictions across all the elements of $H_{t}$. Another obstacle is that the number of covariance terms increases with $n^{2}$, where $n$ is the dimension of the system. This becomes computationally challenging unless $n$ is relatively small. Multivariate GARCH models have been introduced to address these issues, see BauwensLaurentRombouts:2006, SilvennoinenTerasvirta:2009, and FrancqZakoian2019 for reviews of this literature. In this paper, we adopt a novel approach that guarantees a positive definite $H_{t}$, while the complexity of the model can be contained with a simple factor model. We model the conditional variances and the conditional correlations separately, similar to the Dynamic Conditional Correlation (DCC) model by Engle2002. See also EngleSheppard:2001, PakelShephardSheppardEngle:2021, Aielli:2013, and EngleLedoitWolf:2019, and see EngleKelly:2012 for the DCC variants known as DECO (Dynamic Equicorrelation Correlation) and Block-DECO.
Our main contributions are the following. First, we develop a new class of multivariate GARCH models that facilitate flexible modeling of the correlation structure, while positive definiteness is assured as an innate property. We refer to these as Multivariate Realized GARCH (MRG) models. The main methodological contribution is the dynamic model for the correlation matrix, which can accommodate a simple factor structure and utilize realized measures of correlations in the modeling. Conveniently, the factor approach can greatly reduce the number of latent variables and parameters to be estimated. Second, we show that this factor structure arises naturally with block correlation matrices, including equicorrelation matrices. A block correlation structure may be motivated by the sector classification for companies, or some other partitioning of assets into clusters. Not only is the block correlation specification equivalent to a simple linear factor model, the log-likelihood function is also readily available and easy to evaluate, even for high dimensions. Third, we demonstrate the usefulness of the framework in an empirical application with nine assets. We find that the MRG model improves empirical fit, both in-sample and out-of-sample, relative to the DCC, DECO, and Block-DECO models. The predicted covariance matrices can be used for portfolio construction, such as variance minimization. In an out-of-sample comparison, we find that the empirical portfolio variance is reduced by a factor of two relative to the equal-weighted portfolio. Fourth, we make an interesting auxiliary empirical observation. We find the vector representation of the realized correlation matrix is approximately Gaussian distributed. This result is analogous to existing results for the logarithmically transformed realized variances, see ABDE:2001 and ABDL:2001. An auxiliary result of our analysis is a framework that makes it possible to estimate the Block-DECO model by maximum likelihood for any number of blocks.
Early multivariate GARCH models relied solely on daily returns to update the conditional covariance matrix. The MRG model, introduced in this paper, incorporates realized measures of volatilities and correlations computed from high frequency data. Realized measures are beneficial because they provide accurate signals for dynamic modeling of conditional variances and correlations. These measures gained prominence following their empirical applications in AndersenBollerslev:1998a and subsequent theoretical results by ABDL:2001, BNS:2002, ABDL:2003, BNS:2004, see also HansenLundeVolForecastingHandbook and references therein. Realized measures were initially used to evaluate the performance of GARCH models, see AndersenBollerslev:1998a. A very natural progression was to incorporate realized measures into GARCH models. This was explored in Engle2002b who found that adding the realized variance as an exogenous variable, leads to significant improvements in the empirical fit. This development was followed by more comprehensive models that specified dynamic processes for the realized measures themselves, including the MEM by engle-gallo:06, the HEAVY model by ShephardSheppard:2010, and the Realized GARCH model by HansenHuangShek:2012. Multivariate extensions of these models were proposed in NoureldinShephardSheppard:2012, HansenLundeVoev:2014, and GorgiHansenJanusKoopman:2019. Another way to incorporate realized measures in multivariate GARCH models is explored in BauwensStortiViolante:2012, who build on the Conditional Autoregressive Wishart model of Golosnoy_Gribisch_Liesenfeld_2012.
Our approach to modeling correlations could, with some modifications, be implemented using only daily returns. However, incorporating realized measures into the modeling offers significant advantages. For example, including realized measures makes a GARCH model more responsive to sudden shifts in volatility, leading to substantial improvements in empirical fit and predictive performance, see HansenHuang:2016. The Realized GARCH framework facilitates the incorporation of realized measures of volatility into modeling. The model proposed in this paper represents the first multivariate generalization of the Realized GARCH framework, facilitating the incorporation of realized measures of correlations without imposing additional restrictions on the covariance structure..
The new class of multivariate GARCH models is based on the vector parametrization, $\gamma_{t}=\gamma(C_{t})$, where the mapping $\gamma=\gamma(C)$ is defined by stacking the $d=n(n-1)/2$ below-diagonal elements of $\log C$ (the matrix logarithm of $C$) into the vector $\gamma$, see ArchakovHansen:Correlation. The mapping $C\mapsto\gamma(C)$ is one-to-one between the set of non-singular correlation matrices and $\mathbb{R}^{d}$. So, any vector $\gamma\in\mathbb{R}^{d}$ will map to a unique positive definite correlation matrix, $C(\gamma)$, without the need for additional restrictions. However, it is easy to impose additional structure on $\gamma$ (i.e. structure on $C(\gamma)$) as we will demonstrate with a factor model. This structure makes it possible to estimate the model with a large number of assets.
We are not first to use the matrix logarithm in this context. For instance, ChiuLeonardTsui:1996, Kawakatsu:2006, and AsaiSo:2015 applied the matrix logarithm to covariance matrices. The transformation has also been used in stochastic volatility models, see IshiharaOmoriAsai:2016, and in reduced-form models of realized covariance matrices, see e.g. BauerVorkink:2011 and Weigand:2014.\footnote{Additional related literature includes the work by Liu_2009, ChiriacVoev:2011, Golosnoy_Gribisch_Liesenfeld_2012, and BauwensStortiViolante:2012.} The logarithmic transformation of conditional covariance matrices is also related to the dynamic eigenvalue model by HetlandPedersenRahbek:2023. Our approach, which applies the matrix logarithm to the correlation matrix, allows us to model each of the conditional variances with univariate GARCH models, which has additional benefits. For instance, it enables us to explicitly model the empirically important leverage effect. HafnerWang:2023 have recently proposed a dynamic model of correlations that also uses the parametrization by ArchakovHansen:Correlation. Their model uses the score-driven framework by CrealKoopmanLucas:2013, whereas we build on the Realized GARCH framework and develop a parsimonious factor model for the correlation structure.
We proceed as follows. In Section (ref), we introduce notation, the modeling framework, and discuss how the factor structure can be imposed on the correlation matrix. Section (ref) details the estimation of the model and how the model can be used for forecasting. An extensive empirical analysis with nine asset returns series from three economic sectors is presented in Section (ref). We compare the new model with existing models using out-of-sample criteria in Section (ref), we conclude in Section (ref). We derive analytical expressions for the derivatives of the log-likelihood estimation in Appendix A. These greatly speed up the estimation of the model. Appendix B has step-by-step directions for maximum likelihood estimation of the model.
In this section, we present the details of the model. The MRG is based on the vector parametrization of the correlation matrix, \[ \gamma=\gamma(C)=\mathrm{vecl}(\log C), \] where $\log C$ represents the logarithmically transformed correlation matrix\footnote{For a nonsingular correlation matrix, we have $\log C=Q\log\Lambda Q^{\prime}$, where $C=Q\Lambda Q^{\prime}$ is the spectral decomposition of $C$, so that $\Lambda$ is a diagonal matrix with the eigenvalues of $C$.} and $\mathrm{vecl}(\cdot)$ extracts and vectorizes the elements below the diagonal. To illustrate this parametrization, consider the following example,
In the bivariate case, $n=2$, we have $\gamma(C)=\frac{1}{2}\log\frac{1-\rho}{1+\rho}$, which is the Fisher transformed correlation. This parametrization was used in HansenLundeVoev:2014 and the model we propose here can therefore be viewed as a natural generalization of the bivariate structure in HansenLundeVoev:2014.\footnote{HansenLundeVoev:2014 accommodated the case $n>2$ by fusing bivariate models to a larger system, which induces a restricted structure on $C_{t}$. }
Our theoretical results also have useful applications for existing models. For instance, they make it straightforward to estimate the Block-DECO model by EngleKelly:2012 by maximum likelihood. In EngleKelly:2012 the correlations within each block were obtained by averaging over estimated correlations, which were based on an auxiliary DCC model for the full dimension. They derived an expression for the log-likelihood for the case with $K=2$, but not for $K>2$. Instead they proposed to use composite likelihood methods when $K>2$. The canonical representation of block matrices by ArchakovHansen:CanonicalBlockMatrix makes it straightforward to evaluate the log-likelihood function for any $K$. The problem is effectively simplified to involve a single, low-dimensional $K\times K$ matrix. Another advantage of the parametrization, $\gamma(C)$, is that it is very simple to specify a low-dimensional factor model that is equivalent to a block structure in $C$. This avoids having to take averages of elements from an auxiliary and high-dimensional DCC model.
We let $r_{t}=(r_{1,t},\ldots,r_{n,t})^{\prime}$ denote a $n$-dimensional vector of returns in period $t$, where $t$ represents a generic unit of time, such as a trading day. The conditional mean is denoted by $\mu_{t}=\mathbb{E}(r_{t}|\mathcal{F}_{t-1})$ and the conditional variance by $H_{t}=\mathrm{var}(r_{t}|\mathcal{F}_{t-1})$, where $\{\mathcal{F}_{t}\}$ is the natural filtration for $(r_{t},\mathrm{RM}_{t})$. Here $\mathrm{RM}_{t}$ denotes an ex-post empirical measure of $H_{t}$, such as the realized covariance matrix, see BNS:2004, or the multivariate realized kernel by BNHLS-MRK:2011.
We decompose the conditional covariance matrix into variances and correlations,
where $\Lambda_{h_{t}}=\mathrm{diag}(h_{1,t},\ldots,h_{n,t})$ with $h_{i,t}=[H_{t}]_{ii}$, $i=1,\ldots,n$. So, $h_{i,t}$ is the conditional variance of $r_{i,t}$, $i=1,\ldots,n$, and $C_{t}=\mathrm{corr}(r_{t}|\mathcal{F}_{t-1})$ is the conditional correlation matrix of $r_{t}$. The DCC structure in ((ref)) enables us to disentangle the dynamic properties of the conditional variances from that of the conditional correlations. We will deviate from the DCC framework in the way we parametrize $C_{t}$ and by incorporating realized measures of variances and correlations into the model.
The central component of a GARCH model is the equation that specifies the dynamic properties of $H_{t}$ and how these are influenced by lagged returns. This equation can be enhanced to include realized measures of volatility. The Realized GARCH model is characterized by measurement equations that specify how the realized measures are related to the contemporaneous conditional moments, i.e. the elements of $H_{t}$.
Here we will model the conditional variances and the conditional correlations separately. This leads to two sets of GARCH and measurement equations, that utilize the appropriate realized measures. From $\mathrm{RM}_{t}\in\mathbb{R}^{n\times n}$, a positive definite realized measure of the covariance matrix in period $t$, we extract the diagonal elements $x_{t}=\mathrm{diag}(\mathrm{RM}_{t})$ and the corresponding correlation matrix denoted by \[ Y_{t}=\Lambda_{x_{t}}^{-1/2}\mathrm{RM}_{t}\Lambda_{x_{t}}^{-1/2}. \] Here $\Lambda_{x_{t}}=\mathrm{diag}(x_{1,t},\ldots,x_{n,t})$ denotes the diagonal matrix with the elements of $x_{t}$ on the diagonal, and it follows that $Y_{t}$ inherits the positive definiteness from $\mathrm{RM}_{t}$, such that $y_{t}=\gamma(Y_{t})\in\mathbb{R}^{d}$ is well-defined. In summary, $x_{t}$ and $y_{t}$ are the observed empirical measures of the latent variables, $h_{t}$ and $\gamma_{t}$, respectively, such that time variation in $x_{t}$ and $y_{t}$ contains information about time variation in the corresponding latent variables. The realized measure, $\mathrm{RM}_{t}$, is typically a consistent estimator of the quadratic variation. But the (ex-post) quadratic variation is not identical to the (ex-ante) conditional variance, $H_{t}$. We should therefore expect a non-trivial measurement errors in $x_{t}$ and $y_{t}$, because neither are perfect measurements of $h_{t}$ and $\gamma_{t}=\gamma(C_{t})$, respectively.
We let $I_{n}$ denote the $n\times n$ identity matrix and $1_{\{\cdot\}}$ the indicator function, which equals one if the expression within the curly brackets is true and zero otherwise. With the required notation in place, we are now ready to introduce the multivariate realized GARCH (MRG) model.
We specify univariate Realized GARCH models for each return series and the corresponding realized measures.
The return equation for each of the returns series takes the form
We have here assumed that $\mathbb{E}(r_{i,t}|\mathcal{F}_{t-1})=\mu_{i}$ is constant, as is often done in GARCH models, and it follows that the standardized return, $z_{i,t}=h_{i,t}^{-1/2}(r_{i,t}-\mu_{i})$, is such that $\mathbb{E}(z_{i,t}|\mathcal{F}_{t-1})=0$ and $\mathrm{var}(z_{i,t}|\mathcal{F}_{t-1})=1$. Note that the standardized returns, $z_{i,t}$ and $z_{j,t}$ for $i\neq j$, are not assumed to be uncorrelated (nor are they likely to be).
The corresponding GARCH and measurement equations are given by
for $i=1,\ldots,n$ and $t=1,\ldots,T$, where $\omega_{i},\beta_{i},\alpha,\xi_{i},\varphi_{i}\in\mathbb{R}$. The measurement errors, $v_{t}=(v_{1,t},\ldots,v_{n,t})^{\prime}$, may be dependent and correlated with the corresponding measurement errors in transformed realized correlations (defined below). The two functions, $\tau_{i}(z)=\tau_{i,1}z+\tau_{i,2}(z^{2}-1)$ and $\delta_{i}(z)=\delta_{i,1}z+\delta_{i,2}(z^{2}-1)$, are leverage functions that capture dependencies between returns and volatility innovations. This dependency is known to be empirically important, and it is typically labelled the leverage effect, see Black:1976, Christie:1982, Engle_Ng_1993.\footnote{The leverage effect is sometimes used to refer to a linear dependence, i.e. the (usually negative) correlation between returns and changes in return volatility.} This quadratic form is motivated by results in HansenHuangShek:2012 and HansenHuang:2016 who found that a second-order Hermite polynomial suffices for capturing the asymmetry dependence between return shocks and volatility shocks. The three equations, ((ref))-((ref)), define a univariate Realized GARCH model for each asset, $i=1,\ldots,n$, where the characteristic feature of a Realized GARCH model is the measurement equation that relates each realized measure to the corresponding conditional moment in the model.
The key novel innovation of the Multivariate Realized GARCH model is the way we model the dynamic conditional correlation matrix, \[ C_{t}=\mathrm{var}(z_{t}|\mathcal{F}_{t-1}). \] We simply model the elements of $\gamma_{t}=\gamma(C_{t})\in\mathbb{R}^{d}$, or linear combinations thereof, in the same way we model $h_{i,t}$, $i=1,\ldots,n$. In its most general specification (Full) we specify a GARCH equation and a measurement equation for each element, $j=1,\ldots,d$,
Here $y_{t}=\gamma(Y_{t})$ is the transformed realized correlation matrix, such that $y_{j,t}$ is the appropriate empirical measurement of $\gamma_{j,t}$, $j=1,\ldots,d$. The measurement error in the transformed realized correlations, $\tilde{v}_{t}=(\tilde{v}_{1,t},\ldots,\tilde{v}_{d,t})^{\prime}$, may have dependent elements and may be correlated with the measurement errors in, $v_{t}$, in ((ref)).
A drawback of modeling all elements of $\gamma_{t}$ is that the number of latent variables in $C_{t}$, $d=n(n-1)/2$, becomes unmanageable unless $n$ is small. While it is possible to estimate the model with $n=9$ ($d=36$), but (at the time of writing this) it is difficult to estimate the full model with dimensions much larger than that. This necessitates additional structure on the model when $n$ is large. Fortunately, it is simple to impose a factor structure, $\gamma_{t}=\varrho(\zeta_{t})$, where $\zeta_{t}\in\mathbb{R}^{r}$ is a lower dimensional vector of factors. The underlying assumption is that the variation in $C_{t}$ is driven by $r<d$ factors.
A natural starting point is the linear factor model, \[ \gamma_{t}=A\zeta_{t}, \] where $A$ is a $d\times r$ matrix. This enables us to reduce the number of GARCH equations and measurement equations from $d$ to $r$, by replacing ((ref)) and ((ref)) with
respectively, where $\check{y}_{t}=(A^{\prime}A)^{-1}A^{\prime}y_{t}\in\mathbb{R}^{r}$. From well-known projection arguments, it follows that $\check{y}_{t}$ is the realized quantity that corresponds to $\zeta_{t}$.\footnote{If $A$ has full column rank, $r$, then there exists an $d\times(d-r)$ matrix, $A_{\bot}$ so that $(A,A_{\bot})$ is a full rank matrix, and $A^{\prime}A_{\bot}=0$. Thus, $\gamma_{t}=A\zeta_{t}$ implies $A_{\bot}^{\prime}\gamma_{t}=0$ and the identity $I_{d}=A(A^{\prime}A)^{-1}A^{\prime}+A_{\bot}(A_{\bot}^{\prime}A_{\bot})^{-1}A_{\bot}^{\prime}$ shows that $\gamma_{t}=A\zeta_{t}\Rightarrow\gamma_{t}=A(A^{\prime}A)^{-1}A^{\prime}\gamma_{t}$. The vector of transformed realized correlations, $y_{t}$, is our empirical “signal” about $\gamma_{t}$, and the identity, \[ y_{t}=[A(A^{\prime}A)^{-1}A^{\prime}+A_{\bot}(A_{\bot}^{\prime}A_{\bot})^{-1}A_{\bot}^{\prime}]y_{t}=A\check{y}_{t}+A_{\bot}(A_{\bot}^{\prime}A_{\bot})^{-1}A_{\bot}^{\prime}y_{t}, \] shows that $\check{y}_{t}$ is the appropriate signal about $\zeta_{t}$, whenever $\gamma_{t}=A\zeta_{t}$.} The measurement errors, $\check{v}_{t}=(\check{v}_{1,t},\ldots,\check{v}_{r,t})^{\prime}$, in ((ref)) may have dependent elements and be correlated with the measurement errors, $v_{t}$, in ((ref)). In the empirical analysis, we define $u_{t}=(v_{t}^{\prime},\tilde{v}_{t}^{\prime})^{\prime}$ and adopt a Gaussian likelihood function with $u_{t}\sim iidN(0,\Sigma)$, where $\Sigma$ is an arbitrary covariance matrix.
The matrix, $A$, is needed for this implementation and $A$ may be known in advance or can be determined empirically. Estimating $A$, including $r$, is an interesting problem that we leave for future research.\footnote{We also leave other generalizations, such as non-linear factor structures, $\gamma_{t}=\varrho(\zeta_{t})$, for future research.} In this paper, we focus on the case where $A$ is known, and we show that any block correlation structure is equivalent to a linear factor structure with a known $A$-matrix.
Note that the linear factor model, $\gamma_{t}=A\zeta_{t}$, is characterized by the subspace spanned by the columns of $A$, not the particular choice for $A\in\mathbb{R}^{n\times r}$. This follows from the fact that $\gamma_{t}=\tilde{A}\tilde{\zeta}_{t}$ and $\gamma_{t}=A\zeta_{t}$ are observationally equivalent whenever $\tilde{A}=A\Phi$ and $\Phi\in\mathbb{R}^{r\times r}$ is invertible. (The relation between the two latent factor variables will be $\tilde{\zeta}_{t}=\Phi^{-1}\zeta_{t}$).
There are situations where the original measurement equation in ((ref)) should be used even when the factor structure is adopted in the GARCH equation ((ref)). For instance, if multiple model specifications are compared in terms of their total log-likelihoods, then all models should have measurement equations for the same realized variables to avoid an apples-to-oranges comparison. Models with different measurement equations may be compared in terms of their partial log-likelihood for returns, which is often the primary objective in multivariate GARCH models. The relevant terms of the log-likelihood for these comparisons are detailed in the next section.
Let $n=n_{1}+\cdots+n_{K}$, where $n_{1},n_{2},\ldots,n_{K}\geq1$. A block correlation matrix is characterized by the structure, \[ C=\left[
\right], \] where all elements within each $n_{i}\times n_{j}$ matrix, $C_{[i,j]}$, are identical and equal to $\rho_{ij}$, except for the the diagonal-block, $C_{[i,i]}$, that have ones along the diagonal and other elements equal to $\rho_{ii}$.
Regardless of the number of blocks and their dimensions, $\log C_{t}$ has the same block structure as $C_{t}$, see ArchakovHansen:CanonicalBlockMatrix.\footnote{This result holds for all block matrices, including non-symmetric block matrices, see ArchakovHansen:CanonicalBlockMatrix.} This is illustrated in the following example:
A block structure arises naturally in applications where the correlation between two variables is defined by their group classification. An $n\times n$ correlation matrix has $n(n-1)/2$ correlations, whereas a block correlation matrix with $K\times K$ blocks has, at most, $K(K-1)/2+K$ distinct correlations.\footnote{The exact number of distinct correlation is $r=K(K-1)/2+\tilde{K}$, where $\tilde{K}$ is the number of blocks that contain two or more elements. The reason for the distinction between $K$ and $\tilde{K}$ is that an $1\times1$ diagonal block does not have a correlation coefficient.} The fact that $\log C$ inherits the block structure of $C$ facilitates a parsimonious modeling of dynamic block correlation matrices. The elements of the transformed matrix can be modeled in an unrestricted way, without compromising the structure a correlation matrix must have. This completely bypasses the non-linear cross restrictions on $C$'s elements ensure positive definiteness. So, a dynamic model of the (below diagonal) elements of $\log C$ is a very convenient implementation of the block structure on $C$.
The block structure offers a useful dimension reduction, but in order to make use of it, one has to specify the clusters that define the block structure. In our empirical analysis, we consider a case with 9 assets, and a correlation structure with $K=1$, $K=3$, and $K=9$. The sector-based clusters, $K=3$, reduces the number of latent correlation variables from $d=36$ to $r=6$.
A correlation matrix with a block structure leads to a very scalable model, because increasing $n$ does not increase the complexity of the second stage estimation. The correlation factors, $\zeta_{t}\in\mathbb{R}^{r}$, can be used to model an arbitrarily large number of assets, so long as all assets can be classified within the $K$ clusters.
The block correlation structure also simplifies many computational aspects of the model. For instance, the inverse correlation matrix is readily available for any $K$, see ArchakovHansen:CanonicalBlockMatrix. Previously, a closed-form expressions for $C^{-1}$ was only available for $K=2$, see EngleKelly:2012. For a block correlation matrix, $C$, with $K\times K$ blocks, we define the (symmetric) $K\times K$ matrix, $B$, whose elements are given by $b_{ii}=1+(n_{i}-1)\rho_{ii}$ and $b_{ij}=\rho_{ij}\sqrt{n_{i}n_{j}}$, for $i\neq j$. It is simple to verify that $(i,j)$-th block of $C$ can be expressed as \[ C_{[i,j]}=b_{ij}P_{[i,j]}+1_{\{i=j\}}(1-\rho_{ii})(I_{n_{i}}-P_{[i,i]}),\qquad\text{for}\quad i,j=1,\ldots,K, \] where all elements of $P_{[i,j]}\in\mathbb{R}^{n_{i}\times n_{j}}$ equal $\frac{1}{\sqrt{n_{i}n_{j}}}$. The determinant of $C$ and the $(i,j)$-th block of the inverse correlation matrix, $C^{-1}$, can be expressed as
respectively, where $b_{ij}^{\#}$ is the $ij$-th element of $B^{-1}$, see ArchakovHansen:CanonicalBlockMatrix. These closed-form expressions for $C^{-1}$ and $\det C$ facilitate simple evaluation of the Gaussian log-likelihood function when $C$ has a block structure.
Estimating the parameters of the MRG model is relatively simple, because the model is observation-driven, with variation in the latent variables, $h_{t}$ and $\gamma_{t}$, being driven by observable variables, $r_{t}$, $x_{t}$, and $y_{t}$.
We can factorize the joint density of $(r_{t},x_{t},y_{t})$, conditional on past observations, into the marginal density for returns and the density for realized variables, conditional on contemporaneous returns. Thus, the joint density is expressed as the product, $f_{t-1}(r_{t},x_{t},y_{t})=f_{t-1}(r_{t})f_{t-1}(x_{t},y_{t}|r_{t})$. The log-likelihood function can therefore be deduced from
and parameters may be estimated by quasi maximum likelihood estimation, by specifying Gaussian likelihood functions for $z_{t}$ and $u_{t}$. More specifically, under the assumption that $C_{t}^{-1/2}z_{t}\sim iidN(0,I_{n})$ and $u_{t}\sim iidN(0,\Sigma)$ are mutually independent. The leverage functions, $\tau_{i}(\cdot)$ and $\delta_{i}(\cdot)$, serve to eliminate certain forms of dependence between $v_{i,t}$ and $z_{i,t}$, making the assumed independence somewhat more realistic.
Let $\theta=(\theta_{1}^{\prime},\theta_{2}^{\prime})^{\prime}$ represent all unknown parameters in the model, where $\theta_{1}$ includes the parameters in the multivariate GARCH-X model for the returns and $\theta_{2}$ represent the parameters in the measurement equations, which define the model for the realized measures. The Gaussian specification and ((ref)) imply that the log-likelihood function is given by \[ \ell(\theta)=\sum_{t=1}^{T}\ell_{r,t}(\theta_{1})+\sum_{t=1}^{T}\ell_{x,y|r,t}(\theta_{2}), \] where
with $c_{n}=n\log2\pi$ and
Here we have used a condensed notation, where $\xi$ is a vector with elements, $\xi_{i}$, $i=1,\ldots,n$, and similar for $\delta$ and $\tilde{\xi}$, and $\Phi$ and $\tilde{\Phi}$ are diagonal matrices with diagonal elements $(\varphi_{1},\ldots,\varphi_{n})$ and $(\tilde{\varphi}_{1},\ldots,\tilde{\varphi}_{d})$, respectively.
For a particular value of $\theta$ it is straightforward to evaluate the log-likelihood function. The latent variables, $\{h_{t},C_{t}\}$, can be computed recursively. Given $h_{t-1}$ and $C_{t-1}$ and the observable $(r_{t},x_{t},y_{t})$, we compute $h_{t}$ and $C_{t}$ with the GARCH equations and infer $z_{t}$ from the return equation, ((ref)), and $u_{t}$ from the measurement equations. This is repeated for period $t+1$, and so forth, and the likelihood function can be evaluated for the full sample. The starting values for the latent variables, $h_{1}$ and $C_{1}$, can be treated as unknown parameters (as part of $\theta$), which we recommend. Alternatively, $h_{1}$ and $C_{1}$, can be assigned particular values, which may be based on appropriate empirical quantities.
The structure of the log-likelihood function allows the maximization problem to be simplified. Given the residuals, $\hat{u}_{t}$, $t=1,\ldots,T$, it can be shown that the maximum likelihood estimator of $\Sigma$ is $\hat{\Sigma}=T^{-1}\sum_{t=1}^{T}\hat{u}_{t}\hat{u}_{t}^{\prime}$, such that \[ \sum_{t=1}^{T}\hat{u}_{t}^{\prime}\hat{\Sigma}^{-1}\hat{u}_{t}=\mathrm{tr}\{\hat{\Sigma}^{-1}\sum_{t=1}^{T}\hat{u}_{t}\hat{u}_{t}^{\prime}\}=\mathrm{tr}\{TI_{n+k}\}=T(n+k), \] where $I_{n+k}$ is the $(n+k)\times(n+k)$ identity matrix. So, the objective to be maximized is (apart from a constant) given by
where the omitted constant is $-\tfrac{T}{2}(c_{n}+c_{d}+n+k)=-\tfrac{n(n+1)}{4}T\log2\pi-T\tfrac{n+k}{2}$. As stated earlier, both the determinant, $\det C_{t}$, and the inverse, $C_{t}^{-1}$, are simple to evaluate when $C_{t}$ has a block structure using ((ref)) and ((ref)). More details about the maximum likelihood estimation can be found in Appendix (ref).
Joint estimation is possible if $n$ is small (say $n<10$), but it tends to be slow. So, it will often be convenient to adopt a two-stage estimation method. In the first stage, we estimate the univariate Realized GARCH models for each of the $n$ return series. From the estimated models we obtain $h_{i,t}$ and $z_{i,t}=h_{i,t}^{-1/2}(r_{i,t}-\mu_{i})$, $i=1,\ldots,n$ and $t=1,\ldots,T$, and in the second stage, we estimate the parameters that relate to the dynamic conditional correlation matrix, $C_{t}=\mathrm{var}(z_{t})$.
So, in the first stage, we maximize \[ -\frac{1}{2}\sum_{t=1}^{T}\left\{ \log h_{i,t}+\log\hat{\sigma}_{v_{i}}^{2}+z_{i,t}^{2}\right\} , \] with respect to $\vartheta_{1,i}=(\omega_{i},\beta_{i},\tau_{i,1},\tau_{i,2},\alpha_{i},\xi_{i},\varphi_{i},\delta_{i,1},\delta_{i,2})^{\prime},$ for each $i=1,\ldots,n$, where $\hat{\sigma}_{v_{i}}^{2}=\tfrac{1}{T}\sum_{t=1}^{T}v_{i,t}^{2}$, and in the second stage, we maximize
with respect $\vartheta_{2}=(\vartheta_{2,1}^{\prime},\ldots,\vartheta_{2,r}^{\prime})^{\prime}$, where $\vartheta_{2,j}=(\check{\omega_{j}},\check{\beta}_{j},\check{\alpha}_{j},\check{\xi}_{j},\check{\varphi}_{j})^{\prime}$, $j=1,\ldots,r$, are the parameters in ((ref)) and ((ref)), with the vector of transformed realized correlation variables given by $\check{y}_{t}=(A^{\prime}A)^{-1}A^{\prime}y_{t}$ and $u_{t}=(v_{t}^{\prime},\check{v}_{t}^{\prime})^{\prime}$. Note that the two-stage estimation involves a different partitioning of the parameters than the one used in ((ref)). The parameters in the first-stage, $(\vartheta_{1,1},\ldots,\vartheta_{1,n})$, include elements from both $\theta_{1}$ and $\theta_{2}$, and the same is the case for the second-stage parameters, $\vartheta_{2}$.
In the special case, without a factor structure imposed on $C_{t}$, (Full), we simply have $r=d$, and the parameters are identical to those in ((ref)) and ((ref)).
In our empirical analysis, we adopt this two-stage estimation method. All model-specifications use the same first-stage, such that the model comparisons concern their ability to capture the dynamic correlation structure without confounding these with aspects that relate to the marginal distributions.\footnote{In an earlier version of this paper, we estimated the model by maximum likelihood (using a shorter sample period). The parameter estimates were very similar, but the estimation was much slower, especially for the most flexible specification (Full), which has 36 latent variables to model $C_{t}$.}
For the estimation problem, we have derived analytic expressions for the gradient vector and the Fisher information matrix. This are presented in the Appendix. The analytical expressions greatly reduce the time it takes to estimate the model, as illustrated in Table (ref). The Table presents second-stage estimation times for the computationally most demanding specification, with an unrestricted correlation matrix. The model is estimated with $T=4,744$ and $n$ ranging from $2$ to $9$. The variables are subsets of the assets using in our empirical analysis).{
}
Without analytical derivatives (which was used in an earlier version of this paper) the estimation of the full model is relatively slow. Using numerical derivatives, it takes about 78 seconds to estimated the model with $n=3$, more than 30 minutes to estimate the model with $n=6$, and more than four hours to estimate the model with $n=9$. In contrast, it only takes 7 minutes to estimate the model with $n=9$ when the analytical derivatives are used. \footnote{The model was estimated with Matlab's optimization function, fminunc, and the computation of numerical derivatives was optimized to use all (24) cores in parallel. Estimation with analytical expressions uses the quasi-newton method with the analytical gradient in the first 25 iterations, after which the trust-region method is used with the analytical gradient and the analytical approximation for the corresponding Hessian.}
The analytical derivatives, derived in the appendix will also be useful for computing standard errors, testing for parameter stability, see Nyblom89, and other statistical applications.
One-step ahead forecasting of the return distributions from the model is straightforward. All dynamic variables are specified with an observation-driven structure, such that forecasts are given from known functions of lagged variables. From the observed variables in period $t$, all the conditional variances and correlations for period $t+1$ can be computed from the GARCH equations. The elements of $H_{t+h}$, are not predetermined beyond horizon $h=1$, because they also depend on future realizations of $z_{t}$ and $u_{t}$. It is nevertheless straightforward to compute a distributional forecasts for $H_{t+h}$ using simulation or bootstrap methods. So, multi-step ahead forecasts can be inferred from the estimated model for any forecasting horizon. Forecasting schemes for the Realized GARCH models of this kind are detailed in LundeOlesenEnergyRealizedGARCH and HansenLundeVoev:2014. In this context, a bootstrap method will typically be preferred because it is more robust to distributional misspecification.
Our empirical analysis spans a sample period from January 2, 2002 to December 31, 2020, which has 4,744 trading days after the removal of holidays and trading days with reduced trading hours. We use daily close-to-close returns and compute realized variances and correlations from high-frequency data.
We include nine stocks in our analysis. Three stocks from the energy sector, CVX, MRO, and OXY, three stocks from the Health Care sector, JNJ, LLY, and MRK, and three stocks from the Information Technology sector, AAPL, MU, and ORCL.
We construct close-to-close daily returns for the individual stocks using closing prices, adjusted for stock splits and dividends, from the CRSP US Stock Database. Intraday transaction data were obtained from the TAQ database, and these were cleaned in accordance with the methodology detailed in BNHLS-RKpractice:2009. From the high-frequency data, we compute the $9\times9$ multivariate realized kernel estimates for each trading day, $\mathrm{RM}_{t}\in\mathbb{R}^{9\times9}$, and these are used to define $x_{t}\in\mathbb{R}^{9}$ and $y_{t}\in\mathbb{R}^{36}$.
We present summary statistics for the nine return series and their corresponding realized variance measures in Table (ref). These statistics are consistent with typical estimates for such time series. We note that Health Care stocks (middle three columns) had the lowest volatility, whereas IT stocks (the last three columns) had the largest average volatility in the sample period. This can be seen from the standard deviations in the second row, and the means and medians (Q-50%) of the Realized volatilities in the lower part of Table (ref).
Summary statistics for the realized correlations are presented in Table (ref), where the shaded regions illustrate the block structure we use in a sector-based factor model. The numbers below the main diagonal are the average realized correlations (the off-diagonal elements of $\bar{Y}=\tfrac{1}{T}\sum_{t=1}^{T}Y_{t}$) and the numbers above the diagonal are the corresponding averages for the transformed quantities (the off-diagonal elements of $\tfrac{1}{T}\sum_{t=1}^{T}\log Y_{t}$). Note that the realized correlations within each of the blocks have similar averages. The three assets from the energy sector are highly correlated, with correlations of about 0.55 on average. The average within-sector correlations for the Health Care sector and Information Technology sector stocks are about 0.39 and 0.31, respectively. The between-sector correlations tend to be smaller, and these range between 0.17 and 0.28. Similar patterns are seen for the logarithmically transformed correlation variables.
The time series of realized correlations are shown in Figure (ref). The left subplots present the within-sector correlations (gray lines) and their average daily correlation (red line) for each of the three sectors. The right subplots present the between-sector correlations (gray lines) and their daily average (red line) for the three sector pairs. The 36 correlation time series are computed from the $9\times9$ multivariate Realized Kernel estimator. Correlations within each of the six categories tend to move together; however, there are notable differences across the six types of within-sector and between-sector correlations, both in terms of their average level and their variation over time. For instance, the average between-sectors correlation for Health Care returns and Information Technology returns does not have a sharp decline in late 2008, as can be seen for the two other between-sectors correlations involving Energy sector returns. Figure (ref) provides additional motivation for exploring a block structure for the correlation matrix.
The empirical results are based on the two-stage estimation procedure, described in Section (ref), and we use three different specifications for the correlation matrix: equicorrelation (Equi), block correlation (Block), and unrestricted (Full). The simplest model is the equicorrelation model that assumes a single common dynamic correlation coefficient that is common for all correlations in $C_{t}$. The most general specification is the fully dynamic correlation matrix with 36 dynamic correlation coefficients. The equicorrelation model has a single latent variable to model the time variation in the correlation matrix, such that $\zeta_{t}$ is univariate ($r=1)$ in this case. The block correlation model employs six latent variables ($r=6)$, whereas the most flexible specification, Full, has $d=36$ latent variables.
Parameter estimates from the first-stage estimation (univariate Realized GARCH models) are reported in Table (ref). The estimates are based on the full sample period from January 3, 2002, to December 31, 2020. Each column in Table (ref) corresponds to one of nine assets in our analysis. We report parameter estimates along with their corresponding standard errors, shown below in parentheses.\footnote{We obtain standard errors numerically by calculating the Hessian matrix and the information matrix via the outer gradient product.} The model implies an AR(1) model for $h_{i,t}$ with $\pi=\beta+\alpha\varphi$ as the autoregressive coefficient. This can be seen by substituting the measurement equation into the corresponding GARCH equation. In the bottom of Table (ref), we report this persistence parameter for each of the conditional variances. These estimates are all close to unity, as expected, since volatility is known to be persistent. Overall, the parameter estimates are in line with results in the existing literature on GARCH models and Realized GARCH models. In the last row of Table (ref), we report the average partial log-likelihoods, $-2T^{-1}\ell_{r_{i}}=\log2\pi+\frac{1}{T}\sum_{t=1}^{T}\left\{ \log h_{i,t}+z_{i,t}^{2}\right\} $, $i=1,\ldots,9$, that summarize how well the estimated univariate Realized GARCH models describe the conditional (marginal) distribution of returns, for each of the nine asset.
Table (ref) presents the parameter estimates from the second-stage estimation of the Multivariate Realized GARCH model, where the model for the dynamic correlation matrix is estimated. We report parameter estimates for the three specifications. The first column is for the equicorrelation model (Equi), the next six columns are for the six dynamic correlations in the 3$\times$3 block-correlation structure (Block), and the last column is for the unrestricted correlation structure (Full), where we present the range of the estimates from the 36 models for $\gamma_{1,t},\ldots,\gamma_{36,t}$. For instance, the estimated intercepts in the GARCH equation, $\tilde{\omega}_{j}$, ranged between $-0.01$ and $0.026$. We also report the persistence parameter, $\tilde{\pi}=\tilde{\beta}+\tilde{\alpha}\cdot\tilde{\varphi}$, for the conditional correlations, and these are quite similar to those of the conditional variances. Thus, both conditional variances and conditional correlations are found to be persistent.
It is not meaningful to compare the total log-likelihood for models with different measurement equations. It is, however, meaningful to compare their log-likelihood for returns, $\ell_{r}$, see ((ref)), which reflects how well the model describes the conditional distribution of returns. We can decompose $\ell_{r}$ as \[ -2\ell_{r}=\underset{-2\sum_{i=1}^{n}\ell_{r_{i}}}{\underbrace{c_{n}+\sum_{i=1}^{n}\sum_{t=1}^{T}\log h_{i,t}+z_{i,t}^{2}}}+\sum_{t=1}^{T}\log\det C_{t}+z_{t}^{\prime}(C_{t}^{-1}-I)z_{t}, \] where the first term is the log-likelihoods for the marginal distributions of returns (which is common for all model-specifications in our analysis), and the last term captures the effect of different correlation models. We therefore report $-2T^{-1}\ell_{r}$ in Table (ref) along with the corresponding BIC that includes a penalty for model complexity. Each latent variables, $\zeta_{j,t}$, adds five parameters, $(\tilde{\omega}_{j},\tilde{\alpha}_{j},\tilde{\beta}_{j},\tilde{\xi}_{j},\tilde{\varphi}_{j})$, to the complexity of the model.
The equicorrelation structure has an average log-likelihood for returns that is much smaller than those of the more flexible specification. The average daily difference is about one unit, which add up to a substantial difference with $T=4,744$ days in the sample. The difference between the block specification and full specification is more modest at 0.0655 units per day, which adds up to about 310 units of $2\ell_{r}$ over the full sample. This improvement is achieved with 150 additional parameters, and the Bayesian information criterion (BIC) prefers the block specification in this case. In the next section, we evaluate these specifications in out-of-sample comparisons, and compare them to DCC-type benchmark models.
The estimated dynamic equicorrelation time series is presented in Figure (ref), along with the daily average realized correlation (red line). The horizontal blue dashed line is the estimated equicorrelation under the assumption that the correlation is constant over the sample period.
It is comforting that the estimated correlation is in line with with the average realized correlation, and the gradual variation in the series strongly suggests that the conditional correlation is time-varying.
In Figure (ref) we present the corresponding results for the estimated block specification. In the left panels, we present the estimated within-sector correlation (black lines) and the corresponding daily average of the realized correlations (red lines), and in the right panels we present the results for between-sectors correlations. The horizontal dashed line in each of the plots is the estimated block correlation, under the assumption that it is constant over the sample period.
ABDE:2001 and ABDL:2001 found that the logarithmic realized variances for stock returns and exchange rate data are approximately normally distributed. ABDE:2001 also found realized correlations to be approximately normally distributed. Here, we find that the elements of $y_{t}=\gamma(Y_{t})$, which are the transformed realized correlations, are also approximately normally distributed.
Panel (a) in Figure (ref) presents Q-Q plots for the empirical distribution of the transformed realized correlations against the normal distribution. The results are based on the sample period from January 2, 2002 to December 31, 2020, which has 4,744 daily observations. The quantiles of their empirical distributions are plotted against the corresponding quantiles of the normal distribution. The left column of panel (a) has Q-Q plots for the three unique series within each of the three diagonal blocks of $Y_{t}$, and the Q-Q plots in the right column of panel (a) are for the nine series in each of the three off-diagonal blocks of $Y_{t}$. The black dots are based on elements of the transformed realized correlations, $y_{t}$, that are used in the Full specification, and the red dots represent within-block variables, $\check{y}_{t}$, that are the variables used in the Block specification. The plots in panel (b) present the corresponding results for the residuals in the measurement equations, where red dots represent $\tilde{v}_{t}$ and black dots represent $\check{v}_{t}$.
The Q-Q plots show that both transformed realized correlations and the corresponding model residuals have empirical distributions that are reasonably well approximated by Gaussian distributions, albeit there are some deviations in the tail regions. The discrepancies are most pronounced for the between-sectors blocks, as can be seen in the right columns in both panels of Figure (ref). The Q-Q plots for the variables in the Block specification appear to approximate Gaussian distributions more closely than those in the Full specification. This is not entirely unexpected because the elements of $\check{y}_{t}$ and $\check{v}_{t}$ are effectively defined as averages over elements of $y_{t}$ and $\tilde{v}_{t}$. The results for $\tilde{v}_{t}$ and $\check{v}_{t}$ provide some justification for adopting a Gaussian specification in the measurement equation.
In addition to Q-Q plots we compute the skewness and excess kurtosis for the variables in the Full specification, $y_{j,t}$ and $\tilde{v}_{j,t}$, $j=1,\ldots,36$, and present these with box-and-whiskers plots in Figure (ref). The boxes cover the inter-quartile range, whiskers the observations that are no more than 3/2 times the interquartile range away from the edge of a box, and circles represent observations outside the whiskers (outliers). The variables have, with one exception, a level of skewness that is fairly close to zero, while the excess kurtosis is about 0.5 for most variables and a handful of variables have excess kurtosis larger than one. These are closer to Gaussian moments than the skewness and kurtosis statistics reported in ABDE:2001 for (log) realized variances.
In this section, we compare the three specifications of the Multivariate Realized GARCH model with several benchmark models based on their out-of-sample performance. Our objective is to compare the different specifications for the correlation matrix out-of-sample and compare the Multivariate Realized GARCH model with natural and suitable benchmark models.
In order not to confound the comparisons with features that relate to other parts of the model, the first-stage estimation will be identical for all model in the comparisons. Specifically, we estimate the univariate Realized GARCH models, ((ref))-((ref)), for each return series, such that $h_{i,t}$ and $z_{i,t}$, $i=1,\ldots,n$, $t=1,\ldots,T$, are common for all model-specifications.
We adopt the Constant Conditional Correlation (CCC) model by Bollerslev:1990 and the Dynamic Conditional Correlation (DCC) model by Engle2002 as benchmark models for $C_{t}$. The former has (as its name suggests) a constant conditional correlation matrix and the latter uses GARCH-type dynamic for updating $C_{t}$. We label these models as CCC$^{+}$ and DCC$^{+}$, respectively, because they are enhanced CCC and DCC models that utilize realized measures of volatility for modeling the univariate conditional variances. The key features of the three types of models are summarized in Table (ref).
We consider three specifications for the correlation matrix: Equi, Block, and Full, for each model type: MRG, CCC, and DCC. The DCC$^{+}$ with equicorrelation is similar to the DECO model by EngleKelly:2012, with the key difference being that DCC$^{+}$ employs Realized GARCH models for each of the nine return series. Similarly, DCC$^{+}$ with block correlation can be viewed as an enhanced version of Block-DECO by EngleKelly:2012.\setlength\extrarowheight{8pt}
\setlength\extrarowheight{0pt}
Estimation of the CCC$^{+}$ with constant correlations simply amounts to maximum likelihood estimation of the correlation matrix for $z_{t}$. For the Full specification this is simply the sample correlation matrix, and the estimation of Equi and Block specifications is simple using the results in ArchakovHansen:CanonicalBlockMatrix. The DCC model with equicorrelation and block-correlation matrices were studied in EngleKelly:2012, and we estimate Equi and Block variants of the DCC$^{+}$ models using the method described in EngleKelly:2012. The dynamic correlation part of the DCC$^{+}$-Full model is estimated as a standard DCC model, see Aielli:2013.
The estimated CCC$^{+}$, DCC$^{+}$, and MRG models are evaluated and compared out-of-sample using log-likelihood criteria and in terms of performance for portfolio construction with a minimum-variance objective.
Of the 19 years of daily data, we use the last nine years (2,226 trading days), from January 3rd, 2012 to December 31st, 2020, for out-of-sample evaluation (testing sample). We estimate the models using (in-sample) data from January 2nd, 2002, to December 30th, 2011 (2,518 trading days), and evaluate and compare the estimated models with data from 2012. By the end of each calendar year, we update the model estimates, such that model evaluation is based on the most recent ten calendar years of data. For instance, out-of-sample comparisons during 2013 are based on models that were estimated with data from the calendar years, 2003 to 2012, and out-of-sample comparisons during 2020 are based on model-specifications that were estimated with data from the ten years from 2010 to 2019.
Multivariate GARCH models aim to describe the conditional distribution of the vector of returns, which can be quantified with the log-likelihood function for the vectors of returns. So, one natural way to compare the different models and specifications is in terms of their log-likelihoods for returns, $\ell_{r}$. In this subsection, we evaluate and compare the specifications in terms of their in-sample and out-of-sample log-likelihoods for returns, $\ell_{r,t}(\hat{\theta})$. Here $\hat{\theta}$ denotes the parameter estimates obtained from in-sample data. This comparison amounts to one-day-ahead density forecasting of the return vector with the predictive log-likelihood as a gain function, see Amisano2007, GewekeAmisano:2010, and references therein.
We compute the average in-sample and average out-of-sample return log-likelihoods, for each of the nine model-specifications. These are presented in Figure (ref) with a bar chart that reports the average log-likelihoods relative to the simplest model-specification, CCC$^{+}$-Equi, which has the smallest log-likelihood, both in-sample and out-of-sample. Good in-sample fit is no guarantee of good out-of-sample fit, but in this application we find that the relative rankings of model-specifications is largely preserved. The best model-specification in-sample, MRG-Full, is also the best model-specification out-of-sample. The three model-specifications with the worst in-sample fit all use equicorrelation structures, and these are also the three model-specifications with the worst out-of-sample fit. This is evidence that equicorrelation structures are too restrictive in this application. The block specifications are substantially better, but the full specifications have the best out-of-sample performance for all three model types. The dynamic correlation models, DCC$^{+}$ and MRG, clearly dominate the static CCC$^{+}$ model, with the MRG model having the best performance across all specifications for $C_{t}$. Once again, these results are found to be true in-sample as well as out-of-sample.
Table (ref) has additional details about the model comparisons and the statistical significance of these. The first two rows of Table (ref) report the numerical values for the bar plot in Figure (ref), and the subsequent rows have the analogous out-of-sample values for each out-of-sample calendar year. Below each of the out-of-sample statistics, we report the corresponding model confidence set (MCS) $p$-values, see HansenLundeNasonMCS. A small MCS $p$-value is evidence that the model-specification is inferior to other model-specifications in the comparison, and we use a bold font to identify the model-specifications that are included in the 95% MCS.
The best out-of-sample performance is achieved with MRG-Full, and all alternative model-specifications, except for MRG-Block, are significantly worse over the entire out-of-sample period. The MCSs for the calendar year tend to include a larger number of model-specifications, which is to expected because the shorter samples offers less information to discriminate between the competitors. MRG-Full is the only model-specification that is included in every MCS. We note that DCC$^{+}$-Block performs particularly well in 2015, but this could be a statistical anomaly, since this model-specification is inferior to MRG-Full in all other years.
The last calendar year, 2020, was a special year with the turbulence surrounding the outbreak of COVID-19 and the subsequent stock market rally. This is also the year where heterogeneous correlations were most beneficial, as evident by the large gap between equicorrelation specifications and more flexible specifications. This likely reflects the heterogeneous impact that the pandemic had on different sectors.
Next, we compare the model-specifications in terms of their ability to produce a low-variance portfolio out-of-sample. At each point in time, we compute the one-period-ahead optimal minimum-variance portfolio weights, as defined by the predicted covariance structure for each model-specification.
Let $H_{(m),t}$ denote the conditional covariance matrix for $r_{t}$, as predicted by the $m$-th model-specification at time $t-1$. We deduce the corresponding global minimum-variance (GMV) portfolio by solving \[
\] where $\iota=(1,\ldots,1)^{\prime}$ is a $n$-dimensional vector of ones. In the absence of leverage constraints, such as no-shortening constraints, the well-known solution is: \[ \omega_{(m),t}^{*}=\frac{H_{(m),t}^{-1}\iota}{\iota^{\prime}H_{(m),t}^{-1}\iota}, \] and the resulting portfolio returns are given by \[ R_{(m),t}^{\mathrm{mv}}=\omega_{(m),t}^{*\prime}r_{t}. \] Different model-specifications yield different $H_{(m),t}$, $m=1,\ldots,9$, and the resulting portfolio weights, returns, and returns-variances, will therefore be different. To this comparison, we also add the simple equal-weighted portfolio as another benchmark portfolio. The returns of the equal-weighted portfolio are given by \[ R_{t}^{\mathrm{ew}}=\frac{1}{n}\iota^{\prime}r_{t}, \] such that each asset is weighted by $1/n$, where $n=9$ in this application.
We compare the empirical variances of the ten portfolios, both in-sample and out-of-sample, and we report the results in units of annualized volatilities in Table (ref). We also report the analogous results for each of the nine years in the out-of-sample period.
A result that stands out from Table (ref), is that the equal-weighted portfolio is inferior to all nine model-specifications. This is not too surprising, because the equal-weighted portfolio does not utilize information about the covariance structure. The model-based portfolios benefit from modeling $H_{t}$, which reduce the portfolio variance by as much as 50%, which translates to a reduction in annualized volatility by a factor of about $\sqrt{2}$.
The MRG-based specifications have the smallest variances, with MRG-Block being the best model-specification in this application -- both in-sample and out-of-sample. MRG-Block reduces the annualized volatility by 50 to 150 basis points, relative to other model-based portfolios.\footnote{In practice, one would also want to account for portfolio turnover, because transaction costs may offset the gain from the reduction in variance.} The parsimonious structure of MRG-Block helps reduce the portfolio variance and MRG-Block is the only model-specification that is included in all MCSs. The MCSs for the individual calendar years are less informative, because there is not enough information to discriminate between competitors, with the exception of the equal-weighted portfolio, which is not included in any MCS.
Figure (ref) is a graphical illustration of the in-sample and out-of-sample performance of the nine model-specifications. The equal-weighted portfolio cannot be seen because its very inferior performance places it outside the range shown in Figure (ref). It is easy to see that the CCC$^{+}$ specifications (green dots) have the worst performance among the model-based portfolios. Interestingly, the most flexible specification, Full, have the worst out-of-sample portfolio performance for all types of models, which is evidence that the most flexible correlation specification, Full, suffers from overfitting in this application.
The different outcomes in the two out-of-sample comparisons highlight at important feature of model diagnostics. The “best” model-specification depends on the intended use of the model, and the empirical objective should be taken into account when applying model diagnostic tests, see HansenDumitrescu:2022. In the first comparison, we saw that the Full specification has better in-sample and out-of-sample performance in terms of the predicted log-likelihood. Whereas in the second comparison, the Block specification was significantly better for constructing minimum-variance portfolios.
In this paper, we have introduced a novel framework for multivariate GARCH modeling. Our key methodological contribution is the dynamic modeling of the conditional correlation matrix $C_{t}\in\mathbb{R}^{n\times n}$ using a vector parametrization, $\gamma_{t}=\gamma(C_{t})$. Importantly, this approach facilitates simple linear factor models of the correlation structure, while utilizing realized measures of correlations. In its most general form, the model does not impose any restrictions on $C_{t}$, while it is guaranteed to produce a valid positive definite correlation matrix. In many situations, it will be desirable to impose additional structure, especially in high-dimensional system, and the factor model for $\gamma_{t}$ serves this purpose by reducing the number of free parameters and latent variables.
Imposing block correlation structures is one way to reduce the number of free parameters. We have shown that a block structure is equivalent to a simple linear factor model, $\gamma_{t}=A\zeta_{t}\in\mathbb{R}^{d}$, where $\zeta_{t}\in\mathbb{R}^{r}$ with $r<d$ and $A$ is a matrix defined by the block structure in $C_{t}$. While the block structure is a useful example, the linear factor structure is not specific to block correlation specifications. Other choices for $A$ can be entertained, and it will be interesting to explore data-dependent choices for $A$ and non-linear factor models, $\gamma_{t}=\varrho(\zeta_{t})$. These are topics we leave for future research.
The MRG model can be seen as the natural generalization of HansenLundeVoev:2014 to higher dimensions, because $\gamma(C)$ is identical to the Fisher transformed correlation when $n=2$, and the Fisher transformation was the transformation used in the bivariate structure proposed in HansenLundeVoev:2014.
We have applied the MRG model to nine assets from three sectors, and we used three different specifications for the correlation matrix. We compared the MRG model to CCC and DCC style models and found the MRG model to be superior in terms of the predicted likelihood of returns and in terms of portfolio construction with a minimum variance objective. For all three types of models, MRG, CCC, and DCC, we also compared different structures on the conditional correlation matrix, Equi, Block, and Full. The equicorrelation structure was found to be too restrictive, with a performance that was uniformly inferior to more flexible structures. For the portfolio construction problem, the sector-based block correlation structure had the best out-of-sample performance, and a close second to the Full specification in terms of predictive log-likelihood. Conveniently, the Block specification, with a fixed number of blocks, is easy to scale to high dimensional applications.
The structure we have developed is also very beneficial for applications of the Block-DECO model by EngleKelly:2012. Specifically, the simplified expressions for the inverse block correlation matrix and its determinant by ArchakovHansen:CanonicalBlockMatrix make it straightforward to evaluate the likelihood function, thus eliminating the need to resort to alternative estimation methods when $K>2$.
We have established a wide range of attractive computational and empirical features of the MRG framework, but several theoretical properties remain unresolved. It would, for instance, be desirable to establish conditions that ensure stationarity, ergodicity, and the existence of moments for time series whose data generating process is MRG. Similarly, it would also be desirable to establish sufficient conditions that ensure the likelihood-based estimators are asymptotically normally distributed. Some results are available for univariate Realized GARCH model, see HansenHuangShek:2012. Moreover, BougerolPicard:1992, CarrascoChen:2002, JensenRahbek:2004, StraumannMikosch:06, MeitzSaikkonen:2008, Kristensen:2009, and FrancqZakoian2019, offer ideas for establishing these results for the MRG model in future research.