Bayesian predictive inference for the aggregate claims distribution and ruin probability

Compound process models are widely used to model the sum of claims against an insurance portfolio or total withdrawals from a bank. The classic compound Poisson-Exponential model is often used, and the standard approach estimates unknown parameters from historical data using maximum likelihood estimation. Unfortunately, this method can significantly underestimate risk. We use a Bayesian approach to compute a predictive distribution for cumulative loss and the ruin probability (i.e., bankruptcy) and obtain closed-form, analytical solutions. We examine the performance of these methods and show, under a range of sample sizes, that our method outperforms traditional approaches, and we provide statistical evidence that it achieves zero coverage error, improving risk assessment.

The compound process model is the workhorse in the insurance industry for modeling future claims where the number and the size of the claims are both random [1]. Of practical interest is the distribution of the aggregate claims at some time T, and the probability that the aggregate claims lead to bankruptcy, or ruin. The standard is the Cramer-Lundberg model where claims arrive according to a Poisson process and, independently, claim sizes have a known distribution [2]. Often the claim-size distribution is assumed to be exponential, which leads to a Poisson-Exponential mixture model. Beyond actuarial science, the Poisson-Exponential compound process appears in diverse fields, including stochastic hydrology [3], reliability engineering [4], and queueing systems [5].

The allure of the Cramer-Lundberg model is that the aggregate claim distribution can, at least in theory, be written in a simple form. However, numerical techniques are still needed to compute it. In the case of exponential claim sizes, there is also a closed-form solution to the so-called ruin probability. Exponential claim sizes are used to model thin-tailed insurance claims, such as auto and basic medical or property claims, but can be combined with an independent heavy-tailed process to model more extreme claims [1, 6]. In this study, we analyze the Poisson-exponential model from a Bayesian perspective to make predictive inferences for the aggregate claims and ruin probability and this limits our immediate results to claim severity distributions that have exponential tails.

The Bayesian approach to modeling future claims and its benefits are not new, Cairns [7] and recent papers highlight that Bayesian methods have become more mainstream in predictive risk modeling [8, 9]. The influential research of Dickson et al. [10] used a Bayesian approach to predict future aggregate claims and highlighted the significant risks posed by traditional methods. However, as we discuss, their approach did not fully capture the benefits of applying Bayesian methods. Since then, work in this area has focused on using numerical methods such as Markov Chain Monte Carlo (MCMC) simulations to model future claims [11–13] under various claim-size distributions. Most recently, Gómez-Déniz et al. [14] modeled the ruin probability using a Bayesian approach for different claim-size distributions, but with the prior restriction that the premiums were sufficient to avoid guaranteed ruin. Bayesian methods have also been used to model queue waiting time [15], which is applicable given the duality with aggregate claims.

Our research has extended this body of work. First, we derived an analytic expression for the aggregate claims distribution by fully marginalizing the weighted convolution of claim sizes, or the conditional aggregate claim distribution. This contrasts with Dickson et al. [10], which used predictive densities as input to traditional numerical methods for computing aggregate claims. Their approach, while an improvement over the so-called plug-in method, does not achieve zero-coverage error (ZCE), an important performance criterion. Our full marginalization of the compound distribution is a key methodological advance that fully incorporates parameter uncertainty, leading to improved prediction coverage. We demonstrate the improvement of our method via simulations and give initial statistical evidence that the Bayesian approach obtains ZCE.

Second, we derived an analytic expression for the ruin probability by marginalizing the known closed-form expression over the model parameters. In contrast to Gómez-Déniz et al. [14], we fully account for the possibility that the chosen safety loading, or equivalently the premium rate, may be insufficient to avoid bankruptcy. This static part of the ruin probability is independent of the capital level and can be significant. We illustrate the credibility of our approach using simulations highlighting the significant improvement in risk assessment vs. other methods. In general, the traditional approaches underestimate risk, and the improved results from the Bayesian method should translate directly into improved capital management decisions, such as premium rate setting and reserve determination.

