In this paper, we introduce a class of multivariate Erlang mixtures and present its desirable properties. We show that a multivariate Erlang mixture could be an ideal multivariate parametric model for insurance modeling, especially when modeling dependence is a concern. When multivariate losses are governed by a multivariate Erlang mixture, many quantities of interest such as joint density and Laplace transform, moments, and Kendall's tau have a closed form. Further, the class is closed under convolutions and mixtures, which enables us to model aggregate losses in a straightforward way. We also introduce a new concept called quasi-comonotonicity that can be useful to derive an upper bound for individual losses in a multivariate stochastic order and upper bounds for stop-loss premiums of the aggregate loss. Finally, an EM algorithm tailored to multivariate Erlang mixtures is presented and numerical experiments are performed to test the efficiency of the algorithm.
