High-Dimensional Sparse Multivariate Stochastic Volatility Models

High-dimensional Sparse Multivariate Stochastic Volatility Models

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

} \fi

{\it Keywords:} Forecasting; Multivariate Stochastic Volatility; \textcolor{black}{Penalized} M-estimation.

{\flushleft {\bf MSC Codes:}} 62F12, 62P20. {\flushleft {\bf JEL Classification:}} C13, C32.

\spacingset{1.5}

Introduction

Over \textcolor{black}{the} past decades, various covariance models have been developed for describing dynamic structures for multivariate economic and financial time series. Within the Multivariate GARCH (MGARCH) family, the dynamic conditional correlation (DCC) model of Engle (2002) and Tse and Tsui (2002), the BEKK model of Baba et al. (1985) and Engle and Kroner (1995), and their variants are commonly used: see the survey of Bauwens, Laurent, and Rombouts (2006), for instance. {As for} the multivariate stochastic volatility (MSV) family, the MSV model of Harvey, Ruiz, and Shephard (1994) was extended, among others, by the factor model of Chib, Nardari, and Shephard (2006) and the dynamic correlation model of Asai and McAleer (2009b): see Ghysels, Harvey, and Renault (1996), Asai, McAleer, and Yu (2006), and Chib, Omori, and Asai (2009) for various univariate and multivariate SV models. Based on a thorough empirical analysis, Chib, Nardari, and Shephard (2006) highlighted that the MSV models usually outperform MGARCH-based models in terms of out-of-sample forecasts.

{Several} methods for estimating \textcolor{black}{the} MSV models have been developed. In their seminal work, Harvey, Ruiz, and Shephard (1994) derived a state space form based on the vector of the logarithm of squared returns. Based on this state space setting, they performed a Kalman-based filtering technique to evaluate and optimize the quasi log likelihood function. In the recent literature, a commonly used method is the Bayesian Markov Chain Monte Carlo \textcolor{black}{(MCMC)}, as described, \textcolor{black}{for example}, in Chib, Omori, and Asai (2009) and Kastner, Fr{$\ddot{\mbox{u}}$}hwirth-Schnatter, and Lopes (2017), among others. An alternative estimation approach is the Monte Carlo Likelihood (MCL) method suggested by Durbin and Koopman (1997, 2001) and applied by Asai, Caporin, and McAleer (2015) and Asai and McAleer (2009a). However, empirical applications in the literature are typically limited to low-dimensional random vectors when \textcolor{black}{methods based on} MCMC or MCL are performed, due to the severe costs in terms of computations, or the intricate choice of suitable priors (for the MCMC case), among others. In the same vein, MGARCH specifications also suffer from the so-called “curse of dimensionality” since the complexity is of order $O(p^2)$ \textcolor{black}{in general}, where $p$ corresponds to the problem dimension, as the specification of a general multivariate dynamic model often induces an explosion of the number of free parameters. Moreover, tricky conditions are required \textcolor{black}{for} the model parameters to satisfy the positive-definiteness of the variance-covariance process.

Another key hurdle of the aforementioned methods is the high non-linearity of the models, which requires the use of likelihood-based estimation techniques. Therefore, strongly reduced versions of such multivariate models are most often considered as soon as $p$ is larger than four or five. {The factor-model-based approach may be a solution to shrink the number of parameters. In particular,} Kastner, Fr{$\ddot{\mbox{u}}$}hwirth-Schnatter, and Lopes (2017) considered factors in their stochastic volatility {framework and provided a joint specification of a large number of covarying time series using a small number of latent factors}. However, this factor-based method requires the identification of the corresponding factors together with the treatment of the rotational \textcolor{black}{indeterminacy} inherent to factor models.

The objective of this \textcolor{black}{study} consists in \textcolor{black}{modeling} high-dimensional variance-covariance matrices within the \textcolor{black}{MSV} framework in a flexible manner and breaking the curse of dimensionality without relying on standard \textcolor{black}{procedures based on} MCMC or MCL. To do so, we introduce a vector autoregressive and moving-average (VARMA) representation for the MSV model {in the same spirit as} Harvey, Ruiz, and Shephard (1994) and apply an OLS-based two-step estimation approach extending the idea of Hannan and Rissanen (1982) and Hannan and Kavalieris (1984). {More precisely, \textcolor{black}{as} a first step, we carry out an OLS estimation of a large dimensional VAR model with a sufficiently large number of lags to approximate the VARMA model. For the purpose of parsimony and to avoid over-fitting, we enforce the nullity of potentially numerous model coefficients using a penalization procedure on the model coefficients}. Our study shares a similar spirit with Poignard and Fermanian (\textcolor{black}{2021}), who provided a framework for high-dimensional variance-covariance within the MGARCH family: they derived some parameterizations to directly generate positive-definite covariance matrices based on multivariate ARCH processes allowing for a linear representation with respect to the parameters. \textcolor{black}{However}, our work differs from theirs in two main respects: our analysis lies within the \textcolor{black}{MSV} family; we consider a general penalization framework {for efficient estimation}, which includes a broad range of potentially non-convex penalty functions.

The main contributions of our method are as follows: using a {penalized} OLS framework, we can directly generate positive-definite variance-covariance matrices without relying on \textcolor{black}{methods like} MCMC/MCL and manage high-dimensional matrix processes; the large sample properties of the two-step estimator are provided; {in particular, we prove the oracle property of the first step estimator with a diverging dimension in the sense of Fan and Li (2001), which ensures the correct identification of the underlying set of nonzero coefficients.}

The remainder of the paper is \textcolor{black}{organized} as follows. In Section (ref), we describe the framework and the new forecasting procedure based on a \textcolor{black}{penalized} OLS estimation framework. Section (ref) contains the large sample properties of the \textcolor{black}{penalized} two-step OLS estimator. Section (ref) reports simulation-based experiment results for in-sample estimates of covariance matrices together with out-of-sample forecasting results based on real financial portfolios. Finally, Section (ref) concludes the paper. All proofs and intermediary results are in the Appendix.