This study is organized as follows. In Section 2, we derive analytic expressions for both the aggregate claims distribution and the ruin probability by fully marginalizing the model parameters. In Section 3, we illustrate the benefit of our approach via a number simulations. We show the plug-in approach can lead to severe underestimates of risk, even for a large amount of historical data. We also provide preliminary statistical evidence that our approach achieves ZCE, which is important for estimating quantiles. Lastly, we illustrate that the Bayesian approach provides a more accurate picture of the ruin probability and compare the results via simulations.

Aggregate claims

Let X 1, ⋯, X n ≡ X 1: n denote the past claims over a time frame of length t and Y 1, ⋯, Y N ≡ Y 1: N denote the future, unknown claims for a time frame T where we have suppressed the dependence of n and N on t and T. We assume the claim sizes X 1: n and Y 1: N are independent with common distribution function f (τ) = α e −ατ and their counts (n, N) are assumed to follow a Poisson process with intensity λ independent of the claim sizes. Both α and λ are unknown.

The aggregate-claims is a Poisson-Exponential (Poi(λ) Exp (α)) compound process given by

and the distribution of the aggregate claims, at time T, can be written as

where p k = P (N = k) are the Poisson probabilities, f is the underlying claim-size density and the exponent k * denotes the k -fold convolution [2], with f 0* (y)≡δ(y) being the Dirac delta function. For an insurance company, the aggregate-claims distribution is a key component used in setting premiums and reserve capital levels.

Using the probability generating function (pgf) for the Poisson distribution, one method for computing the aggregate-claims distribution [16] is to evaluate the pgf at the Laplace transform of the claim-size density and then compute the inverse Laplace transform (L - 1)

In practice, we obtain a numerical solution by discretizing the claim-size density and using the fast Fourier transform (FFT) [17] or, alternatively, Panjer's recursion [18]. Equations 2, 3 are the conditional compound distribution given that the parameters of the size and count densities are known.

A traditional method for determining a predictive distribution is to use the historical data, X 1: n, to estimate the parameters and plug them into Equation 3. This method is often called the plug-in approach. For example, the maximum likelihood (ML) estimates for the Poisson parameter, λ, are λ ^ M L = n / t, and the exponential parameter, α, is α ^ M L = 1 / X ̄, with X ̄ being the sample mean of X 1: n.

Although the plug-in method may be reasonable for large n, it does not incorporate uncertainty in the model parameters. The Bayesian method, however, is ideal for this [7], and the Bayesian predictive distribution conditioned on the past data is

where f (λ, α| x 1: n) is the posterior of the model parameters given the observed data and y is independent of x 1: n.

Once the predictive distribution is determined, statistics of interest, such as quantiles, are computed. The 100 × q 'th quantile, η q is defined by

and the predictive distribution is used to estimate η ^ q. Given the important role of quantiles, a key performance metric is the coverage error [19]. Given a quantile estimate formed from the predictive compound distribution, η ^ q = F ^ S (T) - 1 (q | x 1: n), a desirable feature is that the estimate provides zero coverage error (ZCE) with P (S (T) > η ^ q) = 1 - q, where P is the joint probability distribution of { X 1: n, Y 1: N }. Another definition of ZCE, and one we will find useful for testing, is that F ^ S (T) (X 1: n, Y 1: N) should be a uniform random variable [20].

As mentioned, one approach is to utilize the ML estimates (λ ^ M L, α ^ M L), as if they were the true parameters, and apply one of the standard numerical methods to compute a predictive distribution from Equation 2. We call this the ML, or plug-in, method. Another approach in Dickson et al. [10] computes the predictive distribution of claim sizes and counts and uses it in Equation 2 to obtain a predictive distribution, which we call the DTZ method. The DTZ method has lower coverage error than the ML method, but neither yields ZCE.

The predictive claim-size distribution, f (y | x 1: n), is computed as

with posterior f (α| x 1: n)∝ f (x 1: n |α) f (α). If a conjugate prior distribution from the Gamma family is used, α~Γ(a, b) with b the rate parameter, the resulting posterior is also Gamma, α| x 1: n ~Γ(n + a, Σ+ b), where Σ ≡ Σ i = 1 n x i. Substituting into Equation 6 results in mixing an exponential distribution with a gamma distribution and produces

which is a Pareto(n + a, Σ+ b) distribution. Compared to the exponential distribution, used in the plug-in method, this is a heavier-tailed claim-size distribution.

