Loss Distribution

First, let us assume we have perfect information, i.e. we know the loss distribution for a certain risk.

To keep it simple, we assume the loss distribution is log-normal with a mean (M) of 65 and a coefficient of variation (CV) of 30%. The corresponding log-normal parameters μ and σ can be derived σ² = log (1 + CV²),  μ = log (M) − σ²/2.

The following chart shows the corresponding probability density curve f(x), cumulative distribution function F(x) and survival function S(x) = 1 − F(x).

In the insurance context, the survival function is
often called exceedance probability function, as it describes the probability of exceeding a certain loss.

Expected Loss

Let’s define X as the random variable that describes the ground up loss. The expected loss (in plain English the probability weighted sum of the losses) of a positive random variable is then E[X] = ∫₀^∞x f(x)dx.

Limited Expected Loss

Let’s further assume that losses are limited by α in the contract. Then the limited expected value (LEV), written as E[X ∧ α], is given as E[X ∧ α] = ∫₀^αx f(x)dx + ∫_α^∞α f(x)dx.

To evaluate this sum of integrals we recall that F(x) = ∫f(x)dx and F(∞) = 1. Hence ∫_α^∞α f(x)dx = α∫_α^∞f(x)dx = α − αF(α).

The integration by parts theorem helps with the first part of the integral ∫₀^αx f(x)dx = αF(α) − ∫₀^αF(x)dx.

Therefore, we retrieve the well-known identity E[X ∧ α] = α − ∫₀^αF(x)dx = ∫₀^α1 − F(x)dx = ∫₀^αS(x)dx.

The integral describes the area above the CDF up to α, or the area under the survival function up to α. In the example below the limit α was set to 100.

Loss Cost of a Layer

Suppose we want to insure the layer of claims between 80 and 100.

The expected loss cost would be the difference between the limited expected values of E[X ∧ 100] − E[X ∧ 80] and can be evaluated numerically

sigma <- sqrt(log(1 + 0.3^2))
mu <- log(65) - sigma^2/2
S <- function(x){ 1 - plnorm(x, mu, sigma) }
(lyr <- integrate(S, 80, 100)$value)

## [1] 2.22814

Increased Limit Factor

On the other hand the ratio of the original loss cost with a limit of 100 to a limit of 80 is called an increased limit factor (ILF):

(ILF <- integrate(S, 0, 100)$value / integrate(S, 0, 80)$value)

## [1] 1.03592

Therefore we would expect the LEV to increase by 3.5% as the limit increases from 80 to 100.

Note ILFs are often used for pricing casualty business.

Exposure curves

From the loss distribution we could read of the expected loss for any layer. However, often we will not have that level of information about the risk.

Instead, we will have to infer information from others risks that share similar loss characteristics as a function of the overall exposure, assuming that the relative loss size distribution is independent of the individual risk characteristic.

Hence, we will require a view on the expected full value loss cost for the overall exposure.

MPL <- 200
ExpectedLoss <- 65
Deductible <- seq(0, MPL, 1)
G <- sapply(Deductible, function(x){
  LEVx <- integrate(S, 0, x)$value
  LEVx/ExpectedLoss
})
plot(Deductible/MPL, G, 
     t="l", bty="n", lwd=1.5,
     main="Exposure Curve",
     xlab="Deductible as % of MPL",
     ylab="% of expected loss paid by insured",
     xlim=c(0,1), ylim=c(0,1))
abline(a=0, b=1, lty=2)

Classic distributions

