This paper develops an unbiased Monte Carlo estimator of the transition density of a multivariate jump-diffusion process. The drift, volatility, jump intensity, and jump magnitude are allowed to be state-dependent and non-affine. It is not necessary that the volatility matrix can be diagonalized using a change of variable or change of time. Our density estimator facilitates the parametric estimation of multivariate jump-diffusion models based on discretely observed data. Under conditions that can be verified with our density estimator, the parameter estimators we propose have the same asymptotic behavior as maximum likelihood estimators as the number of data points grows, but the observation frequency of the data is kept fixed. In a numerical case study of practical relevance, our density and parameter estimators are found to be highly accurate and computationally efficient.