Notations. Throughout this paper, we denote the cardinality of a set $E$ by $\text{card}(E)$. For a vector $\mbox{\boldmath$v$} \in {\mathbb R}^d$, the $\ell_p$ norm is $\|\mbox{\boldmath$v$}\|_p = \big(\sum^p_{k=1} |\mbox{\boldmath$v$}_k|^p \big)^{1/p}$ for $p > 0$, and $\|\mbox{\boldmath$v$}\|_{\infty} = \underset{i}{\max}|\mbox{\boldmath$v$}_i|$. Let the subset \textcolor{black}{be} ${\mathcal A} \subseteq \{1,\cdots,d\}$; then, $\mbox{\boldmath$v$}_{{\mathcal A}} \in {\mathbb R}^{\text{card}({\mathcal A})}$ is the vector $\mbox{\boldmath$v$}$ restricted to ${\mathcal A}$. ${\mathcal M}_{m \times n}({\mathbb R})$ denotes the space of $m \times n$ matrices with coefficients in ${\mathbb R}$. For a matrix $A$, $\|A\|_{\textcolor{black}{F}}$ is the \textcolor{black}{Frobenius} norm. We write $A^\top$ (resp. $\mbox{\boldmath$v$}^\top$) to denote the transpose of the matrix $A$ (resp. the vector $\mbox{\boldmath$v$}$). We write $\text{vec}(A)$ to denote the vectorization operator that stacks the columns of $A$ on top of one another into a vector. We denote by $\text{vech}(A)$ the $p(p+1)/2$ vector that stacks the columns of the lower triangular part of the square and symmetric matrix $A$. \textcolor{black}{$\lambda_{\min}(A)$ (resp. $\lambda_{\max}(A)$) denotes the minimum (resp. maximum) eigenvalue of $A$.} \textcolor{black}{We write $\mbox{\normalfonttr}(A)$ to denote the trace of the square matrix $A$}. The $I_p$ matrix is the $p$-dimensional identity matrix. For a function $f: {\mathbb R}^d \rightarrow {\mathbb R}$, we denote \textcolor{black}{the gradient or subgradient of $f$ by $\nabla f$ and the Hessian of $f$ by $\nabla^2 f$}. We denote by $(\nabla^2 f)_{{\mathcal A} {\mathcal A}}$ the Hessian of $f$ restricted to the block ${\mathcal A}$. We write ${\mathcal A}^c$ to denote the complement of the set ${\mathcal A}$.

Penalized OLS framework for MSV

Framework

We consider a $p$-dimensional vectorial stochastic process $(y_t)_{t=1,\cdots,T}$ \textcolor{black}{and denote the vector of its model parameters by $\theta$. We then consider an MSV} decomposition given as

where $\Gamma$ is a $p \times p$ correlation matrix, $\varepsilon_t = (\varepsilon_{1t},\ldots,\varepsilon_{pt})^\top$ is a $p \times 1$ random vector, which is independently and identically distributed (i.i.d.), centered with variance-covariance $\Gamma$, $h_t = (h_{1t},\ldots,h_{pt})^\top$ is a $p \times 1$ vector of log-volatility, $D_t = \mbox{diag}\big(\exp( h_{1t}/2), \ldots, \exp( h_{pt}/2) \big)$ is a diagonal matrix of volatility, $\mu =(\mu_1,\ldots,\mu_p)^\top$ is a $p \times 1$ vector, $\Phi$ is a $p \times p$ matrix, and $\Sigma_{\eta}$ is a $p \times p$ covariance matrix of $\eta_t$. The MSV model ((ref)) and ((ref)) reduces to the MSV model of Harvey, Ruiz, and Shepard (1994) when $\Phi$ is diagonal and $\varepsilon_{it}$ follows a $t$ distribution.

\textcolor{black}{Subsequently}, we define $\textcolor{black}{y_t^{\ell}} = (\log( y_{1t}^2),\ldots, \log( y_{pt}^2))^\top$. Following Harvey, Ruiz, and Shepard (1994), the MSV model can be formulated as a state space model:

where $c=(c_1,\ldots,c_p)^\top$, $\zeta_t = (\zeta_{1t},\ldots,\zeta_{pt})^\top$, and $\alpha_t = h_t -\mu$ with $c_i = \mu_i + {\mathbb E}[\log (\varepsilon_{it}^2)]$ and $\zeta_{it} = \log( \varepsilon_{it}^2) - {\mathbb E}[\log (\varepsilon_{it}^2)]$. Assuming a $t$ distribution for $\varepsilon_{it}$, Harvey, Ruiz, and Shepard (1994) specified the covariance matrix of $\zeta_t$ as $\Sigma_{\zeta}$. Note that ${\mathbb E}[\zeta_t]=0$ by definition. Based on the state space form, these authors suggested a quasi-maximum likelihood estimation of the MSV model using the Kalman filter. Alternative methods were proposed such as the Bayesian MCMC technique of Chib, Nardari, and Shephard (2006) \textcolor{black}{and} the Monte Carlo Likelihood (MCL) method of Durbin and Koopman (1997, 2001). A significant drawback of these methods is the computational cost and, thus, the curse of dimensionality: most of the applications are restricted to small vector sizes and/or reduced forms are fostered.