Similarly, the predictive count distribution, P (k | x 1: n), is computed as

with posterior f (λ| n)∝ f (n |λ) f (λ). If a conjugate prior distribution from the Gamma family is used, λ ~ Γ (ã, b ~), the resulting posterior is also Gamma, λ | n ~ Γ (n + ã, t + b ~). Substituting into Equation 8 results in mixing a Poisson distribution with a gamma and produces a negative binomial predictive count distribution

As with the predictive claim-size distribution, the negative binomial has more dispersion than the Poisson used in the plug-in method. In addition to being a conjugate prior, Gamma priors are motivated by the fact that Jeffreys priors [21], an invariant method of prior construction, are in the Gamma family. In particular, the Jeffreys priors are f (α)∝1/α, or α~Γ(0, 0), and f (λ) ∝ 1 / λ although Jeffreys argued for using f (λ)∝1/λ [22]. As such, we will adopt the improper Γ(0, 0) priors along with the assumption that P ({ n = 0}∪{Σ = 0}) = 0 insuring proper posteriors.

Armed with the predictive distributions in Equations 7, 9, Dickson et al. [10] utilize them to implement Equation 3 to compute a predictive distribution. In essence, the DTZ method is based on Equation 2 where the densities are replaced by their predictive counterparts,

where the hat, ˆ, indicates a predictive density with parameters marginalized. The term, f ^ k * (y) is the k -fold convolution of the predictive size density. As shown in Appendix A, however, the Bayesian method would be to marginalize the compound process, which results in

where we marginalize the k -fold convolution of the claim-size density and, in general, the two approaches are not equal. The predictive count distribution p ^ k remains a negative binomial as in Equation 9. The k -fold convolution of the claim-size density is the density of the sum of k independent sizes, Y k, and for an underlying exponential with k >0 is Y k |α, x 1: n ~Γ(k, α) and the predictive posterior is, therefore, a Γ(k, α) Γ(n, Σ) mixture

Performing this integration results in a scaled beta prime distribution for k >0

where B (k, n) is the standard beta function and f 0 * ^ (y) ≡ δ (y). Using Equation 13 in Equation 11 shows the continuous portion of the Bayesian predictive compound distribution is a mixture of a negative binomial and a generalized beta prime distribution, and we can write

which can be written as

where 2 F 1 (·) denotes the hypergeometric function [23] (See Appendix B for details).

The ruin probability

In addition to computing the aggregate-claims distribution, for a time T in the future, an important risk measure for setting reserve capital levels is the ruin probability. The risk-process is defined as the amount of capital at time T given initial capital u and premium rate c

where S (T) are the aggregate claims at time T.

The probability of ruin at time T is Ψ (T, u) = P (U (t, u) < 0 | t ≤ T) and ultimate ruin is Ψ(u) = Ψ(∞, u). We will simply refer to Ψ(u) as the ruin probability. Solvency (Ψ(u) < 1) requires the safety loading, ρ = c α λ - 1 > 0, which is the percent of premiums collected over the average of claims paid per unit time.

For the Poi(λ)Exp(α) process the ruin probability can be written in closed form [24] as

The plug-in method uses λ ^ M L and α ^ M L in Equation 17 with c = (1 + ρ ^) λ ^ M L / α ^ M L based on a chosen ρ ^ > 0. We use the notation ρ ^ since the true values of λ and α are unknown and, therefore, the true safety loading is unknown. There is always the chance that c is too low and ρ is actually ≤ 0, implying Ψ(u |λ, α) = 1.

The Bayesian ruin probability is derived by marginalizing the conditional distribution over the posterior.

The integral is the sum of two terms. One for α−λ/ c ≤ 0 which implies Ψ(u |λ, α) = 1 and the other for α−λ/ c >0. The first term, which we call the static ruin probability, for ρ ≤ 0, is

and is that part of the ruin probability due to choosing too low a c, independent of u. Since Ψ(u |λ, α) = Ψ α (u |α)Ψ λ (u |λ) is separable, the second term can be written

These two terms, with ρ ~ = 1 + ρ ^, can be written (See Appendix C)

and are summed to compute the ruin probability, Ψ(u | x 1: n), for n >1.

