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.
73,715 characters · 14 sections · 86 citation commands
Factor Multivariate Stochastic Volatility Models of High Dimension
\affil[a]{\it Faculty of Science and Technology, Keio University and Riken-AIP, Japan. [email removed]} \affil[b]{\it Faculty of Economics, Soka University, Japan. [email removed]}
\baselineskip = 7mm
{\bf{Keywords:}} Factor Model; Forecasting; Multivariate Stochastic Volatility.
The generalized autoregressive conditional heteroskedasticity (GARCH) and the stochastic volatility (SV) models are two popular families to specify and analyze the uncertainty of economic and financial time series. Multivariate GARCH (MGARCH) and Multivariate SV (MSV) models can capture time-varying covariance structures and are thus useful for forecasting variance-covariance matrices. Within MGARCH models, the dynamic conditional correlation (DCC) models of engle2002 and tse2002, the BEKK model of baba1985 and engle1995, and their variants are commonly used: see the surveys of bauwens2006 and boudt2019. The MSV process of harvey1994 serves as the foundational framework and has been extended in various ways, including the factor model of chib2006 and the dynamic correlation model of asai2009. For a broader discussion on MSV models and Bayesian Markov chain Monte Carlo (MCMC) techniques for their estimation, see chib2009 and kastner2017.
A major drawback of MGARCH and MSV models is the rapid increase in the number of parameters as the dimension of the observed vector grows. This not only complicates estimation in high dimensions but also degrades predictive performance due to overfitting. While enforcing parsimony in parameterization can mitigate this issue, overly restrictive specifications may fail to capture essential data features. Consequently, there exists a trade-off between model flexibility and parsimony.
A solution to the curse of dimensionality in MGARCH and MSV dynamics is to incorporate a factor-based structure, following the approaches of engle1990 and jacquier1999, respectively. These factor models for time-varying covariance processes align with asset pricing theory in finance. \textcolor{black}{While alexander2001 and weide2002 develop MGARCH models based on principal component analysis, barigozzi2017 consider general dynamic factor models with conditional volatilities.}
In the MSV literature, various factor MSV (fMSV) models have been proposed to promote parsimonious parameterizations; see, e.g., aguilar2000, chib2006, jacquier1999, liesenfeld2003, lopes2007, and pitt1999. Notably, most of these studies rely on Bayesian MCMC methods, with the exception of liesenfeld2003, which applies a maximum likelihood approach based on efficient importance sampling. However, their fMSV model includes only a single factor.
poignard2023_jtsa developed \textcolor{black}{an} OLS-based technique for estimating high-dimensional MSV models. Unlike existing MSV models, this approach allows for interactions between log-volatilities of different returns across different time periods, making it more flexible than the basic MSV model of harvey1994. Moreover, since the method is non-Bayesian, it avoids \textcolor{black}{the} pitfalls of MCMC-based techniques, such as prior specification challenges \textcolor{black}{and} numerical instability. Our main idea is to extend this OLS framework to the estimation of fMSV models. Specifically, we consider the following problem: given $T$ observations of a $p$-dimensional random vector $y_t$, we specify a factor decomposition \textcolor{black}{where} the common factors follow an MSV model, as in jacquier1999. The estimation procedure consists of two stages. First, we estimate the factor model \textcolor{black}{by maximum likelihood under suitable restrictions for identification following the framework of bai2012}. Then, using the estimated factor scores, we estimate the MSV model for the common factors. \textcolor{black}{Although mucher2025 recently developed a similar two-stage estimation approach based on bai2012, our method differs in the second stage. While mucher2025 employ the efficient method of moments, which is a computationally intensive technique, our approach offers computational efficiency.}
Our main contributions can be summarized as follows: we propose a factor model-based MSV estimation procedure \textcolor{black}{that does not rely on simulated-based techniques}; we establish consistency results of the estimators; the relevance of the method is illustrated by simulated experiments and real data applications.
The remainder of the paper is organized as follows. In Section (ref), we describe the fMSV model and the procedure to obtain the forecasts of covariance matrices. Section (ref) details the two-stage estimation method. The simulation experiments and the real data-based out-of-sample applications are provided in Section (ref) and Section (ref), respectively. All the asymptotic results and secondary technical details \textcolor{black}{are deferred to the Appendices}. The implementation of the fMSV model is available in the Github repository \url{https://github.com/Benjamin-Poignard/fMSV}.
Notations. Throughout this paper, we denote the cardinality of a set $E$ by $|E|$. For $\mbox{\boldmath $v$} \in {\mathbb R}^{d}$, the $\ell_p$ norm is $\|\mbox{\boldmath $v$}\|_p = \big(\sum^{\text{d}}_{k=1} |\mbox{\boldmath $v$}_k|^p \big)^{1/p}$ for $p > 0$, and $\|\mbox{\boldmath $v$}\|_{\infty} = \max_i|\mbox{\boldmath $v$}_i|$. We write $A^\top$ (resp. $\mbox{\boldmath $v$}^\top$) to denote the transpose of the matrix $A$ (resp. the vector $\mbox{\boldmath $v$}$). For a symmetric matrix $A$, $\lambda_{\min}(A)$ (resp. $\lambda_{\max}(A)$) is the minimum (resp. maximum) eigenvalue of $A$, and $\text{tr}(A)$ is the trace operator. For a matrix $B$, $\|B\|_s = \lambda^{1/2}_{\max}(B^\top B)$ and \textcolor{black}{$\|B\|_F=\{\text{tr}(B^\top B)\}^{1/2}$} are the spectral and Frobenius norms, respectively. 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 $d(d+1)/2$ vector that stacks the columns of the lower triangular part of $A \in {\mathbb R}^{d\times d}$. The matrix $I_d$ denotes the $d$-dimensional identity matrix. We denote by $O \in {\mathbb R}^{k \times l}$ the $k \times l$ zero matrix. For two matrices $A,B$, $A \otimes B$ is the Kronecker product. For two matrices $A,B$ of the same dimension, $A \odot B$ is the Hadamard product. For $f: {\mathbb R}^{d} \rightarrow {\mathbb R}$, we denote by $\nabla f$ the gradient or subgradient of $f$ and by $\nabla^2 f$ the Hessian of $f$.
We consider a $p$-dimensional vectorial stochastic process $(y_t)_{t=1,\cdots,T}$ following the fMSV structure:
where $f_t \in {\mathbb R}^m$ is the factor variable, $\Lambda \in {\mathbb R}^{p \times m}$ is the loading matrix, $\varepsilon_t = (\varepsilon_{1,t},\ldots,\varepsilon_{p,t})^\top \in {\mathbb R}^p$ is a random vector, which is independently and identically distributed (i.i.d.) with \textcolor{black}{a diagonal} covariance matrix $\Sigma_{\varepsilon}$, $h_t = (h_{1,t},\ldots,h_{m,t})^\top \in {\mathbb R}^m$ is a vector of log-volatilities, $D_t \in {\mathbb R}^{m \times m}$ is a diagonal matrix of volatilities, $\zeta_t = (\zeta_{1,t},\ldots,\zeta_{m,t})^\top \in {\mathbb R}^m$ is an i.i.d. random vector with covariance matrix $I_m$, $\mu =(\mu_1,\ldots,\mu_m)^\top \in {\mathbb R}^m$, $\Phi \in {\mathbb R}^{m \times m}$ \textcolor{black}{is non-diagonal}, $\Sigma_{\eta} \in {\mathbb R}^{m \times m}$ is the covariance matrix of $\eta_t$, and $(\varepsilon_t)_{t=1,\cdots,T}$, $(\zeta_t)_{t=1,\cdots,T}$, and $(\eta_t)_{t=1,\cdots,T}$ are mutually independent. \textcolor{black}{The empirical findings of andersen2001 supports the Gaussian assumption on $\eta_t$ in the context of realized volatility. Therefore, we assume $\eta_t \sim {\mathcal N}_{{\mathbb R}^m}(0, \Sigma_{\eta})$.}
We assume that all the eigenvalues of $\Phi$ are strictly smaller than one in modulus, which guarantees that $h_t$ and $y_t$ are strictly stationary processes. \textcolor{black}{We may relax this condition: see the discussion at the end of Subsection (ref).} Note that a non-diagonal $\Phi$ allows for cross-effects from other past variables, while the existing literature usually specifies a diagonal $\Phi$, except for poignard2023_jtsa. \textcolor{black}{We do not consider higher-order specifications for $h_t$, as empirical evidence favors the parsimonious MGARCH(1,1) model over higher-order alternatives, despite their theoretical interest.}
Suitable restrictions for the identification of the factor model can be found in, e.g., Table 1 in bai2012. In addition to restricting the factor model as $\Lambda_{kl}=0$ for $l>k$ and $\Lambda_{kk}=1$ for $k \leq m$, pitt1999 and chib2006 heavily rely on Bayesian MCMC \textcolor{black}{for estimation purposes}. This is a key difference with our proposed procedure: our first stage devoted to the factor model estimation builds upon the work of \textcolor{black}{bai2012, which lies within the likelihood framework. To ensure the identification of the factor model, we rely on condition IC2 of bai2012: $$\text{IC2:}\;\;M_f:={\mathbb E}[f_tf^\top_t] \; \text{is diagonal with distinct coefficients}, \;\; p^{-1}\Lambda^\top \Sigma^{-1}_{\varepsilon} \Lambda=I_m.$$ By equation ((ref)), ${\mathbb E}[f_tf^\top_t] = {\mathbb E}[{\mathbb E}[D^{1/2}_t\zeta_t\zeta^\top_t D^{1/2}_t|h_t]] = {\mathbb E}[D_t] \neq I_m$, which renders identification conditions IC3 and IC5 of bai2012 -- both imposing $M_f=I_m$ -- inappropriate for our setting. While their condition IC1 leaves $M_f$ unrestricted and their condition IC4 requires $M_f$ to be diagonal, both conditions impose a structure on $\Lambda$, similar to pitt1999 and chib2006. Furthermore, under condition IC2, ${\mathbb V}(f_{k,t})\neq {\mathbb V}(f_{l,t})$ for $k \neq l$. Moreover, ${\mathbb E}[f_t\varepsilon^\top_t]={\mathbb E}[D^{1/2}_t\zeta_t\varepsilon^\top_t]=O\in {\mathbb R}^{m \times p}$. By equation ((ref)), the unconditional variance-covariance ${\mathbb V}(y_t)$ of $y_t$ is ${\mathbb V}(y_t)=\Lambda M_f\Lambda^\top + \Sigma_{\varepsilon}$. Note that bai2012 assume that the covariance $\Sigma_{\varepsilon}$ of $\varepsilon_t$ is diagonal. The time-varying covariance matrix of $y_t$ conditional on $h_t$ in model ((ref))-((ref)) is defined by
}
Our procedure to obtain the forecasts of $H_t$, which will be referred as fMSV hereafter, can be broken down as \textcolor{black}{follows}:
Stage 1(a). Obtain $(\widehat{\Lambda},\widehat{\Sigma}_\varepsilon)$ via \textcolor{black}{bai2012 under IC2}.
\textcolor{black}{Stage 1(b)}. Given, $(\widehat{\Lambda},\widehat{\Sigma}_\varepsilon)$, \textcolor{black}{compute} the generalized least squares (GLS) estimator of $f_t$ as
See Section 14.7 of anderson2003introduction for its derivation.
Stage 2(a). Obtain the forecasts $\widehat{D}_{T+l}$ of $D_{T+l}$ using $\widehat{f}_t$ $(t=1,\ldots,T)$, modifying the work of poignard2023_jtsa.
Stage 2(b). Obtain the forecasts $\widehat{H}_{T+l}$ of $H_{T+l}$ as $\widehat{H}_{T+l} = \widehat{\Lambda} \widehat{D}_{T+l} \widehat{\Lambda}^\top + \widehat{\Sigma}_{\varepsilon}$ for $l=1,2,\ldots,H$.
In other words, the estimation of the fMSV model is decomposed into two stages. The first stage concerns the estimation of the factor model parameters, \textcolor{black}{denoted by $\theta_{\text{FM}}$, as described in Subsection (ref). The second stage relates to the estimation of the common-factor MSV model parameters, denoted by $\theta_{\text{MSV}}$, as detailed in Subsection (ref). Throughout this work, we assume that the number of factors $m$ is fixed and known. The value $m$ can be selected using, e.g., the method of bai2002 or onatski2010. } Note that poignard2023_jtsa consider the minimum mean square linear estimator (MMSLE) to obtain the forecasts of the log-volatilities.
One point in the procedure of poignard2023_jtsa requires a modification. The latter work relies on the logarithm of the squared return, $\log y_{i,t}^2$ $(i=1,\ldots,p)$, \textcolor{black}{which excludes the case $y_{i,t}=0$}. In the context of the fMSV model, \textcolor{black}{as it allows $f_{i,t}=0$, we use $\log (f_{i,t}^2 + c_i) - c_i/(f_{i,t}^2 + c_i)$, where $c_i = 10^{-4} s_i^2$ and $s_i^2$ is the sample mean of $f_{i,t}^2$. As discussed in Section 9.3 of fuller1996, we may use $c_i = 0.02 s_i^2$ if we consider more robust estimation against non-normality. However, our estimation technique does not require normal approximation for $\log f_{i,t}^2$. Instead, we need a slight change for $\log f_{i,t}^2$ to allow the case $f_{i,t}=0$. Our simulation and empirical results are robust to the setting $c_i = 10^{-5} s_i^2$. To perform Stage 2, $f_t$ will be replaced with $\widehat{f}_t$.} The next section details the two-stage estimation method.
\textcolor{black}{We denote by $\theta=(\theta^{\top}_{\text{FM}},\theta^{\top}_{\text{MSV}})^\top$ the vector of the fMSV model parameters, where $\theta_{\text{FM}}$ and $\theta_{\text{MSV}}$ correspond to the factor model parameters and stochastic volatility parameters, respectively. The estimation of $\theta_{\text{FM}}$ is detailed in Subsection (ref). Subsection (ref) describes the estimation of $\theta_{\textcolor{black}{\text{MSV}}}$.}
\textcolor{black}{The factor model parameter vector of equation ((ref)) is defined by $\theta_{\text{FM}}= (\theta^\top_{\Lambda},\theta^\top_{\Sigma_{\varepsilon}})^\top$, with $\theta_{\Lambda} = \text{vec}(\Lambda) \in {\mathbb R}^{pm},\; \theta_{\Sigma_{\varepsilon}}=\text{vech}(\Sigma_{\varepsilon}) \in {\mathbb R}^{p(p+1)/2}$. It belongs to the parameter space $\Theta_{\text{FM}} := \Theta_{\Lambda}\times\Theta_{\Sigma_{\varepsilon}} \subseteq {\mathbb R}^{p m}\times {\mathbb R}^{p(p+1)/2 }$.} Estimation methods for factor models can be divided into two: likelihood-based ones (see, e.g., bai2012,bai2016,bailiao2016,poignard2020,poignard_terada_2025) and PCA-based ones (see, e.g., stock2002a,fan2011,fan2013). Likelihood-based methods estimate $\Lambda,\Sigma_{\varepsilon}$ \textcolor{black}{and $M_f$}, and in a second stage, estimate the factors by, e.g., \textcolor{black}{the generalized least squares method}. PCA procedures usually estimate both $\Lambda$ and $f_t,t=1,\ldots,T$, and then obtain the estimator of $\Sigma_{\varepsilon}$.
\textcolor{black}{Throughout this paper, we will employ the maximum likelihood method under condition IC2 of bai2012 and under the diagonal assumption on $\Sigma_{\varepsilon}$, so that $\theta_{\Sigma_{\varepsilon}} \in {\mathbb R}^p$. Under IC2, $M_f:={\mathbb E}[f_t f^\top_t] \in {\mathbb R}^{m \times m}$ is diagonal with distinct coefficients, and $p^{-1}\Lambda^\top \Sigma^{-1}_{\varepsilon}\Lambda=I_m$. The factor model parameters $\Lambda,M_f,\Sigma_{\varepsilon}$ are estimated by Gaussian MLE, and the corresponding asymptotic properties are provided by bai2012. Now, following Section 8 of bai2012, the maximization of the likelihood function is performed by the EM algorithm and under condition IC3 given by: $$\text{IC3:}\;\; M_f=I_m,\;\; \text{and}\;\; p^{-1}\Lambda^\top\Sigma^{-1}_{\varepsilon}\Lambda\; \text{diagonal with distinct elements (arranged in decreasing order)}.$$ This implies that the EM algorithm will provide a solution for the variance ${\mathbb V}(y_t)= \Lambda\Lambda^\top+\Sigma_{\varepsilon}$, where the corresponding solution $(\widehat{\Lambda}^{\text{IC3}},\widehat{\Sigma}^{\text{IC3}}_{\varepsilon})$ satisfies IC3. More precisely, $\widehat{\Lambda}^{\text{IC3}} = \widetilde{\Lambda}V$, $\widehat{\Sigma}^{\text{IC3}}_{\varepsilon}=\widetilde{\Sigma}_{\varepsilon}$, with $\widetilde{\Lambda}, \widetilde{\Sigma}_{\varepsilon}$ being the values obtained at the final iteration of the EM algorithm, and $V$ is the orthogonal matrix containing the eigenvectors of $p^{-1}\widetilde{\Lambda}^\top\widetilde{\Sigma}^{-1}_{\varepsilon}\widetilde{\Lambda}$ corresponding to decreasing eigenvalues. To obtain the estimators of $\Lambda,\Sigma_{\varepsilon}$ under IC2, we compute $\widehat{\Lambda}^{\text{IC2}}=\widehat{\Lambda}^{\text{IC3}}\big(p^{-1}\widehat{\Lambda}^{\text{IC3}\top}\big(\widehat{\Sigma}^{\text{IC3}}_{\varepsilon}\big)^{-1}\widehat{\Lambda}^{\text{IC3}}\big)^{-1/2}$, and $\widehat{M}^{\text{IC2}}_{f}=p^{-1}\widehat{\Lambda}^{\text{IC3}\top}\big(\widehat{\Sigma}^{\text{IC3}}_{\varepsilon}\big)^{-1}\widehat{\Lambda}^{\text{IC3}}$. Under the diagonal assumption on $\Sigma_{\varepsilon}$, we have $\widehat{\Sigma}^{\text{IC2}}_{\varepsilon}=\widehat{\Sigma}^{\text{IC3}}_{\varepsilon}=\widetilde{\Sigma}_{\varepsilon}$. Using the estimator $\widehat{\theta}_{\text{FM}}=(\widehat{\theta}^\top_\Lambda,\widehat{\theta}^\top_{\Sigma_{\varepsilon}})^\top$ of $\theta_{\text{FM}}$ under IC2, we compute the GLS estimator $\widehat{f}_t = (\widehat{\Lambda}^\top \widehat{\Sigma}_\varepsilon^{-1} \widehat{\Lambda})^{-1} \widehat{\Lambda}^\top \widehat{\Sigma}_\varepsilon^{-1} y_t$ of $f_t$, for all $t$. By Theorem 6.1 of bai2012, under large $p$, $\widehat{f}_t$ is asymptotically normally distributed.}
\textcolor{black}{In Appendix (ref), we show the uniform consistency of $\widehat{f}_t$ under suitable moment conditions and scale conditions on $(p,T)$. This property will facilitate the proofs in the asymptotic analysis of the MSV parameter $\theta_{\text{MSV}}$ detailed in Appendix (ref): it will allow us to control the estimation error resulting from the estimation of the factor variable} when plugged in the loss functions employed for the estimation of $\theta_{\text{MSV}}$.
\textcolor{black}{As for the choice of the number of factors $m$, one may rely on} onatski2010's method based on an eigenvalue threshold estimator using the sample variance-covariance matrix of $y_t$. To efficiently reduce dimensionality, it is generally preferable to have a small number of factors, typically \textcolor{black}{$m\leq 5$}. As noted in deNard2021, there is no consensual selection procedure of the number of factors, making it important to assess the sensitivity of the factor-model-based MSV performance with respect to $m$. Therefore, in our experiments, we set $m \in \{1,2,3\textcolor{black}{, 4, 5}\}$.
\textcolor{black}{We now consider the estimation of $\theta_{\text{MSV}}$, which can be partitioned into two subvectors $\theta_1$ and $\theta_2$ corresponding to the two-stage procedure detailed hereafter in Step 1 and Step 2. }
By transforming $f_t$ to $f_t^\ell =(\log f_{1,t}^2,\ldots,\log f_{m,t}^2)^\top$, we obtain the state space form:
where $\nu=(\nu_1,\ldots,\nu_m)^\top$, $\xi_t = (\xi_{1,t},\ldots,\xi_{m,t})^\top$, and $\alpha_t = h_t - \mu$ with $\nu_i = \mu_i + {\mathbb E}[\log \zeta_{i,t}^2]$ and $\xi_{i,t} = \log \zeta_{i,t}^2 - {\mathbb E}[\log \zeta_{i,t}^2]$. We may apply the Kalman filter to estimate the \textcolor{black}{parameters} in ((ref)), as suggested by Harvey et al. (1994). \textcolor{black}{However, since we assume that $\Phi$ is non-diagonal, estimation by the Kalman filter can be a challenging task, even in a low-dimensional setting. Therefore}, in the same spirit as in poignard2023_jtsa, we estimate $\theta_{\text{MSV}}$ based on the following steps:
Step 1. Since $f_t^\ell$ is the sum of a VAR(1) process, $\alpha_t$, and an i.i.d. error $\zeta_t$, case (i) in Section 3 of granger1976 suggests that $f_t^\ell$ has a VARMA(1,1) representation, which leads a VAR($\infty$) model that will be employed to approximate the unobserved error. More precisely, \textcolor{black}{consider $x_t = (\log (f_{1,t}^2 + c_1) - c_1/(f_{1,t}^2 + c_1),\ldots, \log (f_{m,t}^2 + c_m) - c_m/(f_{m,t}^2 + c_m))^\top$, as detailed in the previous section, in order to obtain the state space form,
where $\nu^*=(\nu_1^*,\ldots,\nu_m^*)^\top$, $\xi_t^* = (\xi_{1,t}^*,\ldots,\xi_{m,t}^*)^\top$,
Since we set $c_i = 10^{-4}s_i^2$, where $s_i^2$ is the sample mean of $f_{it}^2$, $c_i e^{-h_{i,t}}$ is a negligible quantity compared to $\zeta_{i,t}^2$. Hence we treat $\xi_t^*$ as a white noise process with mean zero and covariance matrix $\Sigma_{\xi} = {\mathbb E}[\xi_t \xi_t^\top]$, that is $\xi_t^* \sim WN(0,\Sigma_{\xi})$. It holds when $c_i=0$. As detailed in Appendix (ref), we have the VARMA(1,1) representation for $x_t$ given by
where $\Upsilon$ and $\Sigma_u$ are defined by Appendix (ref).} \textcolor{black}{From the representation, we deduce the VAR($\infty$) form of $x_t$ that we approximate by the VAR($q$) model} $x_t = \psi^* +\sum_{k=1}^q \Psi_k x_{t-k} + u^{(q)}_t$, with $u^{(q)}_t = u_t +\sum_{k>q}\Psi_k x_{t-k}$. The term $\sum_{k>q}\Psi_kx_{t-k}$ \textcolor{black}{corresponds to the remainder term resulting} from the truncation of the VAR($\infty$). Under suitable parameter decay conditions, it can be shown that ${\mathbb E}[\|u^{(q)}_t-u_t\|^r_2]$ is negligible when $q$ is large enough, assuming the existence of the $r$-th moment of $u_t$, as in chang2002 for univariate processes and chang2006 for multivariate processes. In our framework, the number of factors $m$ is fixed, \textcolor{black}{$q$ may potentially diverge with the sample.} We apply a \textcolor{black}{penalty function} to suitably discard the irrelevant past variables. It aims to estimate the sparse support and identify an optimal $q$ \textcolor{black}{lag} from which the past variables can be considered as negligible. \textcolor{black}{Define $\theta_1 := \text{vec}(\underline{\Psi}) \in \Theta_{1} \subseteq {\mathbb R}^{d_1}$, $d_1:=m+qm^2$, with $\underline{\Psi} =(\psi^*,\Psi_1,\cdots,\Psi_q)$ the parameters of the VAR($q$) model and $\Theta_1$ the parameter space of $\theta_1$. Define} the non-penalized loss ${\mathbb L}_T: {\mathbb R}^{m(q+1)T} \times \Theta_{1} \rightarrow {\mathbb R}$. \textcolor{black}{Let $\mathbf{F} = (F^\top_1,\ldots,F^\top_T)^\top \in {\mathbb R}^{m(q+1)T}$, $F_t = (f^\top_t,\ldots,f^\top_{t-q})^\top \in {\mathbb R}^{m(q+1)}$ and let $f_0,\ldots,f_{1-q}$ be initial values (chosen equal to zero). \textcolor{black}{Define $Z_{q,t-1}=(1,x^\top_{t-1},\cdots,x^\top_{t-q})^\top\in {\mathbb R}^{1+mq}$, $x_{j,t} = \log (f_{j,t}^2 + c_j) - c_j/(f_{j,t}^2 + c_j)$ for $j=1,\ldots,m$.} When replacing $f_t$ by its estimator $\widehat{f}_t$, \textcolor{black}{the variables} $\mathbf{F}, F_t, \textcolor{black}{Z_{q,t-1},x_t}$ are denoted by $\widehat{\mathbf{F}}$, $\widehat{F}_t,\textcolor{black}{\widehat{Z}_{q,t-1},\widehat{x}_t}$.} The loss ${\mathbb L}_T(\widehat{\mathbf{F}};\theta)$ is associated to a continuous function $\ell : {\mathbb R}^{m(q+1)} \times \Theta_{1} \rightarrow {\mathbb R}$ that can be written as
The true parameter is denoted by $\theta_{01}$ and is assumed sparse with true support $s := \text{card}({\mathcal S})$, with ${\mathcal S}:=\{k=1,\cdots,d_1: \theta_{01,k}\neq 0\}$. To estimate the latter support, we rely on the penalized M-estimation criterion given by
The penalty function \textcolor{black}{is} the convex adaptive LASSO penalty of zou2006adaptive, where $\lambda_T$ is the tuning parameter and $\tau(\widetilde{\theta}_{1,k})$ is a random weight defined by $\tau(\widetilde{\theta}_{1,k})=|\widetilde{\theta}_{1,k}|^{-\gamma}$, with $\widetilde{\theta}_{1,k}$ a consistent first step estimator. It is defined as the solution:
Theorems (ref) and (ref) in Appendix (ref) establish the consistency of $\widetilde{\theta}_1$ and $\widehat{\theta}_1$, respectively. Theorem (ref) provides the \textcolor{black}{“sparsistency” property of $\widehat{\theta}_1$, that is the recovery of the true zero entries with probability tending to one}: see lam2009 for further details on sparsistency.
In our experiments, the selection of an optimal $\lambda_T$ is performed by \textcolor{black}{a cross-validation method that accounts for the dependent nature of the data. We will set $q=10$ when solving ((ref)) and ((ref)). In the latter problem, we set $\gamma=1$. All the implementation details are provided in Appendix (ref)}.
\noindentStep 2. Define the approximated error $\widehat{u}^{(q)}_t = x_t - \widehat{\psi}^*-\sum^q_{k=1}\widehat{\Psi}_k x_{t-k}$ obtained in Step 1. Given $\widehat{\theta}_1$, we consider the regression \[ x_t = c^* + \Phi x_{t-1}+ \Xi \widehat{u}^{(q)}_{t-1} + v_t, \] where the vector of parameters is $\theta_2:=(c^{*\top},\theta^\top_{\Phi},\theta^\top_{\Xi})^\top \in \Theta_2 \subseteq {\mathbb R}^{d_2}$, $d_2 := m(1+2m)$, and since we replace $u_{t-1}$ by $\widehat{u}^{(q)}_{t-1}$, $(v_t)$ is the error term for this auxiliary regression.
\textcolor{black}{The estimation of $\theta_2$} relies on the second step loss ${\mathbb G}_T : {\mathbb R}^{m(q+1)T} \times \Theta_1 \times \Theta_2 \rightarrow {\mathbb R}$, where ${\mathbb G}_T(\widehat{\mathbf{F}};\widehat{\theta}_1;\theta_2)$ is the loss associated to a continuous function $g: {\mathbb R}^{m(q+1)} \times \Theta_1 \times \Theta_2 \rightarrow {\mathbb R}$. \textcolor{black}{As in Step 1, we replace the latent $x_t$ with $\widehat{x}_t$, so that the second step loss is}
where $\widehat{K}_{t-1}(\widehat{\theta}_1) = (1,\textcolor{black}{\widehat{x}}^\top_{\textcolor{black}{t-1}},\textcolor{black}{\widehat{\overline{u}}}^{(q)\top}_{t-1})^\top \in {\mathbb R}^{1+2m}$ and $\Gamma = (c^*,\Phi,\Xi) \in {\mathbb R}^{m\times (1+2m)}$. The dependence on the first step estimator is through \textcolor{black}{$\widehat{\overline{u}}^{(q)}_t=\widehat{x}_t - \widehat{\psi}^*-\sum^q_{k=1}\widehat{\Psi}_k \widehat{x}_{t-k}$}. The problem of interest is
Theorem (ref) in Appendix (ref) establishes the consistency of $\widehat{\theta}_2$.
\noindentStep 3. From the decomposition of the unconditional variance-covariance of $x_t$ given by
where $\Sigma_x = {\mathbb E}[(x_t-\gamma)(x_t-\gamma)^\top]$, $\Sigma_{\xi} = {\mathbb E}[\xi_t \xi_t^\top]$, and $\Sigma_{\alpha} = {\mathbb E}[\alpha_t \alpha_t^\top]$, the estimators of $\Sigma_{\zeta},\Sigma_{\alpha}$ are deduced as
where $0 < r < 1$ and $S_x$ is the sample covariance matrix of \textcolor{black}{$\widehat{x}_t$}. This ad hoc method aims to treat the positive-definiteness of the estimators and to deal with the high-dimensional issue of $\widehat{\Sigma}_{\textcolor{black}{\xi}}$. While we consider a naive decomposition based on equation ((ref)) for the former, we set $r = (\pi^2/2)(m^{-1} \mathrm{tr} (S_{x}))^{-1}$ in ((ref)). Here, $ \pi^2/2$ is the value of ${\mathbb E}[\textcolor{black}{\xi}_{i,t}^2]$ when \textcolor{black}{$\zeta_{i,t}$} follows the standard normal distribution. The ad hoc estimators yield $\mathrm{tr} (\widehat{\Sigma}_{\textcolor{black}{\xi}}) = r \mathrm{tr}(S_x) = m \pi^2/2$ and $\mathrm{tr} (\widehat{\Sigma}_{\alpha}) = \mathrm{tr}(S_x) - m \pi^2/2$. Using this approach, we are able to estimate $\mathrm{tr}({\Sigma}_{\textcolor{black}{\xi}})$ and $\mathrm{tr}({\Sigma}_{\alpha})$ with accuracy and consistency, respectively. More importantly, the computational cost is negligible, compared with alternative estimators (e.g., the GMM type method) that would require a numerical optimization with constraints on the positive-definiteness of $S_{\xi}$ and $S_x-S_{\xi}$, where $S_{\xi}$ is an alternative estimator.
\textcolor{black}{We dot not impose the restriction $\rho(\Phi)<1$ apriori. For economic and financial data, we need to control the case where several eigenvalues of $\widehat{\Phi}$ exceed 1, using the random walk model. By replacing such eigenvalues of $\widehat{\Phi}$ by one’s, and setting corresponding values of $\widehat{c}$ to be zero, we may construct a state space model for $x_t$ to obtain out-of-sample and in-sample forecasts by the Kalman filter and smoother, respectively.}
\textcolor{black}{ In Appendix (ref), we compare the computational efficiency of our two-stage procedure with the MCMC approach proposed by kastner2017. The empirical results show that our method is computationally more efficient than the MCMC benchmark. }
This section presents numerical experiments on the in-sample comparison between a true $p \times p$ dynamic variance-covariance $(H_t)$ and the conditional forecast $(\widehat{H}_t)$ obtained from a specific model.
The simulated $p$-dimensional observations $(y_t)$ satisfy:
where $t(\gamma)$ is a centered Student distribution with $\delta$ degrees of freedom, $\delta >2$. Hereafter, $\delta=3$. The true variance-covariance $(H_t)$ is based on two data generating processes (DGPs):
We set the sample size $T=2000$ and dimension \textcolor{black}{$p=20, 100, 500$} and we consider $100$ independent batches for $(y_t)$. Once a series is simulated, we estimate the model using the full sample and compare this estimator with the true variance-covariance to measure its accuracy.
\textcolor{black}{We measure the in-sample statistical accuracy of the variance-covariance models based on the following matrix distances proposed by laurent2012:
where $H_t$ is the true variance-covariance and $\widehat{H}_t$ is the estimated variance-covariance matrix at time $t$. $\text{D}_E$ and $\text{D}_F$ are the Euclidean distance, the squared Frobenius norm, respectively, and are closely related. $\text{D}_E$ evaluates how close the individual covariance entries are, whereas $\text{D}_F$ is a matrix-based MSE distance. $\text{D}_S$ is the Stein loss: it is the scale-invariant Gaussian likelihood, which is asymmetric with respect to over and/or under predictions, where under predictions are significantly penalized. In the same spirit, the measure $\text{D}_b$ proposed by laurent2012 is also asymmetric with respect to over and/or under predictions, but significantly penalizes over predictions. The coefficient of asymmetry $b$ is set as $b=3$.}
Equipped with these DGPs, we generate the observations $y_t,t=1,\ldots,T$ and estimate the variance-covariance using the models:
\textcolor{black}{We set $m=1,2,3, \textcolor{black}{4,5},$ for $\text{fGARCH}_m$, $\text{fMSV}_m$ to assess the sensitivity of the method with respect to the factor dimension. }
\textcolor{black}{Tables (ref)–(ref) report the in-sample accuracy measures $\text{D}_E,\text{D}_F,\text{D}_S,\text{D}_3$, averaged over 100 independent simulation batches, for DGP 1 and DGP 2, respectively. These matrix distances emphasize different aspects of covariance forecast accuracy and therefore highlight distinct strengths and weaknesses of the competing models. The distances $\text{D}_E$ and $\text{D}_F$ measure entry-wise covariance accuracy. Under DGP 1, the sBEKK model performs best according to these criteria, indicating an excellent fit of individual covariance entries. The proposed fMSV model becomes increasingly competitive as the dimension grows. However, the sBEKK model exhibits substantial distortions in the inverse covariance structure, as emphasized by the Stein loss $\text{D}_S$, for which the DCC model (resp. factor-based models) achieves the best performance when $p \leq 100$ (resp. $p=500$). In terms of the asymmetric $\text{D}_3$, both sBEKK and fMSV avoid systematic volatility overestimation, whereas the fGARCH specifications tend to inflate dominant eigenvalues. Under DGP 2 (Table (ref)), the sBEKK model becomes severely misspecified across all metrics. In contrast, factor-based models such as fMSV and fGARCH provide a substantially better fit with the true variance-covariance, particularly when $m \geq 2$. In terms of $\text{D}_S$, both fGARCH and fMSV models with $m \geq 2$ deliver the best performances. Moreover, the $m=2$ factor specifications yield the smallest values of $\text{D}_3$, indicating superior control of volatility overprediction. Interestingly, when $m=1$, both fGARCH and fMSV perform poorly across all criteria, underscoring the importance of selecting a sufficiently rich factor structure to capture the key features of the data. Overall, the fMSV model exhibits the most balanced performance across the different accuracy measures and DGPs.}
In this section, the relevance of our method is compared with existing competing models through out-of-sample forecasts based on real financial data. The performances will be assessed through measures of statistical accuracy and economic performances. The accuracy of the forecasts will be ranked by the model confidence set (MCS) procedure of hansen2011, which provides a testing framework for the null hypothesis of forecast equivalence across subsets of models.
We consider hereafter the stochastic process $(y_t)_{t\in{\mathbb Z}}$ in ${\mathbb R}^p$ of the log-stock returns, where $y_{j,t} = 100\times\log(P_{j,t}/P_{j,t-1}), 1 \leq j \leq p$ with $P_{j,t}$ the stock price of the $j$-th index at time $t$. We study three financial portfolios: a low-dimensional portfolio of daily log-returns composed of the MSCI stock index based on the sample December 1998 -- March 2018, which yields a sample size $T = 5006$, and for 23 countries\footnote{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, the United-States}; a mid-dimensional portfolio of daily log-returns listed in the S&P 100, where we selected $94$ firms\footnote{S&P 100 indices excluding AbbVie Inc., Dow Inc., General Motors, Kraft Heinz, Kinder Morgan and PayPal Holdings} that have been continuously included in the index over the period February 2010 -- January 2020, i.e., with a sample size $T=2500$; a high-dimensional portfolio of daily log-returns listed in the S&P 500, where we selected $480$ firms\footnote{S&P 500 indices excluding Etsy Inc., Solaredge Technologies Inc., PayPal, Hewlett Packard Enterprise, Under Armour (class C), Fortive, Lamb Weston, Ingersoll Rand Inc., Ceridian HCM, Linde PLC, Moderna, Fox Corporation (class A and B), Dow Inc., Corteva Inc., Amcor, Otis Worldwide, Carrier Global, Match Group, Viatris Inc.} that have been continuously included in the index over the period September 2014 -- January 2022, i.e., with a sample size $T=1901$. Hereafter, the portfolios are denoted by MSCI, S&P 100 and S&P 500, respectively\footnote{The MSCI, S&P 100 and S&P 500 data can be found on https://www.msci.com/, https://finance.yahoo.com and https://macrobond.com, respectively. MSCI and S&P 100 data are publicly available in \url{https://github.com/Benjamin-Poignard/fMSV}; the S&P 500 data requires a license and is thus not publicly available}.
We measure the statistical accuracy of the variance-covariance models based on the losses defined in ((ref)). Since $H_t$ is the underlying variance-covariance of $y_t$ and is unobserved, we replace it with a proxy for the realized covariance. This proxy is defined as $(1-a) y_ty^\top_t + a T^{*-1} \sum_{s=t-T^*+1}^t y_s y^\top_s$ ($t>T^*$), where $a=0.01$ and $T^*$ denotes the length of the in-sample period. This proxy represents a weighted average of $y_ty^\top_t$, which serves as a consistent (though noisy) estimator of $H_t$, and the rolling sample covariance matrix up to time $t$. The variance-covariance $\widehat{H}_t$ is estimated in-sample and the losses $\text{D}_E, \text{D}_F$, \textcolor{black}{$\text{D}_S$, $\text{D}_3$} are evaluated over the out-of-sample periods. The economic performances are assessed through the Global Minimum Variance Portfolio (GMVP) investment problem. The latter problem at time $t$, in the absence of short-sales constraints, is defined as
where $\widehat{H}_t$ is the $p \times p$ one-step ahead forecast of the conditional variance-covariance matrix built at time $t-1$ and $\iota$ is a $p \times 1$ vector of $1$'s. The explicit solution is given by $\omega_t = \widehat{H}^{-1}_t \iota/\iota^{\top}\widehat{H}^{-1}_t\iota$: as a function depending only on $\widehat{H}_t$, the GMVP performance essentially depends on the precise measurement of the variance-covariance matrix. The following out-of-sample performance metrics (annualized) will be reported: AVG as the average of the out-of-sample portfolio returns, multiplied by $252$; SD as the standard deviation of the out-of-sample portfolio returns, multiplied by $\sqrt{252}$; IR as the information ratio computed as $\textbf{AVG}/\textbf{SD}$. The key performance measure is the out-of-sample SD. The GMVP problem essentially aims to minimize the variance rather than to maximize the expected return (that is high IR). As an alternative to the GMVP, we examine the performance of the risk parity portfolio (RPP), described in Appendix (ref). The RPP aims to equalize risk contributions across all asset classes in order to maximize diversification benefits. We evaluate its performance, expecting a high IR.
The out-of-sample analysis evaluates statistical accuracy using $\text{D}_E$, $\text{D}_F$, $\text{D}_S$, and $\text{D}_3$ as well as GMVP and RPP performances for the competing models described in Subsection (ref). In addition, we consider \textcolor{black}{four} benchmark procedures: \textcolor{black}{the equally weighted portfolio $1/p$, the sample covariance matrix (SCov)}, the geometric-inverse shrinkage (GIS) estimator of ledoit2022, a nonlinear shrinkage estimator based on the symmetrized Kullback-Leibler loss, and the linear shrinkage towards one-parameter covariance estimator (Cov1Para) of ledoit2004. \textcolor{black}{The last three models are estimated using the in-sample period.} The MCS test for ranking the statistical accuracy takes $\text{D}_E,\text{D}_F$, $\text{D}_S$, and $\text{D}_3$ as loss functions. The \textcolor{black}{GMVP} performances are ranked out-of-sample by the MCS test which takes the distance based on the difference of the \textcolor{blue}{squared} returns of portfolios $i$ and $j$, defined as $u_{ij,t} = \big(\mu_{i,t}-\overline{\mu}_i\big)^2 - \big(\mu_{j,t}-\overline{\mu}_j\big)^2$, with $\mu_{k,t}=w^{\top}_{k,t} y_t$ the portfolio return at time $t$, where $w_{k,t}$ represents the GMVP weight deduced from variance-covariance model $k$, and $\overline{\mu}_k$ is the average portfolio return over the period. For the RPP performances, we set $u_{ij,t} = - \mu_{i,t}/s_i + \mu_{j,t}/s_j$ using the RPP weights with $s_k$ the standard deviation of the RPP return over the period. The MCS test is evaluated at the $10\%$ level based on the range statistic and with block bootstrap with $10,000$ replications: see hansen2003 for further technical procedures. We report below the out-of-sample periods used for the test accuracy, together with the number of factors $\widehat{m}$ selected in-sample by the variance-covariance eigenvalue-based procedure of onatski2010 on an indicative basis:
Table (ref) reports the average annualized values of $\text{D}_E$, $\text{D}_F$, \textcolor{black}{$\text{D}_S$, and $\text{D}_3$}, along with the MCS results. \textcolor{black}{The $\text{fMSV}_1$ model achieves the lowest distances for $\text{D}_E$, $\text{D}_F$, and $\text{D}_3$, whereas the results for $\text{D}_S$ vary, favoring the DCC for the low/mid dimensional portfolio and the $\text{fMSV}_5$ model for the high-dimensional portfolio. As detailed in Subsection (ref), $\text{D}_S$ and $\text{D}_3$ are asymmetric with respect to over- and under-predictions. $\text{D}_3$ yields rankings similar to those obtained from $\text{D}_E,\text{D}_F$, whereas the Stein loss $\text{D}_S$ leads to different conclusions. This indicates that relative and inverse-covariance errors drive the differences across models. For the MSCI and S&P 500 portfolios, the selection of fMSV$_1$ is in line with the number of factors determined by the method of onatski2010. Regarding the S&P 100 portfolio, while the approach of onatski2010 chose $m=3$, the distance measures tend to favor simpler models for prediction. This discrepancy implies that more parsimonious models are favored for forecasting as they exclude overfitting noise in the estimation period. In several cases, fMSV models with $m \geq 2$ are included in the MCS.}
Table (ref) displays the GMVP performance results, which can be summarized as follows (note that, unless stated otherwise, the results refer to SD). Overall, the fMSV models provide the best economic performance. Specifically, the $\text{fMSV}_5$ model achieves the lowest standard deviations for the MSCI and S&P 100 portfolios, whereas the S&P 500 portfolio favors the $\text{fMSV}_3$ model. \textcolor{black}{In the latter portfolio, among the benchmark procedures, the GIS is competitive with GARCH-class models but exhibits a higher SD than the fMSV models.} Furthermore, $\text{fMSV}_m$ $(m \geq 2)$ models are included in the MCS for the S&P 500 portfolio. Unlike the statistical accuracy results in Table (ref), these findings indicate that factor models with more factors yield better performances. In portfolio construction, the flexibility inherent in factor models with several factors may have contributed to optimizing asset allocation, particularly during periods of high market risk. Focusing on the $\text{D}_S$ measure and GMVP performance for the fMSV models, $\text{D}_S$ and SD select the same model ($\text{fMSV}_5$) for both MSCI and S&P100 in Tables (ref) and (ref). Regarding S&P500, while Table (ref) shows that $\text{fMSV}_3$ and $\text{fMSV}_5$ are included in the MCS under $\text{D}_S$, Table (ref) indicates that fMSV$_3$ achieves the smallest SD. Overall, these results suggest that $\text{D}_S$ and SD tend to select similar models.
The RPP results are displayed in Table (ref), and indicate that most of IRs are lower than those of the GMVP in Table (ref). In other words, the GMVP outperforms the RPP for the datasets. Reflecting the result, there are no major differences among the models.
Tables (ref) and (ref) illustrate that fMSV models generally outperform competing models, suggesting several key insights: (i) the stochastic nature of fMSV models captures market dynamics more effectively than the deterministic structure of GARCH class models; (ii) fMSV models with a larger number of factors provide the necessary flexibility to capture dynamic inter-asset relationships; and (iii) the strong performance of DCC models in certain cases puts a stress on the importance of incorporating dynamic correlations in covariance estimation.
The paper considers a two-stage estimation where the factor model is estimated in the first stage, while the second stage estimates the MSV processes of the estimated factors. We provide some asymptotic results for the second stage estimators. Monte Carlo results using the DGPs based on MGARCH models indicate that the forecasts of the fMSV model perform well compared with non-factor models. The empirical results based on real data show that the performances of the forecasts of the fMSV are better than those of the competing MGARCH models.
Several directions relating to the theoretical and the empirical analysis on fMSV models can be considered. The first topic concerns the variance-covariance matrix of the idiosyncratic errors. We may extend it to allow stochastic volatility, as in chib2006. The second one is to accommodate realized covariance and asymmetric effects, as in asai2015. \textcolor{black}{The third direction is to incorporate dynamic factors, as explored by barigozzi2017 and lam2012. Fourth, while our current framework is based on an additive factor structure, we could consider multiplicative volatility factors by extending the frameworks of ding2025 and ray2000.} We shall leave these issues for the future research.