We use recent advances in polynomial diffusion processes to develop a continuous-time joint mortality model for the actuarial valuation and risk analysis of life insurance liabilities. The model considers the stochastic nature of future mortality improvements and introduces a common subordinator for the marginal survival processes, resulting in a nontrivial dependence structure between the survival of pairs of individuals. Polynomial diffusion processes can be used to derive closed-form formulae for standard actuarial quantities. The model fits well with a classic dataset provided by a Canadian insurer and can be used to evaluate products issued to multiple lives, as shown through numerical applications.