Motivated by empirical evidence from the joint behavior of realized volatility time series, we propose to model the joint dynamics of log-volatilities using a bivariate fractional Ornstein-Uhlenbeck (2fOU) process. This is a mean reverting Gaussian process with fractal features living in R2. It is the solution of the Langevin equation with the multivariate fractional Brownian motion, in the sense of Amblard et al. (2012), as the driving term. This model is a multivariate version of the Rough Fractional Stochastic Volatility (RFSV) model proposed by Gatheral, Jaisson, and Rosenbaum (2018).

We discuss the main features of the process and propose different estimation procedures to identify its parameters. The first is a two-step method that takes as given the parameters governing the univariate marginals, whereas the second method identifies all the parameters at once. Regarding the first method, the estimation of the univariate process is well documented in the literature, but it often presents a bias in the speed of mean-reversion parameter, an issue that we try to overcome. We derive the asymptotic properties of the estimators and compare the performance of the finite sample behavior to the asymptotic
theory with Monte Carlo experiments.

Finally, an empirical investigation is carried out on 7 realized volatility time series overlapping in time and available on a long sampling period. Our results show how realized volatility time series are strongly correlated and present different degrees of asymmetry in their cross-covariance structure, which can be linked to what are known as spillover effects. These features can all be well captured by our model. Moreover, in accordance with the existing literature, we observe behaviors close to nonstationarity and roughness in the trajectories. A forecasting exercise provides further evidence of the consistency of our
model with observed data.