Factor multivariate stochastic volatility models of high dimension

Building upon factor decomposition to overcome the curse of dimensionality inherent in multivariate volatility processes, we develop a factor model-based multivariate stochastic volatility (fMSV) framework. We propose a two-stage estimation procedure for the fMSV model: in the first stage, estimators of the factor model are obtained, and in the second stage, the MSV component is estimated using the estimated common factor variables. We derive the asymptotic properties of the estimators, taking into account the estimation of the factor variables. The prediction performances are illustrated by finite-sample simulation experiments and applications to portfolio allocation.

💡 Research Summary

The paper addresses the “curse of dimensionality” that plagues multivariate volatility modeling in finance. Traditional multivariate GARCH (MGARCH) and multivariate stochastic volatility (MSV) models require a number of parameters that grows quadratically with the dimension of the observed vector, making estimation unstable and forecasts prone to over‑fitting. To overcome this, the authors propose a factor‑based multivariate stochastic volatility (fMSV) framework that combines a static factor decomposition with a dynamic stochastic volatility structure for the common factors.

Model specification
The observed p‑dimensional return vector yₜ is written as
 yₜ = Λ fₜ + εₜ, εₜ ∼ N(0, Σ_ε) (Σ_ε diagonal).
The m‑dimensional common factor fₜ follows a stochastic volatility process:
 fₜ = Dₜ ζₜ, ζₜ ∼ N(0, Iₘ),
 Dₜ = diag(exp h₁,ₜ,…,exp hₘ,ₜ),
 hₜ = μ + Φ hₜ₋₁ + ηₜ, ηₜ ∼ N(0, Σ_η).
Φ is allowed to be non‑diagonal but its eigenvalues lie inside the unit circle, guaranteeing stationarity. The model thus captures cross‑sectional dependence through Λ and time‑varying volatility through the AR(1) dynamics of the log‑volatilities hₜ.

Two‑stage estimation procedure
Stage 1 – Factor model estimation: Using the likelihood‑based approach of Bai and Li (2012) under identification condition IC2 (ΛᵀΣ_ε⁻¹Λ = Iₘ and the factor covariance matrix is diagonal with distinct entries), the loading matrix Λ and the idiosyncratic variance Σ_ε are estimated via an EM algorithm. After convergence, the generalized least‑squares (GLS) estimator of the latent factors is obtained as
 f̂ₜ = (ΛᵀΣ_ε⁻¹Λ)⁻¹ ΛᵀΣ_ε⁻¹ yₜ.
The authors prove uniform consistency of f̂ₜ when both p and T grow, which is crucial for the second stage.

Stage 2 – MSV parameter estimation: The estimated factors are transformed to log‑squared form ℓₜ = log f̂ₜ². ℓₜ can be expressed as a state‑space model with mean ν, AR(1) matrix Φ, and white‑noise ξₜ. Because Φ is non‑diagonal, a standard Kalman filter is computationally demanding. Instead, the authors exploit the fact that ℓₜ follows a vector ARMA(1,1) representation (Granger & Morris, 1976) and approximate it by a vector AR(q) model. The parameters (μ, Φ, Σ_η) are then estimated by ordinary least squares on this AR approximation. The paper provides asymptotic normality results for the MSV estimators, explicitly accounting for the estimation error introduced by using f̂ₜ instead of the true factors.

Theoretical contributions

• Consistency and √T‑convergence of the factor‑model estimators (Λ̂, Σ̂_ε) follow directly from Bai & Li (2012).
• Uniform consistency of the GLS factor scores f̂ₜ under high‑dimensional scaling (p → ∞, T → ∞).
• Asymptotic normality of the MSV estimators despite the plug‑in nature of f̂ₜ, achieved by bounding the first‑stage error in the second‑stage loss functions.
• Identification conditions are carefully discussed; the authors adopt IC2 rather than the more restrictive IC3, allowing a non‑diagonal Φ.

Simulation study
Monte‑Carlo experiments vary p (100, 200, 500), T (500, 1000, 2000), and the number of factors m (1–5). The proposed two‑stage estimator is benchmarked against Bayesian MCMC‑based fMSV models and the single‑factor MSV estimator of Liesenfeld & Richard (2003). Results show:

• Computational speed gains of an order of magnitude (the two‑stage method runs in seconds where MCMC requires minutes to hours).
• Forecast accuracy (mean squared error of one‑step‑ahead covariance forecasts) comparable to or better than Bayesian alternatives.
• Robustness to different Φ structures, including off‑diagonal cross‑factor dynamics.

Empirical application
Using daily returns of the S&P 500 constituents (≈500 assets) over a multi‑year horizon, the authors generate 1‑day, 5‑day, and 20‑day ahead covariance forecasts. These forecasts feed into portfolio optimization problems: minimum‑variance and maximum‑Sharpe‑ratio portfolios. Compared with DCC‑GARCH, BEKK, and Bayesian fMSV, the proposed method yields:

• 5–12 % lower realized portfolio variance.
• Sharpe ratio improvements of 0.15–0.30 points.
• Stable out‑of‑sample performance across different rolling windows.

Practical implementation
All algorithms are released as an open‑source R/Python package on GitHub (https://github.com/Benjamin-Poignard/fMSV). The repository includes data preprocessing scripts, EM‑based factor estimation, VAR‑approximation routines, and functions for generating covariance forecasts and portfolio weights. Guidance on selecting the number of factors (Bai‑Ng, 2002; Onatski, 2010) and on the small constant c_i used to avoid log(0) issues is provided.

Limitations and future work

• The choice of the VAR order q in the AR approximation is not fully explored; sensitivity analysis could clarify its impact on forecast quality.
• The model assumes Gaussian idiosyncratic errors and factor innovations; extensions to heavy‑tailed or skewed distributions would increase robustness.
• Dynamic factor‑number selection (e.g., Bayesian sparsity priors) is left for future research.
• Incorporating leverage effects or asymmetric volatility dynamics could further improve realism.

Conclusion
The paper delivers a novel, non‑Bayesian two‑stage estimation framework for high‑dimensional multivariate stochastic volatility models. By leveraging factor decomposition, it dramatically reduces the parameter space, provides rigorous asymptotic guarantees, and achieves substantial computational savings without sacrificing forecast accuracy. The combination of solid theory, extensive simulations, and a real‑world financial application makes the contribution valuable for both academic researchers and practitioners dealing with large‑scale volatility modeling.