We first tested the performance of the different methods using simulations to estimate the aggregate-claims distribution over a time horizon T given past claim experience over a time period t and, second, to estimate the ruin probability. In particular, because an insurance company would be interested in estimating a high quantile to help set premiums and capital reserves against future losses, we focus on quantile estimation performance. We chose α = 0.2 and λ = 5 so that, on average, five claims occur per unit time period, with an average claim size of 5. The parameter choice is most applicable for a large insurance portfolio where p 0 = P (N = 0)≈0. We chose different combinations of t and T ranging from 10 to 100 and, for each combination, ran 10,000 simulations.

To begin, we consider the plug-in method, where for each simulation run we estimate the parameters via maximum likelihood and then use those estimates to directly estimate the aggregate claim distribution via the FFT method [17]. Specifically, once α ^ M L is estimated, the exponential density is sampled, and its Fourier transform is computed. We sampled the density at a fine spacing of λ T /1600 with 2 15 samples to ensure accuracy and avoid any numerical/aliasing issues at the cost of computational efficiency. This transform is input to the pgf of the Poisson distribution, along with λ ^ M L, and then the inverse Fourier transform is performed. Once the aggregate claim distribution is estimated, we calculate quantiles which we can compare to the actual quantiles. For each simulation run, we generate future claims for time period T, and then compute the percentage of simulation runs that result in the quantile estimate being exceeded.

Figure 1 shows the performance of the plug-in method for t = T = 10. This is a fairly stressful scenario where we have only 50 past claims, on average, to estimate the total of 50 future claims, on average. On the left, we show the actual quantile level along with the average estimate from the simulation, and the ML quantile estimates appear unbiased. On the right, we show the percent of the simulation runs that produced a future sample in one of the ten 1% prediction intervals (quantile buckets) from 90-91% to 99-100%. If the method is ZCE, this should result in a uniform distribution (i.e., the horizontal line in the graph). While the ML estimate appears unbiased, the plug-in approach shows significant coverage error, particularly for high quantiles. Importantly, the coverage error is positive (i.e., over), and the plug-in method underestimates the quantile, resulting in excess exceedances. In Johnston [25], it was shown that, under mild conditions on the underlying cdf, an unbiased quantile estimator yields coverage error.

To see the impact of additional past claims on estimating future aggregate claims and coverage error, we varied t from 10 to 100 for T = 10 and T = 25. Figure 2 shows the total coverage error above the 90% quantile, and as t increases, the coverage error declines. That said, even for t = 100, or 500 past claims on average, there is still about 1% (over) coverage error when predicting aggregate claims over T = 10, or 50 claims on average. As expected, increasing T while holding t constant increases the coverage error. We noted that the coverage error was approximately constant as a function of t / T, because the predictive counting process for the plug-in method is Poi(nt / T). The ML method requires a significant amount of past data, with the ratio of t / T >10 to achieve minimally acceptable coverage error.

Figure 3 shows the coverage probabilities for all three methods for t = T = 10. We can see that the DTZ approach, by including parameter uncertainty, reduces the coverage error although it does not attain ZCE. The Bayes approach, with marginalized convolution, appears to achieve ZCE even with the small sample size of t = 10. While ZCE can be shown for an exponential distribution using Bayes methods [25], the compound-distribution case is more difficult. In particular, Lawless and Fredette [26] show that for the discrete Poisson distribution a pivotal quantity does not exist, so that it is merely asymptotically ZCE.

We illustrate that the Bayesian method is potentially ZCE using a statistical argument. Given the discretization shown in Figure 3, we can test if the quantized data comes from a discrete uniform distribution using a Chi-Squared test. With J quantile buckets with M j samples, and M total samples within the J buckets, we compute p ^ M (j) = M j / M and the statistic

which, under the uniform hypothesis (H 0), is approximately chi-square with J −1 degrees of freedom. For the simulation shown in Figure 3, the p -values were virtually zero for the ML and DTZ methods, indicating we would readily reject the uniform hypothesis and, therefore, ZCE. Conversely, the p -value for our Bayes approach was 0.62, highlighting consistency with ZCE. We emphasize that the failure to reject the uniform hypothesis, via the chi-square test, does not prove ZCE and that care needs to be taken using a p -value, P (D ≥ D obsv / H 0), to validate the ZCE hypothesis.

