  In this paper we introduce a bivariate distribution on $\mathbb{R}_{+} \times
\mathbb{N}$ arising from a single underlying Markov jump process. The marginal
distributions are phase-type and discrete phase-type distributed, respectively,
which allow for flexible behavior for modeling purposes. We show that the
distribution is dense in the class of distributions on $\mathbb{R}_{+} \times
\mathbb{N}$ and derive some of its main properties, all explicit in terms of
matrix calculus. Furthermore, we develop an effective EM algorithm for the
statistical estimation of the distribution parameters. In the last part of the
paper, we apply our methodology to an insurance dataset, where we model the
number of claims and the mean claim sizes of policyholders, which is seen to
perform favorably. An additional consequence of the latter analysis is that the
total loss size in the entire portfolio is captured substantially better than
with independent phase-type models.
