SummaryThe simultaneous analysis of several financial time series is salient in portfolio setting and risk management. This paper proposes a novel alternating expectation conditional maximisation (AECM) algorithm to estimate the vector autoregressive moving average (VARMA) model with variance gamma (VG) error distribution in the multivariate skewed setting. We explain why the VARMA‐VG model is suitable for high‐frequency returns (HFRs) because VG distribution provides thick tails to capture the high kurtosis in the data and unbounded central density further captures the majority of near‐zero HFRs. The distribution can also be expressed in normal‐mean‐variance mixtures to facilitate model implementation using the Bayesian or expectation maximisation (EM) approach. We adopt the EM approach to avoid the time‐consuming Markov chain Monto Carlo sampling and solve the unbounded density problem in the classical maximum likelihood estimation. We conduct extensive simulation studies to evaluate the accuracy of the proposed AECM estimator and apply the models to analyse the dependency between two HFR series from the time zones that only differ by one hour.