Shown in the Table 1 below, for different t and T = 10, are the mean p -values and the number of p -values below 5% for each method for which we ran the full simulation (10,000 future claim samples) 100 times. The mean p -value for the Bayes method was consistently 0.5, and the number of times ZCE was rejected at the 5% level was consistent with the Chi-Squared test conditioned on ZCE being true. Conversely, the ML method performs poorly. Even with t = 100, rejection occurs 71 out of 100 times. The DTZ method improves but still provides no convincing evidence of ZCE.

These results provide initial evidence that by not using the Bayesian method, unless t >> T, high quantiles for the aggregate claims can be underestimated using the plug-in or DTZ methods. While additional simulation studies are needed, varying model parameters and underlying distributions, the Bayesian method appears to offer promise and potentially be ZCE.

The ruin probability can also be significantly understated when parameter uncertainty is ignored. The plug-in method uses λ ^ M L and α ^ M L, in Equation 17 and c, the premium rate, is chosen based on a desired ρ ^ > 0 with c = (1 + ρ ^) λ ^ M L / α ^ M L. However, since we do not know the true λ, α, the resulting c may be too low and, if so, the actual ρ is ≤ 0. Incorporating this risk is a key improvement of our approach, compared to Gómez-Déniz et al. [14], leading to improved capital management decisions.

Figure 4 shows the static ruin probability, Ψ (ρ ≤ 0 | ρ ~), ρ ~ = 1 + ρ ^, which is independent of the initial capital u, for n = 500. If a ρ ^ = 10 % was chosen, the static ruin probability is 6.6%, which would never dissipate, and a ρ ^ ≥ 25 % would be required to make the static ruin probability ≈0. In other words, even with, on average, 500 past claims to estimate future claim experience (e.g., λ = 5, t = 100), we would need to charge a premium rate at least 25% higher than the historical claim rate. In some applications, using only 500 past claims may seem restrictive, although for stop-loss insurance or modeling nonstationarity this would be a generous sample size.

Figure 5 shows the total ruin probability for two safety loading settings, ρ ^ = 10 % and ρ ^ = 20 %, where we chose t = 100, n = 500 and a premium rate c = 10%. For the plug-in method, this would result in λ ^ = n / t = 5 and α ^ = (1 + ρ ^) λ ^ / c = 50 (1 + ρ ^) and we can see the significant gap between the ruin probability estimate using the plug-in vs. the Bayesian method. For example, with a chosen safety loading of ρ ^ = 20 % and initial capital u = 100, one would be inclined to believe the ruin probability is virtually zero with a large margin for error. However, the Bayesian method indicates the ruin probability is close to 1%.

We also performed an extensive simulation of future claims and computed an empirical estimate for the ruin probability. We used a cutoff T = 10000, computed U (t /λ, u), t = 1, ⋯λ T to ensure a fine sampling of the capital level, and performed 5000 simulation runs to compute an empirical ruin probability. The simulation results are also shown in Figure 5 as the asterisks (*), and the empirical estimates are nearly identical to the Bayesian prediction. While additional simulation studies are required, these initial results indicate the promise in using the Bayesian methodology and the potential improvement in capital reserve and premium rate setting for an insurance company.

In this research, we derived a complete analytic Bayesian solution for determining the aggregate claim distribution and the ruin probability for the Poisson-exponential compound process. This improves on traditional and other approaches and is important for risk assessment, as well as for properly setting capital reserves and premium rates to avoid bankruptcy. We illustrate the benefits of a Bayesian methodology via simulations, which also highlight the inherent risks of the traditional plug-in approach.

There are a number of avenues for future research. First, additional simulations and analytical results are needed to satisfactorily state that the Bayesian method is asymptotically ZCE. This includes not only varying the model parameters but also deviations from the Poisson-exponential process. Comparing the numerical efficiency and stability of the methods is important, particularly for a large number of claims. Ultimately, using real claim data would be useful, and comparing our method's implications for capital reserving and premium rates would be valuable. Finally, we look forward to extending this research by including heavier-tailed claim-size distributions and allowing for dynamic model parameters.