We propose a family of copula-based multivariate distributions with g-and-h marginals. After studying the properties of the distribution, we develop a two-step estimation strategy and analyze via simulation the sampling distribution of the estimators. The methodology is used for the analysis of a seven-dimensional data set containing 40 871 operational losses. The empirical evidence suggests that a distribution based on a single copula is not flexible enough, and thus we model the dependence structure by means of vine copulas. We show that the approach based on regular vines improves the fit. Moreover, even though losses corresponding to different event types are found to be dependent, the assumption of perfect positive dependence is not supported by our analysis. As a result, the value-at-risk of the total operational loss distribution obtained from the copula-based technique is substantially smaller at high confidence levels with respect to the one obtained using the common practice of summing the univariate value-at-risks.