We introduce a compound multivariate distribution designed for modeling insurance losses arising from different risk sources in insurance companies. The distribution is based on a discrete-time Markov Chain and generalizes the multivariate compound negative binomial distribution, which is widely used for modeling insurance losses.We derive fundamental properties of the distribution and discuss computational aspects facilitating calculations of risk measures of the aggregate loss, as well as allocations of the aggregate loss to individual types of risk sources. Explicit formulas for the joint moment generating function and the joint moments of different loss types are derived, and recursive formulas for calculating the joint distributions given. Several special cases of particular interest are analyzed. An illustrative numerical example is provided.