In this paper, we aim \textcolor{black}{to tackle} this issue for \textcolor{black}{the} MSV models using a penalized OLS estimation method. Although the MSV model ((ref)) and ((ref)) might be a basic model, the following advantages with respect to \textcolor{black}{the} MGARCH models can be highlighted: (i) relatively stable estimates and forecasts for variance-covariance matrices; (ii) simpler restrictions for stationarity conditions; \textcolor{black}{and} (iii) no intricate matrix parameterization and/or parameter restrictions to generate positive-definite matrices. \textcolor{black}{Regarding (i), see the theoretical comparison of Taylor (1994, Section 5) and the empirical results of Dan{\'i}elsson (1998) and Ding and Vo (2012), for instance. As for (ii) and (iii), see Bauwens, Laurent, and Rombouts (2006) and Chib, Omori, and Asai (2009) for \textcolor{black}{the} MGARCH and MSV models, respectively.}

\textcolor{black}{ As in Harvey, Ruiz, and Shephard (1994) and Kim, Shephard, and Chib (2002), we consider the log of squared returns. Harvey, Ruiz, and Shephard (1994) suggested the quasi-maximum likelihood (QML) estimation based on the Kalman filter, by treating the distribution of $\log \chi^2(1)$ as a normal distribution. Ruiz (1994) analyzed the asymptotic properties of the QML estimator for the univariate case. Since their QML estimation depends on the numerical optimization algorithm, Shephard (1993) and So, Li, and Lam (1997) developed simulated and standard expectation-maximization algorithms, respectively. However, the inefficiency of the QML estimator comes from the fact that $\log \chi^2(1)$ is highly right-skewed. To fix this issue, Kim, Shephard, and Chib (2002) approximated $\log \chi^2(1)$ by a mixture of normal distributions to carry out a Bayesian MCMC estimation. Instead, our estimation procedure improves the efficiency \textcolor{black}{using} the \textcolor{black}{penalized} OLS regression, as previously detailed.}

Our proposed approach

\textcolor{black}{In this section, we propose a new procedure for estimating high-dimensional stochastic volatility models. Our approach starts from the measurement equation ((ref)). Instead of the state space form, we derive the VARMA representation of $\textcolor{black}{y_t^{\ell}}$ \textcolor{black}{to} apply the ideas of Hannan and Rissanen (1982) and Hannan and Kavalieris (1984) under the framework of a penalized OLS estimation. Our approach consists of four steps and can be summarized as \textcolor{black}{follows:}}

\textcolor{black}{Let us now \textcolor{black}{detail} this four step procedure.} Since $\textcolor{black}{y_t^{\ell}}$ is the sum of a VAR(1) process and an i.i.d. noise by ((ref)), the discussion of Granger and Morris (1976) suggests that $\textcolor{black}{y_t^{\ell}}$ has a VARMA(1,1) representation. By equations ((ref)) and ((ref)), we obtain \[ \textcolor{black}{y_t^{\ell}} = (I-\Phi)c + \Phi \textcolor{black}{y_{t-1}^{\ell}} + (\zeta_t + \eta_{t-1}) - \Phi \zeta_{t-1}, \] which {can alternatively be written as

with $(u_t)$ \textcolor{black}{being} a $p$-dimensional white noise vector with moments ${\mathbb E}[u_t]=0, \; \text{Var}(u_t) = \Sigma_u, \; {\mathbb E}[u_t u_s^\top] = \mathbf{0} \mbox{ for } t \neq s$,} where $\Xi$ and $\Sigma_u$ are obtained by matching moments of $w_t = (\zeta_t + \eta_{t-1}) - \Phi \zeta_{t-1}$ and $w_t^* = u_t + \Xi u_{t-1}$. Using ${\mathbb E}[w_t w_t^\top]={\mathbb E}[w_t^* w_t^{*\top}]$ and ${\mathbb E}[w_t w_{t-1}^\top]={\mathbb E}[w_t^* w_{t-1}^{*\top}]$, the relationship between $(\Xi, \Sigma_u)$ and other parameters is given as follows:

{$\Sigma_{\eta}$ and $\Sigma_{\zeta}$ can be deduced from $\Phi$, $\Xi$, and $\Sigma_u$ based on equations ((ref)) and ((ref)). \textcolor{black}{Let $(x_t)$ denote the mean-subtracted process $x_t=y_t^{\ell} - {\mathbb E}[y_t^{\ell}]$.} Assuming a stable and invertible model, $(x_t)$ has an AR($\infty$) representation:}

\textcolor{black}{with $u_t$ defined in equation ((ref)).} Based on a penalized OLS estimation, we can obtain \textcolor{black}{an approximation of $u_t$ in the first step, denoted as $\widehat{u}^{(m)}_t$. The latter approximation depends on $m$: \textcolor{black}{we} empirically need to specify $m$ sufficiently large as a surrogate of $\infty$ in the summation in ((ref)). Thus, for the sake of parsimony and to avoid the over-fitting issue, we assume sparsity among the $\Psi_i$\textcolor{black}{s}.} In the second step, we calculate the OLS estimator of $(\widehat{c}^*,\widehat{\Phi},\widehat{\Xi})$ by regressing $x_t$ on a constant, $x_{t-1}$, and \textcolor{black}{$\widehat{u}^{(m)}_{t-1}$}. For the third step, we start from the decomposition of the unconditional variance-covariance matrix of $x_t$, which is given by

where $\Sigma_x = {\mathbb E}[\textcolor{black}{x_t x_t^\top}]$, $\Sigma_{\zeta} = {\mathbb E}[\zeta_t \zeta_t^\top]$, and $\Sigma_{\alpha} = {\mathbb E}[\alpha_t \alpha_t^\top]$ with \[ \mbox{\normalfontvec} (\Sigma_{\alpha}) = [ I_{p^2} - (\Phi \otimes \Phi)]^{-1} \mbox{\normalfontvec} (\Sigma_{\eta}). \] {Denoting the sample covariance matrix of $x_t$ and \textcolor{black}{$\widehat{u}^{(m)}_t$} by $S_x$ and $S_{\widehat{u}^{\textcolor{black}{(m)}}}$, respectively, we obtain an estimator of $\Sigma_{\zeta}$ as} \[ S_{\zeta} = - \frac{1}{2} \left[ \widehat{\Phi}^{-1} \widehat{\Xi} S_{\widehat{u}^{\textcolor{black}{(m)}}} + S_{\widehat{u}^{\textcolor{black}{(m)}}} \widehat{\Xi}^\top \widehat{\Phi}^{\top -1} \right], \] by \textcolor{black}{the sample analogous of the mean of $\Sigma_{\zeta}$ obtained by equation ((ref)) and its transpose.} As there is no guarantee for $S_{\zeta}$ and $S_x-S_{\zeta}$ to be \textcolor{black}{positive-definite} by the approach, we consider \textcolor{black}{ad hoc} estimators for $\Sigma_{\zeta}$ and $\Sigma_{\alpha}$ based on decomposition ((ref)). Finally, in the fourth step, we estimate $\Gamma$ by a correlation matrix of $y_t$. \\ \textcolor{black}{To summarize, our procedure can be broken down as follows:}

When the tuning parameter $\lambda_T$ shrinks to zero, Steps 1 and 2 reduce to the standard OLS estimation for low-dimensional VARMA models considered by Hannan and Rissanen (1982) and Hannan and Kavalieris (1984). Step 1 corresponds to a multivariate version of the AR($\infty$) representation of a log-GARCH model. Although Harvey, Ruiz, and Shephard (1994) applied the Kalman filter, its computational cost is non-negligible for large $p$, since the cost {evolves according to $O(T p^2)$} for storing covariance matrices of a $p \times 1$ state vector for all $t=1,\ldots,T$. For the estimators in Step 3, we may improve them by considering moment-matching methods using equations ((ref)), ((ref)), and $(\ref{eq:deco})$ with restrictions on the positive-definiteness of the estimators of $\Sigma_{\zeta}$ and $\Sigma_{\alpha}$. However, we use the above fast {and efficient} method described in Step 3 without the need of a numerical optimization procedure. Finally, the fourth step can easily be adapted to a sparse correlation matrix setting, especially when the size $p/T$ is not negligible.

We now introduce our setting for generating the volatility process. For a low-dimensional case, we can calculate the minimum mean square linear estimator (MMSLE) of $\alpha_t$ based on the full sample \textcolor{black}{$\mbox{\boldmath$ y$}^\ell=(y_1^{\ell \top},\ldots,y_T^{\ell \top})^\top \in {\mathbb R}^{pT}$} by the state space smoothing algorithm. In the high-dimensional case, we consider the multivariate version of Harvey (1998)\textcolor{black}{'s approach} with the vector form of ((ref)) as follows: \[ \textcolor{black}{\mbox{\boldmath$ y$}^{\ell}} = \mbox{\boldmath$ c$}^\dagger + \mbox{\boldmath$ \alpha$} + \mbox{\boldmath$ \zeta$}, \] where $\mbox{\boldmath$ c$}^\dagger = (\iota_T \otimes c) \in {\mathbb R}^{pT}$, $\mbox{\boldmath$ \alpha$}=(\alpha_1^\top,\ldots,\alpha_T^\top)^\top \in {\mathbb R}^{pT}$, and $\mbox{\boldmath$ \zeta$}=(\zeta_1^\top,\ldots,\zeta_T^\top)^\top\in {\mathbb R}^{pT}$. By the model structure, the covariance matrix of $\mbox{\boldmath$ x$}$ is given by \[ V_x = V_{\alpha} + V_{\zeta}, \] where \[ V_{\alpha} = \left(

\right), \] and $V_{\zeta} = (I_T \otimes \Sigma_{\zeta})$. Then\textcolor{black}{,} the MMSLE \textcolor{black}{can be calculated as follows:}

As in Harvey (1998), the covariance matrix is deduced from the relationship $H_t = D_t \Gamma D_t$ such that the sample variance of the standardized variable of $y_{it}$ equals to one. We consider the estimator as $\widetilde{H}_t = \widetilde{D}_t \widehat{\Gamma} \widetilde{D}_t$, where

for $i=1,\ldots,p$. The standardized variables are defined as $\widetilde{z}_{it} = y_{it} / \widetilde{d}_{it}$, which implies $T^{-1} \sum_{t=1}^T \widetilde{z}_{it}^2 = 1$ by definition. We call our proposed parameterization “penalized OLS-MSV”.

Volatility forecasting

We now provide the forecasts for variance-covariance based on our proposed method. The MMSLE for the $l$th-step-ahead forecast of $\alpha_T$ is given by

where $R_l = \left[ \Phi^{T+l-1} \Sigma_{\alpha} \;\; \Phi^{T+l-2} \Sigma_{\alpha} \; \cdots \; \Phi^{l} \Sigma_{\alpha} \right]$. The $l$th-step-ahead forecast of the covariance matrix is given by

where for $i=1,\ldots,p$, \[ \widehat{D}_{T+l} = \mbox{diag} \big(\widehat{d}_{1,T+l},\ldots, \widehat{d}_{p,T+l} \big), \quad \widehat{d}_{i,T+l} = \bar{d}_i \exp \left( \widehat{x}_{i,T+l}/2 \right). \] \textcolor{black}{ By the structure of $R_l$ and $V_x$, the inconsistency on the off-diagonal elements of the third-step estimator may affect the forecasts. We assess its applicability via the Monte Carlo experiments performed in Section (ref).}

Using the covariance forecasts, we can construct time-varying portfolios for the asset returns, in which \textcolor{black}{the} portfolio weights are determined by \textcolor{black}{past information}. \textcolor{black}{To assess the relevance of the volatility models in terms of forecasts, we can obtain the value-at-risk (VaR) threshold of the portfolio return. This threshold is the negative of the 100$q$-th percentile of the portfolio return distribution, with $q$ small such as $q=0.01$, and may be used in the test procedure of Candelon et al. (2011). More precisely, the VaR threshold at time $t+1$ is given by $- \tau_q \sqrt{\widehat{h}_{t+1}^p}$, where $\widehat{h}_{t+1}^p=w_{t+1}^\top \widehat{H}_{t+1} w_{t+1}$, $w_{t+1}$ is the $p\times 1$ vector of portfolio weights at time $t+1$, and $\tau_q$ is the 100$q$-th percentile of the standard normal distribution or a historically simulated distribution of $\{ (w_t^\top y_t)/\sqrt{\widehat{h}_t^p}\}$. Alternatively,} we can compare forecasting models directly based on the minimum-variance portfolio via the test suggested by Engle and Colacito (2006).

Asymptotic properties

In this section, we provide the asymptotic properties of the penalized two-step estimator. We show that the first step estimator satisfies the oracle property for the SCAD and MCP cases and when the number of parameters diverges with the sample size. Conditional on this first step sparse estimator, we derive the conditions for consistency and asymptotic normality of the second step estimator.

First step penalized estimator $\widehat{\Psi}_{1:m}$

\textcolor{black}{In Step 1, we estimate the parameter $\theta = \text{vec}(\Psi_{1:m})$ with $d=mp^2$, the dimension that can diverge with the sample size $T$. Consequently, both dimension $d=d_T$ and parameter $\theta=\theta_T$ are indexed \textcolor{black}{hereafter} by $T$ to highlight the dependence of \textcolor{black}{$d$ and, thus, $\theta$}, with respect to $T$. More formally, we consider a sequence of parametric models ${\mathcal P}_T:=\{{\mathbb P}_{\theta_T}, \,\theta_T \in \Theta_{1,T}\}$, $\Theta_{1,T}\subset {\mathbb R}^{d_T}$. We \textcolor{black}{denote the non-penalized loss function} by ${\mathbb G}_T: {\mathbb R}^{pT} \times \Theta_{1,T} \rightarrow {\mathbb R}$: the value ${\mathbb G}_T(\underline{y};\theta_T)$ with $\theta_T \in \Theta_{1,T}$ evaluates the quality of the “fit” for the realizations of $y_t$ for every $t=1,\cdots,T$ and under ${\mathbb P}_{\theta_T}$. The loss ${\mathbb G}_T(\underline{y};\theta)$ is associated to a continuous function $\ell : {\mathbb R}^{pT} \times \Theta_{1,T} \rightarrow {\mathbb R}$ that can be written as

where $x_t$ corresponds to the vector of continuous transforms of $\log(y^2_{it})$ (vector of mean-subtracted series), $\Psi_{1:m}=(\Psi_1,\cdots,\Psi_m)\in {\mathcal M}_{p \times pm}({\mathbb R})$, and $Z_{m,t-1}=(x^\top_{t-1},\cdots,x^\top_{t-m})^\top\in {\mathbb R}^{pm}$. For every $T$, we assume there exists a unique pseudo-true parameter value $\theta_{0,T}$: for every $T$, the function $\theta_T \mapsto {\mathbb E}[\ell(y_s,s\leq t;\theta_T)]$ is uniquely minimized on $\Theta_{1,T}$ at $\theta_T = \theta_{0,T}$ and the first-order conditions are satisfied, \textcolor{black}{that is}, ${\mathbb E}[ \nabla_{\theta_T}{\mathbb G}_T(\underline{y};\theta_{0,T})] = 0$. In light of the possibly explosive number of parameters for a given $T$, $\theta_{0,T}$ is assumed sparse so that the size of the true support $k_T = \text{card}({\mathcal A}_T)$, with ${\mathcal A}_T:=\{i=1,\cdots,d_T: \theta_{0,i,T}\neq 0\}$, also diverges with $T$. To estimate the latter support, we rely on the penalty function $\textbf{pen}(\frac{\lambda_T}{T},.)$, which is assumed coordinate-separable, \textcolor{black}{that is}, $\textbf{pen}(\frac{\lambda_T}{T},\theta_{T}) = \sum^{d_T}_{i=1}\mbox{\boldmath$p$}(\frac{\lambda_T}{T},|\theta_{i,T}|)$. Then, the penalized problem becomes

} For the penalty function, we consider the convex penalty LASSO $\mbox{\boldmath$p$}(\lambda,|\theta|) = \lambda |\theta|$ of \textcolor{black}{Tibshirani (1996)} and the non-convex penalties SCAD and MCP. The SCAD of Fan and Li (2001) is defined as

where $a>2$. The MCP due to Zhang et al. \textcolor{black}{(2010) is defined} for $b>0$ as

All assumptions we relied on for the large sample analysis are reported in Section (ref) of the Appendix. In particular, the sparsity assumption states that the true parameter vector is sparse, that is\textcolor{black}{,} the cardinality of the true sparse support \textcolor{black}{${\mathcal A}_T$} is of size $k_T < d_T$. We assume stability of the VARMA(1,1) process $(x_t)$ to apply the large sample theory for stationary processes. Finally, we assume suitable regularity conditions for both the non-penalized loss and the penalty function.

We first show the existence of the penalized estimator $\widehat{\theta}_T$ for the three aforementioned penalty cases.

Our second result is dedicated to the oracle property: we show that the penalization procedure in problem ((ref)) asymptotically recovers the true underlying sparse subset ${\mathcal A}_{\textcolor{black}{T}}$ and the nonzero estimated coefficients are normally distributed. We prove the oracle property for the SCAD and MCP only: these penalty functions are non-convex, a key property that enables to relax the incoherence/irrepresentable condition and/or avoid the specification of adaptive weights. \textcolor{black}{The} incoherence/irrepresentable condition - see inequality (3) of Zou (2006) regarding the irrepresentable condition; see Loh and Wainwright (2017) regarding the incoherence condition - is necessary to prove the oracle property for the LASSO: such condition is nontrivial and difficult to empirically verify. Rather than assuming the incoherence/irrepresentable condition, Zou (2006) proposed the adaptive LASSO: stochastic weights are specified in the LASSO penalization to alter the convergence rate of the regularization parameter $\lambda_T$; such weights depend on a first step $\sqrt{T/d_T}$-consistent estimator, typically an \textcolor{black}{non-penalized} OLS estimator: the adaptive LASSO is consequently a two-step procedure. In the same vein, Poignard (2020) specified adaptive weights in the Sparse Group LASSO penalty - $\ell_1+\ell_1/\ell_2$ penalty - since the convexity of the $\ell_1$ and $\ell_1/\ell_2$ norms prevents from satisfying the oracle property. The key advantage of non-convex penalization is the relaxation of the incoherence/irrepresentable condition and avoids a two-step procedure as in the adaptive LASSO.

{

}

{The first step estimator is deduced from the truncated VAR process ((ref)): a VAR($m$) is fitted to obtain $(\widehat{u}^{(m)}_t)$ as an approximation of $(u_t)$. In this context, the specification of a diverging number of parameters is relevant to correctly approximate $(u_t)$. \textcolor{black}{When} $p$ is fixed and $m:=m_T \rightarrow \infty$, our scaling condition $d^3_T=o(T)$ for the oracle property is identical to condition (15.2.5) of Proposition 15.1. (result of Lewis and Reinsel, 1985) of L\"{u}tkepohl (2006). However, our setting does not enable to simultaneously distinguish $m_T \rightarrow \infty$ \textcolor{black}{and} $p:=p_T \rightarrow \infty$. Furthermore, we may derive a finite sample and explicit upper bound for the approximation error $\frac{1}{T}\overset{T}{\underset{t=1}{\sum}}\|\widehat{u}^{(m)}_t-(u_t+\underset{i > m}{\sum}\Psi_ix_{t-i})\|_2$ in the same spirit as in Proposition 2.4. of Wilms, Basu, Bien, and Matteson (2021), who considered an approximation for VARMA(p,q) processes and relied on a LASSO penalization. The derivation of such bound would require additional assumptions on the non-convex penalty functions - such as the $\mu$-amenable assumption as in Loh and Wainwright (2017) - and the derivation of an exponential bound over $\nabla_{\theta_{\textcolor{black}{T}}}{\mathbb G}_T(\underline{y};\theta_{\textcolor{black}{T}})$. We leave this topic for future research.}

Second step estimator $(\widehat{c}^*,\widehat{\Phi},\widehat{\Xi})$

We consider the large sample properties of the second step estimator $\widehat{\gamma} = (\widehat{c}^{*\top},\text{vec}(\widehat{\Phi})^\top,\text{vec}(\widehat{\Xi})^\top)^\top$, which is of size $d_2 = p(1+2p)$ assumed fixed. \textcolor{black}{Conditional on} $\widehat{\theta}_{\textcolor{black}{T}}$, we consider a second step loss function ${\mathbb L}_T$ from ${\mathbb R}^{pT} \times \Theta_2$ to ${\mathbb R}$ with $\Theta_2 \subset {\mathbb R}^{p(1+2p)}$ compact, and ${\mathbb L}_T(\underline{y};\widehat{\theta}_{\textcolor{black}{T}},\gamma)$ is the empirical loss associated to a continuous function $f: {\mathbb R}^{pT} \times \Theta_{1,\textcolor{black}{T}} \times \Theta_2 \rightarrow {\mathbb R}$, \textcolor{black}{that is},

where $K_{t-1}(\widehat{\theta}_{\textcolor{black}{T}}) = (1,\textcolor{black}{y_{t-1}^{\ell\top}},{\widehat{u}^{(m)\top}_{t-1}})^\top \in {\mathbb R}^{1+2p}$ and $\Gamma = (c^*,\Phi,\Xi) \in {\mathcal M}_{p\times (1+2p)}({\mathbb R})$. The dependence with respect to the first step estimator is through {$\widehat{u}^{(m)}_t$}. The problem of interest is

The function $\gamma \rightarrow {\mathbb E}[f(y_s,s\leq t; \theta_{0,\textcolor{black}{T}},\gamma)]$ is assumed to be uniquely minimized at $\gamma = \gamma_0$\textcolor{black}{; the} true parameter vector $\gamma_0 = (c^{*\top}_0,\text{vec}(\Phi_0)^\top,\text{vec}(\Xi_0)^\top)^\top$.

We now derive the asymptotic distribution of the second step estimator $\widehat{\gamma}$ \textcolor{black}{conditional on} $\widehat{\theta}_{\textcolor{black}{T}}$ \textcolor{black}{whose elements belong to ${\mathcal A}_{\textcolor{black}{T}}$}.

The third step estimator is accurate in the sense that $p^{-1}\mbox{\normalfonttr}(\widehat{\Sigma}_{\zeta})$ always takes the true value ${\mathbb E}[p^{-1}\mbox{\normalfonttr}(\widehat{\Sigma}_{\zeta})]=\pi^2/2$ under the Gaussian assumption. To improve the estimator, we can consider the structure $\Sigma_{\zeta} = \sigma_{\zeta}^2 P_{\zeta}$, where $P_{\zeta}$ is a correlation matrix and $\sigma_{\zeta}^2$ is the variance for non-Gaussian assumption. Neglecting the computational costs, we may estimate \textcolor{black}{positive-definite} $\Sigma_{\zeta}$ and $\Sigma_{\alpha}$ given $\widehat{\Phi}$ and $\widehat{\Xi}$ under the restrictions discussed below equation ((ref)).

Empirical analysis

Simulation experiment

In this section, we empirically investigate the ability of the proposed penalization method to better capture complex variance-covariance processes. We simulate the $p$-dimensional stochastic process $(\epsilon_t)$ based on two data generating processes (DGP\textcolor{black}{s}): the multivariate ARCH and the BEKK processes. For the multivariate ARCH with $q^*$ lags - M-ARCH($q^*$) in the rest of the paper - case, we consider the DGP

where $q^*$ is the number of lagged matrices being functions of $\epsilon_{t-k}$ and the $p^2 \times p^2$ square matrices $A_k$ satisfy the stationarity conditions of Theorem 2 of Boussama (2006) together with the positivity condition given by Gouri\'eroux (1997). We generate the diagonal elements of $A_k$ from a uniform distribution ${\mathcal U}([0.01,0.05])$ and the off-diagonal ones from ${\mathcal U}([-0.01,0.01])$ under the ordering constraint $\forall k \geq 2, \forall i,j, \, |A_{k,ij}| \leq |A_{k-1,ij}|$. Note that these coefficients are more constrained (i.e.\textcolor{black}{,} closer to zero) when the dimension $p$ increases. As for the matrix $\Omega$, the diagonal and off-diagonal elements are simulated from ${\mathcal U}([0.1,0.2])$ and ${\mathcal U}([-0.01,0.01])$, respectively. As for the BEKK process, the DGP is based on

where $A,B$ are $p \times p$ matrices, satisfying the stationarity constraint $\|D^+_p\{ \big(A \otimes A \big) + \big(B \otimes B \big)\} D_p\|_s < 1$, \textcolor{black}{and $D_p$ and $D^+_p$ are the duplication matrix and elimination matrix, respectively} (see Subsection 11.3 “Stationarity of VEC and BEKK Models" of Francq and Zako\"ian (2010) for the stationarity condition and Remark 11.1 for the definition of the latter matrices). The entries of $A$ and $B$ are generated from the uniform distribution ${\mathcal U}([-0.8,0.8])$. The matrix $\Omega$ is generated as in the M-ARCH($q^*$) case. Unlike the M-ARCH($q^*$) case, the BEKK dynamic includes an autoregressive component through $B$, which motivated the use of larger lags when estimating our proposed parameterizations. In both proposed dynamics, we initialize the observations $(\epsilon_{k},\cdots,\epsilon_{1})$ with centered and unit variance multivariate Gaussian distribution, where $k = q^*$ in the M-ARCH model and $k = 1$ in the BEKK \textcolor{black}{model}. \textcolor{black}{Further}, conditional on the past $k$ observations, we generate $H_t$ and, thus, $\epsilon_t$ according to a centered multivariate Gaussian distribution with variance-covariance $H_t$.

We consider the problem sizes, $p = 15, 50, 100$, and $T = 800$ observations for each of them. For the M-ARCH($q^*$)-based data generating process, we considered $q^*=2$ when $p=15$ and $q^*=1$ when $p=50, 100$. \textcolor{black}{Subsequently}, we propose to compare the true variance-covariance processes - BEKK and M-ARCH($q^*$) - and the estimated ones through our proposed MSV model and the scalar DCC together with the constant correlation model (CCC). The estimation of the DCC model is based on the classic two-step Gaussian QMLE, where the marginal conditional volatility processes are specified as GARCH(1,1) and a correlation targeting procedure is applied in the second step, providing an estimated trajectory $\widehat{H}^{\text{dcc}}_t$. The CCC is estimated thanks to a joint estimation of the GARCH(1,1) parameters and correlation parameters through a Gaussian QML, which provides an estimated process $\widehat{H}^{\text{ccc}}_t$. More details on the DCC and CCC can be found in Appendix (ref).

Regarding our proposed variance-covariance dynamic, the penalized OLS-MSV, denoted as $\widehat{H}^{\text{ols,scad}}_t$ for the SCAD OLS-MSV, $\widehat{H}^{\text{ols,mcp}}_t$ for the MCP OLS-MSV, and the non-penalized version of the OLS-MSV denoted as $\widehat{H}^{\text{ols}}_t$, we considered two settings depending on the DGP. In the M-ARCH($q^*$) case, we set the number of lags in Step 1 in ((ref)) as $m = 10$ when $p=15$ and set \textcolor{black}{it} as $m=5$ for a dimension $p=50, 100$. In the BEKK case, due to the autoregressive nature of the latter dynamic, more lags were specified: we selected $m=30$ (resp. $m=15$, resp. $m=5$) when $p=15$ (resp. $p=50$, resp. $p=100$). For both DGPs, the correlation matrix of Step 4 is estimated as the sample correlation matrix estimator. In the SCAD and MCP cases, the coefficients $a$ and $b$ are set as $3.5$ - a value close to the optimal one as in Fan and Li (2001) - and $3$\textcolor{black}{, respectively}.

We compare the true variance-covariance and the estimated variance-covariance processes through the aforementioned models. To do so, we specify a matrix distance, namely, the Frobenius norm, defined as $ || A - B ||_F := \sqrt{\mbox{\normalfonttr}((A - B)^\top(A - B))} $. We compute the previous norm for each $t$ and for $A = H_t$ and $B \in \big\{\widehat{H}^{\text{dcc}}_t,\widehat{H}^{\text{ccc}}_t,\widehat{H}^{\text{ols}}_t,\widehat{H}^{\text{ols,scad}}_t,\widehat{H}^{\text{ols,mcp}}_t \big\}$. We take the average of those quantities over $T=800$ periods of time. Since we repeat this experiment $100$ times, this provides an average gap for all those simulations.

By a cross-validation (CV) procedure - see, \textcolor{black}{for example}, Hastie \textcolor{black}{et} al. (2015, Chap. 2) - we selected the regularization parameter and emphasize that the standard CV developed for i.i.d. data \textcolor{black}{cannot} be used in our time series framework. To fix this issue, we used the hv-CV procedure devised by Racine (2000), which consists in leaving a gap between the test sample and the training sample, on both sides of the test sample.

The average difference results are reported in Table (ref) for the M-ARCH($q^*$)\textcolor{black}{-based} DGP and Table (ref) for the BEKK-based DGP. First, our proposed method provides better in-sample results in terms of accuracy compared to standard MGARCH models. The results are closer to each other in the BEKK\textcolor{black}{-based} DGP case, essentially due to the presence of an autoregressive term, which is a priori in \textcolor{black}{favor} of the DCC/CCC model. Interestingly, our results \textcolor{black}{emphasize} the gain in considering a penalized MSV, especially when the dimension increases.

Application to real data

To assess the relevance of the proposed penalized method, we propose a real data experiment, where we focus on direct out-of-sample evaluation methods, which allow for \textcolor{black}{pair-wise} comparisons. They test whether some of the variance-covariance models provide better forecasts in terms of portfolio volatility behavior. Following the methodology of Engle and Colacito (2006), we develop a mean-variance portfolio approach to test the $H_t$ forecasts. Intuitively, if a conditional covariance process is misspecified, the minimum variance portfolio should emphasize such a shortcoming, compared to other models. \textcolor{black}{Here}, consider an investor who allocates a fixed amount between $p$ stocks, according to a minimum-variance strategy and independently at each time $t$:

where $w_t$ is the $p \times 1$ vector of portfolio weights chosen at (the end of) time $t-1$, $\iota$ is a $p \times 1$ vector of $1$, and $H_t$ is the estimated conditional covariance matrix of the asset returns at time $t$. The solution of ((ref)) is given by the global minimum variance portfolio $w_t = H^{-1}_t \iota/\iota^{\top}H^{-1}_t\iota$. Engle and Colacito (2006) \textcolor{black}{showed} that the realized portfolio volatility is the smallest when the variance-covariance matrices are correctly specified. \textcolor{black}{Consequently}, if wealth is allocated using two different dynamic models $i$ and $j$, whose predicted covariance matrices are $(H^{i}_t)$ and $(H^{j}_t)$, the strategy providing the smallest portfolio variance will be considered as the best. To do so, we consider a sequence of minimum variance portfolio weights $(w_{i,t})$ and $(w_{j,t})$, depending on the model. \textcolor{black}{Further}, we consider a distance based on the difference of the squared returns of the two portfolios, defined as $u_{ij,t} = \left\{w^{\top}_{i,t} \epsilon_t\right\}^2 - \left\{w^{\top}_{j,t} \epsilon_t\right\}^2$. The portfolio variances are the same if the predicted covariance matrices are the same. Thus, we test the null hypothesis ${\mathcal H}_0: \; \mathbb{E}\left[u_{ij,t}\right] = 0$ by the Diebold and Mariano (1995) test. It consists of a least squares regression using HAC standard errors, given by $u_{ij,t} =\alpha + \epsilon_{u,t}$, $\mathbb{E}[\epsilon_{u,t}]=0$, and we test $ \mathbf{H0}: \alpha = 0$. If the mean of $u_{ij,t}$ is significantly positive (resp. negative), the forecasts given by the covariance matrices of model $j$ (resp. $i$) are preferred. \textcolor{black}{Following Engle and Colacito (2006), we compute the test statistic

with $h$ \textcolor{black}{being} the number of one-period ahead forecasts and $\sqrt{\widehat{\text{Var}}(\sqrt{h}\,\overline{u}_{ij})}$ is a \textcolor{black}{heteroscedasticity} and autocorrelation consistent estimator of the asymptotic variance of $\sqrt{h}\,\overline{u}_{ij})$. In particular, $\widehat{\text{DM}}_{ij} \overset{d}{\longrightarrow} {\mathcal N}_{{\mathbb R}}(0,1)$ under $\textbf{H0}$.}

We run the latter test to compare the scalar DCC (DCC), the Orthogonal GARCH (O-G), the BEKK (BEKK), and our OLS-MSV method (MSV) together with its penalized counterpart (denoted as MSV-SCA and MSV-MCP for the SCAD and MCP, respectively). We also consider the adaptive LASSO, denoted as MSV-AL, as an additional convex penalization technique developed by Zou (2006). In that case, we selected $\delta=3$ (resp. $\delta=4$) for the low-dimensional (resp. high-dimensional) portfolio as the power entering the stochastic weights and the first step estimator is the \textcolor{black}{non-penalized} OLS estimator, following Zou (2006). The definitions of the BEKK and O-GARCH processes are in Appendix (ref). No variance-targeting was applied in the BEKK.

We consider two different data sets. \textcolor{black}{First is} a low-dimensional portfolio of daily financial returns composed of the MSCI stock index for the following 23 countries over the period December 1998 - March 2018: Australia, Austria, Belgium, Canada, Denmark, Finland, France, Germany, Greece, Hong Kong, Ireland, Italy, Japan, Netherlands, New Zealand, Norway, Portugal, Singapore, Spain, Sweden, Switzerland, the United-Kingdom, \textcolor{black}{and} the United-States. The second portfolio corresponds to the daily stock returns of the S&P 100, where we considered firms that have been continuously included in the index from December 2015 until January 2020, excluding AbbVieInc., Dow Inc., Facebook, Inc., General Motors, Kraft Heinz, Kinder Morgan, and PayPal Holdings\textcolor{black}{. This} leaves $94$ assets and, thus, corresponds to the high-dimensional portfolio. The matrix $H_t$ in ((ref)) is deduced from the aforementioned dynamics that have been estimated on the sub-sample December 1998 - November 2015 (resp. February 2010 - January 2018), \textcolor{black}{that is 4000 (resp. $2000$) observations} for the MSCI (resp. S&P 100) portfolio. Once the latter process is estimated in-sample, out-of-sample predictions \textcolor{black}{have been} plugged into the program ((ref)) between December 2015 and March 2018 (resp. February 2018 - January 2020) for the MSCI (resp. S&P 100) portfolio. For both data sets, the MSV-based models are estimated with $m=10$ lags, \textcolor{black}{and as an alternative, with $m=20$ lags. The calibration of $m$ is set in light of the scaling $d^3_T=o(T)$ in Theorem (ref) to satisfy the oracle property, or equivalently $d_T = O(T^c)$ with $0<c<1/3$\textcolor{black}{; that is}, there exists $L>0$ a finite constant such that $d_T \leq L \, T^c$. For $p$ fixed and $T$ the in-sample size, then, $T^{c} \approx 16$ (resp. $13$) with $c=1/3-\epsilon$, $\epsilon\rightarrow 0$ for the MSCI (resp. S&P) portfolio, providing a calibration order for the number of lags. The matrix forecast comparisons are provided in Tables (ref) and (ref). First, the results \textcolor{black}{emphasize} that the proposed penalized OLS-MSV method outperforms the MGARCH-based competitors in both portfolio cases. For the low-dimensional portfolio, no clear-cut results are in favor of penalized MSV over non-penalized MSV. \textcolor{black}{However}, interestingly, fostering sparsity yields much better forecasting performances for the high-dimensional portfolio in the adaptive LASSO and SCAD (with $20$ lags), at least. Furthermore, the calibration $m=20$ for the low-dimensional portfolio does not provide better forecasts. The results are different in the high-dimensional portfolio: for each penalized MSV model, a larger $m$ with sparse estimation provides better performances; there is a gain in penalizing the MSV process over the non-penalized version of the MSV.}

The results based on the Diebold-Mariano test are limited since they are \textcolor{black}{pair-wise} comparisons. \textcolor{black}{It} is not possible to ensure that an optimal test is clearly identified. To tackle this issue, Hansen, Lunde, and Nason (2003, 2011) proposed the Model Confidence Set (MCS) method, which is a testing framework for the null hypothesis of equivalence across subsets of models. Starting with a full set of candidate models, the MCS method sequentially trims the elements of this set, thus reducing the number of viable models. To be more precise, this approach performs an iterative selection procedure testing the null hypothesis of equal forecasting ability among all models included in a set ${\mathcal M}$ (the starting set containing all candidate models) for a given loss function. The null hypothesis is $\mathbf{H0:} \; {\mathbb E}[u_{ij,t}]=0$, $i>j$ for any $i,j \in {\mathcal M}$. To test $\mathbf{H0}$, Hansen, Lunde, and Nason (2003) proposed the following two statistics:

\textcolor{black}{where} $h$ is the number of one-period ahead forecasts and $\sqrt{\widehat{\text{Var}}(\overline{u}_{ij})}$ is the bootstrap estimate of the variance of $\overline{u}_{ij}$. The $p$-values of the test statistics are obtained using a bootstrap method. For a given confidence level, if $\mathbf{H0}$ is rejected, the worst performing model is excluded from the set ${\mathcal M}$, where such a model is identified using the \textcolor{black}{following} rule:

where the variance is obtained again using a bootstrap approach. Table (ref) reports the MCS results for both MSCI and S$\&$P 100 portfolios when applied to the loss function $u_{ij,t}$. We used the statistic $t_{SQ}$ to compute the $p$-values for three confidence levels ($5\%, \; 10\%$, and $20\%$). These $p$-values inform about the included/excluded models for a given confidence level. If a $p$-value is larger than the fixed confidence level, the corresponding model is included in the MCS test of statistically equivalent models. The higher the $p$-value, the better the model is in terms of prediction accuracy. For the MSCI portfolio results in Table (ref), we can draw the following remarks: our MSV specifications are always included for any confidence level, contrary to the standard MGARCH models; among the MSV specifications, the penalized processes provide better forecasting performances. For the high-dimensional S$\&$P 100 \textcolor{black}{portfolio} in Table (ref), only the adaptive LASSO (with $m=10$) and SCAD (with $m=20$) penalized models are included in the test for all confidence levels. These results support our findings in the Diebold-Mariano test.

Conclusion

The focus of this \textcolor{black}{study} was \textcolor{black}{on} the estimation of high-dimensional MSV models. Our main contribution consisted in proposing an estimation framework that does not rely on standard MCMC/MCL methods but instead on a penalized OLS framework for state-space estimation. The corresponding large sample properties of the two-step estimator are derived. In particular, we considered a sparse first step estimator for a broad range of penalty functions when the number of parameters is diverging. We derived an explicit convergence rate of the second step estimator and its large sample distribution. The performances of our proposed method compared to standard MGARCH models are illustrated through simulated experiments together with an out-of-sample analysis for prediction accuracy, where our method clearly outperformed the competing MGARCH models. These results also emphasized the gain of penalization, which manages \textcolor{black}{the over-fitting problem}.

Various issues and extensions can be further considered. Our proposed model could be extended to accommodate a factor structure, long memory, and/or asymmetry, as discussed in Asai, McAleer, and Yu (2006) and Chib, Omori, and Asai (2009). \textcolor{black}{Besides} the factor setting considered by Chib, Nardari, and Shephard (2006) and So and Choi (2009), we can consider the rotation of the variables as in Harvey, Ruiz, and Shephard (1994) and Hafner and Preminger (2009). For the long memory property, So and Kwok (2006) extended Harvey, Ruiz, and Shephard (1994)'s model\textcolor{black}{; hence}, we can consider applying their work. For including asymmetric effects, we may extend Harvey and Shephard (1996)\textcolor{black}{'s approach}; see Asai and McAleer (2009a) for the multivariate case. Another direction \textcolor{black}{would include} modeling directly the variance-covariance matrix $H_t$ without relying on the decomposition $D_t\Gamma D_t$. To do so, a $\log$-type dynamic on $H_t$ could be considered and the estimation could be managed through the development of a suitable state-\textcolor{black}{space-based} setting.

\spacingset{0.7}

References

\spacingset{1}