Risk aggregation is a popular method used to estimate the sum of a collection of financial assets or events, where each asset or event is modelled as a random variable. Applications include insurance, operational risk, stress testing and sensitivity analysis. In practice, the sum of a set of random variables involves the use of two well-known mathematical operations: n-fold convolution (for a fixed number n) and N-fold convolution, defined as the compound sum of a frequency distribution N and a severity distribution, where the number of constant n-fold convolutions is determined by N, where the severity and frequency variables are independent, and continuous, currently numerical solutions such as, Panjer’s recursion, fast Fourier transforms and Monte Carlo simulation produce acceptable results. However, they have not been designed to cope with new modelling challenges that require hybrid models containing discrete explanatory (regime switching) variables or where discrete and continuous variables are inter-dependent and may influence the severity and frequency in complex, non-linear, ways. This paper describes a Bayesian Factorisation and Elimination (BFE) algorithm that performs convolution on the hybrid models required to aggregate risk in the presence of causal dependencies. This algorithm exploits a number of advances from the field of Bayesian Networks, covering methods to approximate statistical and conditionally deterministic functions to factorise multivariate distributions for efficient computation. Experiments show that BFE is as accurate on conventional problems as competing methods. For more difficult hybrid problems BFE can provide a more general solution that the others cannot offer. In addition, the BFE approach can be easily extended to perform deconvolution for the purposes of stress testing and sensitivity analysis in a way that competing methods do not.