Polynomial approximations for bivariate aggregate claims amount probability distributions

A numerical method to compute bivariate probability distributions from their Laplace transforms is presented. The method consists in an orthogonal projection of the probability density function with respect to a probability measure that belongs to a Natural Exponential Family with Quadratic Variance Function (NEF-QVF). A particular link to Lancaster probabilities is highlighted. The procedure allows a quick and accurate calculation of probabilities of interest and does not require strong coding skills. Numerical illustrations and comparisons with other methods are provided. This work is motivated by actuarial applications. We aim at recovering the joint distribution of two aggregate claims amounts associated with two insurance policy portfolios that are closely related, and at computing survival functions for reinsurance losses in presence of two non-proportional reinsurance treaties