As already mentioned, the exposure curve function of X is defined as the ratio of the limited expected value and the expectation. Since the loss X is a positive random variable, the exposure curve is $$ G_X(d) = \frac{\int_0^d (1-F_X(x))dx }{\int_0^1 (1-F_X(x))dx }. $$ Note that the exposure curve is a concave function for d ∈ ]0, 1[.

There is a direct link between the distribution function and the exposure curve. Since $$ F_X(x) = \left(1- \frac{G_X'(x)}{G_X'(0)}\right)1\!\!1_{[0,1[}(x) + 1\!\!1_{[1,+\infty[}(x), $$ defining the exposure curve or the distribution function is equivalent. The exposure curve is also a concave increasing function, see e.g. (Antal 2003).

One-inflated distributions

The MBBEFD distribution, first parametrization

We denote this first parametrization by MBBEFD(a, b). We define the parameter domain 𝒟_a, b as Let us note that this domain includes two particular sets 𝒟_a, 1 = {(a, 1), a + 1 > 0} and 𝒟_0, b = {(0, b), b > 0}.

The MBBEFD distribution, second parametrization

Fitting methods: MLE and TLMME

We consider two fitting methods, namely the maximum likelihood estimation (MLE) and the total-loss-moment matching estimation. Both methods are implemented in the fitDR function of the package.

As its name suggest, the maximum likelihood function consists in maximizing the likelihood function defined as $$ \mathcal L(\theta, x_1,\dots,x_n)= \prod_{i=1}^n f_X(x_i), $$ where f_X denotes the density of the loss random variable X and θ is the parameter vector. When X is a continuous non-inflated distribution, f_X is the usual density function (obtained by differentiating the distribution function) when X is a mixed-type distribution with a non-null probability P(X = 1), f_X is the product of P(X = x)1  1_x = 1 and F_X′(x)1  1_x ≠ 1.

For example when considering the beta distribution, the likelihood function to be maximized is $$ \mathcal L((a,b), x_1,\dots,x_n)= \prod_{i=1}^n \frac{x_i^{a-1}(1-x_i)^{b-1}}{\beta(a,b)}. $$

The total-loss-moment matching estimation consists in matching the total loss and the first moments of a distribution. That is we solve the system $$ \left\{\begin{array}{ll} P(X = 1;\theta) = tl_n \\ E(X; \theta) = \bar x_n\\ \vdots \\ E(X^k; \theta) = m_{n,k}\\ \end{array}\right. \text{ where } \left\{\begin{array}{ll} tl_n = \frac{1}{n}\sum_{i=1}^n 1\!\!1_{x_i =1},\\ \bar x_n = \frac{1}{n}\sum_{i=1}^n x_i,\\ m_{n,k} = \frac{1}{n}\sum_{i=1}^n (x_i)^k. \end{array}\right. $$ The number of equation k + 1 equals the number of parameters. This methods is relevant only for one-inflated distributions as for classic continuous distributions such as the beta distribution we know that P(X = 1) = 0.

The function of the package returns an object of class inheriting from the class for which generic functions are available . There is also a bootstrap procedure implemented in the function.

Examples

Let us fit some distributions with a simulated dataset based on a mixture of two distributions (regular beta and mbbefd).

set.seed(123456)
x <- c(rbeta(50, 3, 1/2), rmbbefd(50, 1/2, 1/10))
f1 <- fitDR(x, "mbbefd", method="mle")
summary(f1)

## Fitting of the distribution ' mbbefd ' by maximum likelihood 
## Parameters : 
##      estimate Std. Error
## a 0.035909278         NA
## b 0.009872438         NA
## Loglikelihood:  -39.69145   AIC:  83.3829   BIC:  88.59324 
## Correlation matrix:
## [1] NA

b1 <- bootDR(f1, niter=20)
summary(b1)

## Parametric bootstrap medians and 95% percentile CI 
##       Median        2.5%      97.5%
## a 0.03603436 0.012784198 0.06813441
## b 0.00930872 0.002569712 0.02273344

The graph correspond respectively to the histogram (and empirical density) against the fitted density and a two-dimensional level curve and scatter plot for bootstrap parameters.

As in the package, we can compare multiple fits on the same graph or with the same statistics. We just have to pass a list of objects.

f2 <- fitDR(x, "oibeta", method="mle")
f3 <- fitDR(x, "oiunif", method="mle")
gofstat(list(f1, f2, f3))

## Goodness-of-fit statistics
##                              1-mle-mbbefd 2-mle-oibeta 3-mle-oiunif
## Kolmogorov-Smirnov statistic    0.1400000    0.1400000    0.2154661
## Cramer-von Mises statistic      0.2715687    0.1414871    1.7856663
## Anderson-Darling statistic            Inf          Inf          Inf
## 
## Goodness-of-fit criteria
##                                1-mle-mbbefd 2-mle-oibeta 3-mle-oiunif
## Akaike's Information Criterion     83.38290     41.84037     82.99270
## Bayesian Information Criterion     88.59324     49.65588     85.59787

par(mfrow=c(1, 2))
cdfcomp(list(f1, f2, f3), leg=c("mbbefd", "oibeta", "oiunif"))
denscomp(list(f1, f2, f3), leg=c("mbbefd", "oibeta", "oiunif"), 
         ylim=c(0,4), xleg="topleft")

Of course, we can also compare the different fitted exposure curves using the function . That is

par(cex=0.8)
eccomp(list(f1, f2, f3), leg=c("mbbefd", "oibeta", "oiunif"), do.points=FALSE)

Computation of premium rate (example)

The following example data is taken from the earlier Swiss Re technical paper above. It is the aim to find the risk premium for a per risk WXL cover (fire only) with the limits CHF 3, 500, 000 xs CHF 1, 500, 000.

The total gross net premium income (GNPI) is CHF 95, 975, 000 in 2004, with an expected loss ratio of 55%. Therefore we have a view on the mean expected loss cost.

The risk profile given by the cedant is from the year 2002. However, instead of re-indexing the historical data to 2004, we back-index the data to 2002 with a factor of 457/550 = 0.83:

CHF 2, 908, 182 xs CHF 1, 246, 364

The data presented below is from the client, with pre-selected exposure curves, in this case the parameter c of the MBBEFD curve.

The policies have been grouped into different exposure bands. The mean MPL is simply the average of the lower and upper band.

In order to use the exposure curves we need the expected loss and the deductible as % of the MPL.

We consider three scenarios in the following subsections.

Premium_rate <- function(Deductible, Limit, MaxMPL, 
                   MeanMPL, GrossPremium, C, ULR){
  DedPerMPL <- Deductible / ifelse(MaxMPL < Deductible, Deductible, 
                                   ifelse(MaxMPL<Limit, MeanMPL, Limit))
  LossShare <- 1 - apply(cbind(DedPerMPL, C), 1, 
                         function(x){ 
                           ecMBBEFD(x[1], 
                                          b=swissRe(x[2])['b'], 
                                          g=swissRe(x[2])['g'])
                           })
  NetPremium <- GrossPremium * ifelse(MaxMPL < Limit, 1, 
                                      Limit/MaxMPL)
  XL_Premium <- NetPremium * LossShare
  # Rate on Line
  sum(XL_Premium)/sum(NetPremium)*ULR
}

pr <- Premium_rate(Deductible = D,
                   Limit = L,
                   MaxMPL = ClientData$MaxMPL,
                   MeanMPL = ClientData$MeanMPLGrossLoss,
                   GrossPremium = ClientData$GrossPremium,
                   C=ClientData$ExposureCurve, ULR=0.55)