Abstract-Multifractal (MF) analysis enables the theoretical study of scale invariance models and their practical assessment via wavelet leaders. Yet, the accurate estimation of MF parameters remains a challenging task. For a range of applications, notably biomedical, the performance can potentially be improved by taking advantage of the multivariate nature of data. However, this has barely been considered in the context of MF analysis. This paper proposes a Bayesian model that enables the joint estimation of MF parameters for multivariate time series. It builds on a recently introduced statistical model for leaders and is formulated using a 3D gamma Markov random field joint prior for the MF parameters of the voxels of spatio-temporal data, represented as a multivariate time series, that counteracts the statistical variability induced by small sample size. Numerical simulations indicate that the proposed Bayesian estimator significantly outperforms current state-of-the-art algorithms.