In spedygiorgio/mbbefd: Maxwell Boltzmann Bose Einstein Fermi Dirac Distribution and Destruction Rate Modelling
library("mbbefd")
#library("knitcitations")
#cleanbib()
#options("citation_format" = "pandoc")
#bib <- read.bibtex("mbbefd.bib")
library(lattice)
my.settings <- canonical.theme(color=FALSE)
my.settings[['fontsize']] = list(text = 8, points = 4)
my.settings[['strip.background']]$col <- "darkgrey"
my.settings[['strip.border']]$col<- "black"
Introduction and general notation
Exposure rating is a tool for insurance pricing that allocates premium to bands
of damage ratios or severity of losses. First ideas were published in [@salzmann1963]. It is often used to price non-proportional
reinsurance contracts, such as excess of loss (XL) reinsurance.
Exposure rating uses the loss experience of a similar portfolio of policies to estimate the expected losses of the portfolio to be covered. The method is frequently used as a benchmark when there is no sufficient credible claims history from the client.
In this vignette, we first present the general notation and concepts of loss modelling.
Secondly, we focus on destruction models implemented in the \pkg{mbbefd} package.
One popular model uses the MBBEFD distribution introduced by [@bernegger97].
Finally, we provide an example of pricing of XL contract.
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 $\mu$ and $\sigma$ can be derived
$$
\sigma^2=\log(1+CV^2),~~
\mu =\log(M) - \sigma^2/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)$.
n <- 100
Loss <- seq(1, 150, length=n)
mu <- 4.13
sigma <- 0.29
CDF <- plnorm(Loss, mu, sigma)
Density <- dlnorm(Loss, mu, sigma)
Survival <- 1 - CDF
dat <- data.frame(Loss=rep(Loss, 3),
Value=c(Density, CDF, Survival),
Type=gl(3, n,
labels=c("PDF: f(x)",
"CDF: F(x)",
"SF: S(x)=1-F(x)"),
ordered=TRUE))
xyplot(Value ~ Loss | Type, data=dat, ylab="",
layout=c(3,1), as.table=TRUE, type="l",
par.settings = my.settings,
par.strip.text=list(col="white", font=2),
scales=list(relation="free", alternating=1))
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] = \int_{0}^\infty x\, f(x) dx.
$$
Limited Expected Loss
Let's further assume that losses are limited by $\alpha$ in the contract. Then
the limited expected value (LEV), written as $E[X \wedge \alpha]$, is given as
$$
E[X \wedge \alpha] = \int_0^\alpha x \, f(x) dx + \int_\alpha^\infty \alpha \, f(x) dx.
$$
To evaluate this sum of integrals we recall that $F(x)=\int f(x) dx$ and $F(\infty)=1$. Hence
$$
\int_\alpha^\infty \alpha \, f(x) dx = \alpha \int_\alpha^\infty f(x) dx = \alpha - \alpha F(\alpha).
$$
The integration by parts theorem helps with the first part of the integral
$$
\int_0^\alpha x \, f(x) dx = \alpha F(\alpha) - \int_0^\alpha F(x) dx.
$$
Therefore, we retrieve the well-known identity
$$
E[X \wedge \alpha] = \alpha - \int_0^\alpha F(x) dx = \int_0^\alpha 1-F(x) dx
= \int_0^\alpha S(x) dx.
$$
The integral describes the area above the CDF up to $\alpha$,
or the area under the survival function up to $\alpha$. In the example
below the limit $\alpha$ was set to 100.
alpha <- 100
fxy <- function(x,y,...){
x1 <- as.numeric(x[x<=alpha])
if(panel.number()<2){
panel.polygon(c(x1, rev(x1)),
rev(c(rep(1, length(x1)),
rev(y[x<=alpha]))),
col="skyblue", border=NA)
}else{
panel.polygon(c(x1, rev(x1)),
c(rep(0, length(x1)),
rev(y[x<=alpha])),
col="skyblue", border=NA)
}
panel.xyplot(x,y,...)
#panel.text(x=20, y=0.4, label=paste("LEV[X]", cex=2)
}
xyplot(Value ~ Loss | Type,
data=subset(dat, !Type %in% "PDF: f(x)"),
ylab="",
panel=fxy,
layout=c(2,1), as.table=TRUE, type="l",
par.settings = my.settings,
par.strip.text=list(col="white", font=2),
scales=list(relation="free", alternating=1))
Loss Cost of a Layer
Suppose we want to insure the layer of claims between 80 and 100.
alpha1 <- 80
alpha2 <- 100
xyplot(Value ~ Loss | Type,
data=subset(dat, !Type %in% "PDF: f(x)"),
ylab="",
panel=function(x,y,...){
x1 <- c(x[x >= alpha1 & x<=alpha2])
if(panel.number()<2){
panel.polygon(
c(x1, rev(x1)),
rev(c(rep(1, length(x[x >= alpha1 & x<=alpha2])),
rev(y[x >= alpha1 & x<=alpha2]))),
col="skyblue", border=NA)
}else{
panel.polygon(
c(x1, rev(x1)),
c(rep(0, length(x[x >= alpha1 & x<=alpha2])),
rev(y[x >= alpha1 & x<=alpha2])),
col="skyblue", border=NA)
}
panel.xyplot(x,y,...)
},
layout=c(2,1), as.table=TRUE, t="l",
par.settings = my.settings,
par.strip.text=list(col="white", font=2),
scales=list(relation="free", alternating=1))
The expected loss cost would be the difference between the
limited expected values of $E[X \wedge 100] - E[X \wedge 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)
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)
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.
Normalising loss experience using TIV and MPL
To make risks more comparable we look at ratio of losses to the
underlying exposure, where the exposure is given as the
sum insured (SI), or better the total insured value (TIV),
or perhaps as the maximum probable loss (MPL).
Other definitions and measures are also popular, such as possible maximum
loss (PML), estimate maximum loss (EML), or maximum foreseeable loss (MFL).
While the TIV and SI are straightforward to understand, the other metrics can
be little more challenging to assess.
Suppose we insure a big production plant with of several buildings against fire.
It is unlikely that a fire will destroy the whole facility, instead it is
believed that the distance between the building will ensure that fires can be
contained in a local area. Hence the MPL might be only the highest value of any
of those buildings.
Shifting the metric from loss amounts to damage ratio or destruction rates (loss
as a % of the exposure metric) allows us to compare and benchmarks risks.
Analysing deductibles
Most insurance policies are written with a deductible, so that the
insured will cover the first $X of the losses.
Thus, the reinsurer needs to understand, by how much the claims burden is reduced
for a given deductible.
From the previous section we have learned how to calculate the limited expected
value, which is the loss cost to the insured.
The exposure curve (or also called first loss curve or original loss cure )
is defined as the proportion of the LEV for a given deductible $d$ compared to
the overall expected claims cost $E[X]$
$$
G(d) = \frac{E[X \wedge d]}{E[X]}.
$$
For our example, using a log-normal distribution and a MPL of 200 we get:
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)
The steepness of the curve is related to the severity of the loss
distribution. The closer the curve is to the diagonal the greater
the proportion of large loss.
If all losses were total losses then the exposure curve would be
identical with the diagonal.
Different perils and exposures will have different exposure curves.
Well known exposure curves are those used by Swiss Re and Lloyd's.
Destruction rate models
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\in]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)\ind_{[0,1[}(x) + \ind_{[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. [@antal03].
Uniform distribution
\label{sec:unif}
The most trivial example of exposure curve is obtained for the uniform distribution on $I$.
We consider $F_X(x)=x$ leading to
$$
G_X(d)
= d(2-d).
$$
Beta distribution
\label{sec:beta}
A more interesting example is obtained for the Beta distribution on $I$
(e.g.[@kotzjohnsonbalak94v2]
for which the density is $f_X(x) = x^{a-1}(1-x)^{b-1}/\beta(a,b)$ for $x\in ]0,1[$ and $a,b>0$
where $\beta(.,.)$ denotes the
beta function, see e.g.[@nist10].
The distribution function is obtained in terms of the incomplete beta ratio function
$F_X(x) = \beta(x;a,b)/\beta(a,b) = I(x;a,b)$.
We get
$$
G_X(d) \beta(a,b;x) - \frac{x^a(1-x)^b}{a\beta(a,b)} + x(1-\beta(a,b;x)) \frac{a+b}{a},
$$
where $\beta(.;.,.)$ denotes the incomplete beta function.
One-inflated distributions
\label{sec:oidistr:generic}
Characterizations
\label{sec:oidistr:generic:charac}
Let us consider a continuous distribution function $F_0$ of a random variable $X_0$.
The corresponding distribution function of the one-inflated random variable $X_1$ is
$$
F_1(x) = (1-p_1) F_0(x) + p_1 \ind_{[1,+\infty[}(x).
$$
There is no density but an improper density $(1-p)F_0'(x)$ and a probability mass $p_1$ at $x=1$.
The one-inflated beta distribution
We consider the one-inflated beta distribution.
Using Section \ref{sec:beta}, we obtain the following distribution function
$$
F_X(x) =
\systL
0 & \text{if } x<0 \
I(x;a,b)(1-p_1) & \text{if } 0\leq x <1 \
1 & \text{if } x\geq 1
\systR
$$
where $I(x;a,b)$ denotes the incomplete beta ratio function.
This leads to a non-null probability at $x=1$, $P(X=1)=p_1$.
The improper density function is
$$
\tilde f_X(x) = (1-p_1) \frac{x^{a-1}(1-x)^{b-1}}{\beta(a,b)}.
$$
The expectation is
$$
E(X)
= p_1 + (1-p_1)\frac{a}{a+b}.
$$
The exposure curve is
$$
G_X(d) =
\frac{(1-p_1) \left(1 - I(d;a,b) \frac{b}{a+b} - \frac{d^a(1-d)^b}{(a+b)\beta(a,b)}\right) + p_1d
}{p_1 + (1-p_1)\frac{a}{a+b}}.
$$
The MBBEFD distribution, first parametrization
We denote this first parametrization by $\mbox{MBBEFD}(a,b)$.
We define the parameter domain $\calD_{a,b}$ as
\begin{equation}
\calD_{a,b} = {(a,b)\in\R^2, a+1>0, a(1-b)>0, b>0}
\cup {(a,b), a=+\infty, b<1}.
\label{eq:paramdomain:mbbefd:ab}
\end{equation}
Let us note that this domain includes two particular sets $\calD_{a,1}={(a,1), a+1>0}$ and $\calD_{0,b}={(0,b), b>0}$.
Characterization by the exposure curve
The MBBEFD distribution is defined by the following exposure curve for $(a,b)\in\calD_{a,b}$
\begin{equation}
\forall x\in I,~
G_X(x) =
\systL
\frac{\ln(\frac{a+b^x}{a+1})}{\ln(\frac{a+b}{a+1})} & \text{if } a(1-b) >0 \
\frac{1-b^x}{1-b} & \text{if } a=+\infty \text{ and } b<1 \
x & \text{if } a=0 \text{ or } b=1. \
\systR
\label{expcurve:mbbefd:ab}
\end{equation}
The two special cases of $a(1-b)=0$
correspond to $\calD_{a,1}$ and $\calD_{0,b}$.
Note that the denominator is a normalizing constant to ensure the term belongs to $[0,1]$.
Distribution, density and quantile functions
Differentiating $G_X$, we obtain the following distribution function
for $(a,b)\in\calD_{a,b}$
\begin{equation}
\forall x\in I,~
F_X(x) =
\systL
\left(1-\frac{(a+1)b^x}{a+b^x}\right)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{if } a(1-b) >0 \
(1-b^x)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{if } a=+\infty \text{ and } b<1 \
\ind_{[1,+\infty[}(x) & \text{if } a=0 \text{ or } b=1. \
\systR
\label{cdf:mbbefd:ab}
\end{equation}
Note that the MBBEFD distribution is a mixed-type distribution with mass probability at $x=1$
\begin{equation}
P(X=1) = \frac{(a+1)b}{a+b} = p_{a,b},
\label{mass0:mbbefd:ab}
\end{equation}
which equals to 1 when $a(1-b)=0$. In other words, for $\calD_{a,1}$ and $\calD_{0,b}$, $X$ has a Dirac distribution at $x=1$.
When $a=+\infty$, the total loss probability is $P(X=1)=b$.
For $a(1-b)>0$, the improper density function is
\begin{equation}
\tilde f_X(x)
\systL
-\frac{a(a+1)b^x\ln(b)}{(a+b^x)^2}\ind_{[0,1[}(x) & \text{if } a(1-b) >0 \
-\ln(b)b^x\ind_{[0,1[}(x) & \text{if } a=+\infty \text{ and } b<1 \
0 & \text{if } a=0 \text{ or } b=1. \
\systR
\label{impd:mbbefd:ab}
\end{equation}
The quantile function is
\begin{equation}
\forall p\in [0,1],~
q_X(p) =
\systL
\frac{\ln\left(\frac{(1-p)a}{a+p}\right)}{\ln(b)} \ind_{[0,1-p_{a,b}[}(p) + \ind_{[1-p_{a,b},1]}(p) & \text{if } a(1-b) >0 \
\frac{\ln(1-p)}{\ln(b)}\ind_{[0,1-b[}(p) + \ind_{[1-b,1]}(p) & \text{if } a=+\infty \text{ and } b<1 \
\ind_{]0,1]}(p) & \text{if } a=0 \text{ or } b=1. \
\systR
\label{q:mbbefd:ab}
\end{equation}
Moments
Using the definition of the exposure curve, we have
$E(X) = 1/G_X'(0).$
The expectation for $\mbox{MBBEFD}(a,b)$ is
$$
E(X)=\frac{\ln(\frac{a+b}{a+1})}{\ln(b)} (a+1).
$$
When $a=0$ or $b=1$, the expectation is simply $E(X)=1$.
The MBBEFD distribution, second parametrization
\label{sec:mbbefd2p}
Parameter domain
For fitting purposes and for verifying parameter constraints,
[@bernegger97] proposed a second parametrization $\mbox{MBBEFD}(g,b)$.
Using the following parameter $g=1/p_{a,b}$, it is possible to reformulate the $\mbox{MBBEFD}(a,b)$.
That is
$$
g= \frac{a+b}{(a+1)b}
\Leftrightarrow
a=\frac{(g-1)b}{1-gb}.
$$
So $g\geq 1$ guarantees that $\frac{a+b}{(a+1)b}\in [0,1]$, in addition to $b>0$.
The special case $g=1$ leading to a Dirac distribution at $x=1$ corresponds to $a(1-b)=0$ in the previous parametrization.
The parameter domain is
$$
\widetilde{\calD}_{g,b}= \left{
(g,b)\in\R^2, b>0, g\geq 1
\right}.
$$
Characterization by the exposure curve
The exposure curve is defined for $x\in[0,1]$ as
\begin{equation}
G_X(x)=
\systL
\frac{\ln(\frac{(g-1)b}{1-b}+\frac{1-gb}{1-b}b^x)}{\ln(gb)} & \text{ if } g>1, b\neq 1, b\neq 1/g \
\frac{\ln(1+(g-1)x)}{\ln( g)} & \text{ if } g>1, b= 1 \
\frac{1-b^x}{1-b} & \text{ if } g>1, bg= 1 \
x & \text{ if } g=1 \text{ or } b= 0
\systR
\label{expcurve:mbbefd:gb}
\end{equation}
Note that the case $g>1, bg=1$ implies $g=1/b, b<1$.
Distribution, density and quantile functions
The resulting distribution function is
\begin{equation}
F_X(x)=
\systL
\left(1- \frac{1-b}{(g-1)b^{1-x}+1-gb}\right)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{ if } g>1, b\neq 1, b\neq 1/g \
\left(1- \frac{1}{1+(g-1)x}\right)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{ if } g>1, b= 1 \
(1-b^x)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{ if } g>1, bg= 1 \
\ind_{[1,+\infty[}(x) & \text{ if } g=1 \text{ or } b= 0
\systR
\label{cdf:mbbefd:gb}
\end{equation}
As in the previous parametrization, there is a non-null probability at $x=1$,
\begin{equation}
P(X=1) = 1/g.
\label{mass0:mbbefd:gb}
\end{equation}
The improper density function is for $x\in ]0,1[$
\begin{equation}
\tilde f_X(x) =
\systL
-\frac{(1-b)(g-1)\ln(b)b^{1-x}}{((g-1)b^{1-x}+1-gb)^2} & \text{ if } g>1, b\neq 1, b\neq 1/g \
\frac{g-1}{(1+(g-1)x)^2} & \text{ if } g>1, b= 1 \
-\ln(b) b^x & \text{ if } g>1, bg= 1 \
0 & \text{ if } g=1 \text{ or } b= 0
\systR
\label{impd:mbbefd:gb}
\end{equation}
The quantile function is
\begin{equation}
\forall p\in [0,1],~
q_X(p) =
\systL
\left(1-\frac{\ln\left(\frac{gb-1}{g-1} +\frac{1-b}{(1-p)(g-1)}\right)}{\ln(b)} \right)
\ind_{[0,1-1/g[}(p) + \ind_{[1-1/g,1]}(p) & \text{ if } g>1, b\neq 1, b\neq 1/g \
\frac{p}{(1-p)(g-1)}\ind_{[0,1-1/g[}(p) + \ind_{[1-1/g,1]}(p) & \text{ if } g>1, b= 1 \
\frac{\ln(1-p)}{\ln(b)}\ind_{[0,1-1/g[}(p) + \ind_{[1-1/g,1]}(p) & \text{ if } g>1, bg= 1 \
\ind_{]0,1]}(p) & \text{ if } g=1 \text{ or } b= 0
\systR
\label{q:mbbefd:gb}
\end{equation}
Moments
Let us compute the first two moments.
$$
E(X) = 1/G_X'(0) =
\systL
\frac{\ln(gb)(1-b)}{\ln(b)(1-gb)} & \text{ if } g>1, b\neq 1, b\neq 1/g \
\frac{\ln(g)}{g-1} & \text{ if } g>1, b= 1 \
\frac{b-1}{\ln(b)} & \text{ if } g>1, bg= 1 \
1 & \text{ if } g=1 \text{ or } b= 0
\systR
$$
Fitting methods and pricing examples
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 \pkg{mbbefd} package.
As its name suggest, the maximum likelihood function consists in maximizing
the likelihood function defined as
$$
\calL(\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
$\theta$ 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)\ind_{x=1}$ and $F_X'(x)\ind_{x\neq 1}$.
For example when considering the beta distribution, the likelihood function to be maximized is
$$
\calL((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
$$
\systL
P(X = 1;\theta) = tl_n \
E(X; \theta) = \bar x_n\
\vdots \
E(X^k; \theta) = m_{n,k}\
\systR
\text{ where }
\systL
tl_n = \frac{1}{n}\sum_{i=1}^n \ind_{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.
\systR
$$
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 \pkg{fitDR} function of the \pkg{mbbefd} package returns an object of class
\pkg{"DR"} inheriting from the \pkg{"fitdist"} class for which generic functions
are available \pkg{print, plot, summary, gofstat}.
There is also a bootstrap procedure implemented in the \pkg{bootDR} 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)
b1 <- bootDR(f1, niter=20)
summary(b1)
par(cex.main=0.8, cex.lab=0.8, cex.axis=0.8, cex=0.8)
denscomp(f1, demp=TRUE, main="Histogram and theoretical densities")
par(cex.main=0.8, cex.lab=0.8, cex.axis=0.8, mar=c(2,2,2,2))
plot(b1, enhance=TRUE, main="Bootstrapped value of parameters")
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 \pkg{fitdistrplus} package, we can compare multiple fits on the
same graph or with the same statistics. We just have to pass a list of
\pkg{"DR"} objects.
f2 <- fitDR(x, "oibeta", method="mle")
f3 <- fitDR(x, "oiunif", method="mle")
gofstat(list(f1, f2, f3))
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 \pkg{eccomp}. 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
$$
\mbox{CHF } 3,500,000 \mbox{ 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$:
$$
\mbox{CHF } 2,908,182
\mbox{ 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.
Client <- scan(textConnection(
'150 75 250 200 400 325 600 500 800 700 1000 900 1250 1125 1500 1375 1750 1625 2000 1875 2500 2250 3000 2750 4000 3500 5500 4750 9000 7250 12500 10750 18000 15250 24000 21000 36000 30000 48000 42000 72000 60000 90000 81000'))
MaxMPL <- Client[seq(from=1, to=length(Client),by = 2)]
MeanMPLGrossLoss <- Client[seq(from=2, to=length(Client),by = 2)]
GrossPremium <- scan(textConnection(
'33434 14568 6324 4584 3341 1405 1169 683 613 554 700 552 1194 1490 4177 3527 3249 2712 2588 1988 657 1918'
))
ExposureCurve <- c(rep(1.5,3), rep(2.0,3), rep(3.0, 4), rep(4.0, 12))
ClientData <- data.frame(MaxMPL, MeanMPLGrossLoss,
GrossPremium, ExposureCurve)
ClientData2 <- ClientData
names(ClientData2) <- c("Max MPL '000", "Mean MPL Gross Loss '000",
"Gross Premium '000", "Exposure Curve Parameter c")
library(pander)
panderOptions('big.mark', ',')
panderOptions('table.split.cells', 8)
pander(ClientData2, justify = rep('right', ncol(ClientData2)))
## Example layer
D <- 1.5e3*457/550
L <- 5e3*457/550
MPLoss <- 3.5e6
retainedDed <- D/MPLoss*1000
retainedLoss <- ecMBBEFD(retainedDed, b=swissRe(4)['b'], g=swissRe(4)['g'])
prem <- 1194 * (1 - retainedLoss)
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.
Losses below deductible
Losses with a maximum MPL below the deductible are not relevant and hence
no reinsurance premium needs to be calculated.
Losses above deductible, but below limit
As an example we use the band that has a maximum MPL of CHF 4m,
with a mean MPL gross loss of CH 3.5m and gross premium of CH 1.194m.
The client would retain $1,246/3,500=35.6\%$ of the
maximum MPL. Using our exposure curve we can read of that this would reduce our
expected claims burden by $35.6\%$.
Therefore the reinsurance premium for this band is 1,194k
$\cdot 20.5\% =244.8k$.
Losses above limit
Here we select the band with a maximum MPL of CHF 90m and a mean gross
MPL of CHF 81m. Because of the limit of CHF $4154.5$ only that
amount matters. From the exposure curve we can read
of again the average proportion of loss paid by the cedant
$0.8\%$.
To calculate the premium we first have to derive the proportion of
premium that is attached to the risk below the limit: 4,155k/8,100k
$\cdot 1,918k=98.4k$.
Multiplying the premium with the reinsurer's share of the loss gives
us a premium of CHF
$23.4k$.
Premium Rate
For the overall portfolio we can calculate the premium rate, using the expected loss ratio of 55%:
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)
Applying the function to our data we get a premium rate of $1.47\%$.
Differences to the Swiss Re paper result from rounding errors.
References
spedygiorgio/mbbefd documentation built on Sept. 2, 2023, 1:55 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library("mbbefd") #library("knitcitations") #cleanbib() #options("citation_format" = "pandoc") #bib <- read.bibtex("mbbefd.bib") library(lattice) my.settings <- canonical.theme(color=FALSE) my.settings[['fontsize']] = list(text = 8, points = 4) my.settings[['strip.background']]$col <- "darkgrey" my.settings[['strip.border']]$col<- "black"
Introduction and general notation
Exposure rating is a tool for insurance pricing that allocates premium to bands of damage ratios or severity of losses. First ideas were published in [@salzmann1963]. It is often used to price non-proportional reinsurance contracts, such as excess of loss (XL) reinsurance.
Exposure rating uses the loss experience of a similar portfolio of policies to estimate the expected losses of the portfolio to be covered. The method is frequently used as a benchmark when there is no sufficient credible claims history from the client.
In this vignette, we first present the general notation and concepts of loss modelling. Secondly, we focus on destruction models implemented in the \pkg{mbbefd} package. One popular model uses the MBBEFD distribution introduced by [@bernegger97]. Finally, we provide an example of pricing of XL contract.
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 $\mu$ and $\sigma$ can be derived $$ \sigma^2=\log(1+CV^2),~~ \mu =\log(M) - \sigma^2/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)$.
n <- 100 Loss <- seq(1, 150, length=n) mu <- 4.13 sigma <- 0.29 CDF <- plnorm(Loss, mu, sigma) Density <- dlnorm(Loss, mu, sigma) Survival <- 1 - CDF dat <- data.frame(Loss=rep(Loss, 3), Value=c(Density, CDF, Survival), Type=gl(3, n, labels=c("PDF: f(x)", "CDF: F(x)", "SF: S(x)=1-F(x)"), ordered=TRUE)) xyplot(Value ~ Loss | Type, data=dat, ylab="", layout=c(3,1), as.table=TRUE, type="l", par.settings = my.settings, par.strip.text=list(col="white", font=2), scales=list(relation="free", alternating=1))
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] = \int_{0}^\infty x\, f(x) dx. $$
Limited Expected Loss
Let's further assume that losses are limited by $\alpha$ in the contract. Then the limited expected value (LEV), written as $E[X \wedge \alpha]$, is given as $$ E[X \wedge \alpha] = \int_0^\alpha x \, f(x) dx + \int_\alpha^\infty \alpha \, f(x) dx. $$
To evaluate this sum of integrals we recall that $F(x)=\int f(x) dx$ and $F(\infty)=1$. Hence $$ \int_\alpha^\infty \alpha \, f(x) dx = \alpha \int_\alpha^\infty f(x) dx = \alpha - \alpha F(\alpha). $$
The integration by parts theorem helps with the first part of the integral $$ \int_0^\alpha x \, f(x) dx = \alpha F(\alpha) - \int_0^\alpha F(x) dx. $$
Therefore, we retrieve the well-known identity $$ E[X \wedge \alpha] = \alpha - \int_0^\alpha F(x) dx = \int_0^\alpha 1-F(x) dx = \int_0^\alpha S(x) dx. $$
The integral describes the area above the CDF up to $\alpha$, or the area under the survival function up to $\alpha$. In the example below the limit $\alpha$ was set to 100.
alpha <- 100 fxy <- function(x,y,...){ x1 <- as.numeric(x[x<=alpha]) if(panel.number()<2){ panel.polygon(c(x1, rev(x1)), rev(c(rep(1, length(x1)), rev(y[x<=alpha]))), col="skyblue", border=NA) }else{ panel.polygon(c(x1, rev(x1)), c(rep(0, length(x1)), rev(y[x<=alpha])), col="skyblue", border=NA) } panel.xyplot(x,y,...) #panel.text(x=20, y=0.4, label=paste("LEV[X]", cex=2) } xyplot(Value ~ Loss | Type, data=subset(dat, !Type %in% "PDF: f(x)"), ylab="", panel=fxy, layout=c(2,1), as.table=TRUE, type="l", par.settings = my.settings, par.strip.text=list(col="white", font=2), scales=list(relation="free", alternating=1))
Loss Cost of a Layer
Suppose we want to insure the layer of claims between 80 and 100.
alpha1 <- 80 alpha2 <- 100 xyplot(Value ~ Loss | Type, data=subset(dat, !Type %in% "PDF: f(x)"), ylab="", panel=function(x,y,...){ x1 <- c(x[x >= alpha1 & x<=alpha2]) if(panel.number()<2){ panel.polygon( c(x1, rev(x1)), rev(c(rep(1, length(x[x >= alpha1 & x<=alpha2])), rev(y[x >= alpha1 & x<=alpha2]))), col="skyblue", border=NA) }else{ panel.polygon( c(x1, rev(x1)), c(rep(0, length(x[x >= alpha1 & x<=alpha2])), rev(y[x >= alpha1 & x<=alpha2])), col="skyblue", border=NA) } panel.xyplot(x,y,...) }, layout=c(2,1), as.table=TRUE, t="l", par.settings = my.settings, par.strip.text=list(col="white", font=2), scales=list(relation="free", alternating=1))
The expected loss cost would be the difference between the limited expected values of $E[X \wedge 100] - E[X \wedge 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)
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)
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.
Normalising loss experience using TIV and MPL
To make risks more comparable we look at ratio of losses to the underlying exposure, where the exposure is given as the sum insured (SI), or better the total insured value (TIV), or perhaps as the maximum probable loss (MPL).
Other definitions and measures are also popular, such as possible maximum loss (PML), estimate maximum loss (EML), or maximum foreseeable loss (MFL).
While the TIV and SI are straightforward to understand, the other metrics can be little more challenging to assess.
Suppose we insure a big production plant with of several buildings against fire. It is unlikely that a fire will destroy the whole facility, instead it is believed that the distance between the building will ensure that fires can be contained in a local area. Hence the MPL might be only the highest value of any of those buildings.
Shifting the metric from loss amounts to damage ratio or destruction rates (loss as a % of the exposure metric) allows us to compare and benchmarks risks.
Analysing deductibles
Most insurance policies are written with a deductible, so that the insured will cover the first $X of the losses.
Thus, the reinsurer needs to understand, by how much the claims burden is reduced for a given deductible.
From the previous section we have learned how to calculate the limited expected value, which is the loss cost to the insured.
The exposure curve (or also called first loss curve or original loss cure ) is defined as the proportion of the LEV for a given deductible $d$ compared to the overall expected claims cost $E[X]$ $$ G(d) = \frac{E[X \wedge d]}{E[X]}. $$
For our example, using a log-normal distribution and a MPL of 200 we get:
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)
The steepness of the curve is related to the severity of the loss distribution. The closer the curve is to the diagonal the greater the proportion of large loss.
If all losses were total losses then the exposure curve would be identical with the diagonal.
Different perils and exposures will have different exposure curves. Well known exposure curves are those used by Swiss Re and Lloyd's.
Destruction rate models
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\in]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)\ind_{[0,1[}(x) + \ind_{[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. [@antal03].
Uniform distribution
\label{sec:unif}
The most trivial example of exposure curve is obtained for the uniform distribution on $I$. We consider $F_X(x)=x$ leading to $$ G_X(d) = d(2-d). $$
Beta distribution
\label{sec:beta}
A more interesting example is obtained for the Beta distribution on $I$ (e.g.[@kotzjohnsonbalak94v2] for which the density is $f_X(x) = x^{a-1}(1-x)^{b-1}/\beta(a,b)$ for $x\in ]0,1[$ and $a,b>0$ where $\beta(.,.)$ denotes the beta function, see e.g.[@nist10]. The distribution function is obtained in terms of the incomplete beta ratio function $F_X(x) = \beta(x;a,b)/\beta(a,b) = I(x;a,b)$. We get $$ G_X(d) \beta(a,b;x) - \frac{x^a(1-x)^b}{a\beta(a,b)} + x(1-\beta(a,b;x)) \frac{a+b}{a}, $$ where $\beta(.;.,.)$ denotes the incomplete beta function.
One-inflated distributions
\label{sec:oidistr:generic}
Characterizations
\label{sec:oidistr:generic:charac} Let us consider a continuous distribution function $F_0$ of a random variable $X_0$. The corresponding distribution function of the one-inflated random variable $X_1$ is $$ F_1(x) = (1-p_1) F_0(x) + p_1 \ind_{[1,+\infty[}(x). $$ There is no density but an improper density $(1-p)F_0'(x)$ and a probability mass $p_1$ at $x=1$.
The one-inflated beta distribution
We consider the one-inflated beta distribution.
Using Section \ref{sec:beta}, we obtain the following distribution function
$$
F_X(x) =
\systL
0 & \text{if } x<0 \
I(x;a,b)(1-p_1) & \text{if } 0\leq x <1 \
1 & \text{if } x\geq 1
\systR
$$
where $I(x;a,b)$ denotes the incomplete beta ratio function.
This leads to a non-null probability at $x=1$, $P(X=1)=p_1$.
The improper density function is
$$
\tilde f_X(x) = (1-p_1) \frac{x^{a-1}(1-x)^{b-1}}{\beta(a,b)}.
$$
The expectation is
$$
E(X)
= p_1 + (1-p_1)\frac{a}{a+b}.
$$
The exposure curve is
$$
G_X(d) =
\frac{(1-p_1) \left(1 - I(d;a,b) \frac{b}{a+b} - \frac{d^a(1-d)^b}{(a+b)\beta(a,b)}\right) + p_1d
}{p_1 + (1-p_1)\frac{a}{a+b}}.
$$
The MBBEFD distribution, first parametrization
We denote this first parametrization by $\mbox{MBBEFD}(a,b)$. We define the parameter domain $\calD_{a,b}$ as \begin{equation} \calD_{a,b} = {(a,b)\in\R^2, a+1>0, a(1-b)>0, b>0} \cup {(a,b), a=+\infty, b<1}. \label{eq:paramdomain:mbbefd:ab} \end{equation} Let us note that this domain includes two particular sets $\calD_{a,1}={(a,1), a+1>0}$ and $\calD_{0,b}={(0,b), b>0}$.
Characterization by the exposure curve
The MBBEFD distribution is defined by the following exposure curve for $(a,b)\in\calD_{a,b}$ \begin{equation} \forall x\in I,~ G_X(x) = \systL \frac{\ln(\frac{a+b^x}{a+1})}{\ln(\frac{a+b}{a+1})} & \text{if } a(1-b) >0 \ \frac{1-b^x}{1-b} & \text{if } a=+\infty \text{ and } b<1 \ x & \text{if } a=0 \text{ or } b=1. \ \systR \label{expcurve:mbbefd:ab} \end{equation} The two special cases of $a(1-b)=0$ correspond to $\calD_{a,1}$ and $\calD_{0,b}$. Note that the denominator is a normalizing constant to ensure the term belongs to $[0,1]$.
Distribution, density and quantile functions
Differentiating $G_X$, we obtain the following distribution function for $(a,b)\in\calD_{a,b}$ \begin{equation} \forall x\in I,~ F_X(x) = \systL \left(1-\frac{(a+1)b^x}{a+b^x}\right)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{if } a(1-b) >0 \ (1-b^x)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{if } a=+\infty \text{ and } b<1 \ \ind_{[1,+\infty[}(x) & \text{if } a=0 \text{ or } b=1. \ \systR \label{cdf:mbbefd:ab} \end{equation} Note that the MBBEFD distribution is a mixed-type distribution with mass probability at $x=1$ \begin{equation} P(X=1) = \frac{(a+1)b}{a+b} = p_{a,b}, \label{mass0:mbbefd:ab} \end{equation} which equals to 1 when $a(1-b)=0$. In other words, for $\calD_{a,1}$ and $\calD_{0,b}$, $X$ has a Dirac distribution at $x=1$. When $a=+\infty$, the total loss probability is $P(X=1)=b$.
For $a(1-b)>0$, the improper density function is \begin{equation} \tilde f_X(x) \systL -\frac{a(a+1)b^x\ln(b)}{(a+b^x)^2}\ind_{[0,1[}(x) & \text{if } a(1-b) >0 \ -\ln(b)b^x\ind_{[0,1[}(x) & \text{if } a=+\infty \text{ and } b<1 \ 0 & \text{if } a=0 \text{ or } b=1. \ \systR \label{impd:mbbefd:ab} \end{equation} The quantile function is \begin{equation} \forall p\in [0,1],~ q_X(p) = \systL \frac{\ln\left(\frac{(1-p)a}{a+p}\right)}{\ln(b)} \ind_{[0,1-p_{a,b}[}(p) + \ind_{[1-p_{a,b},1]}(p) & \text{if } a(1-b) >0 \ \frac{\ln(1-p)}{\ln(b)}\ind_{[0,1-b[}(p) + \ind_{[1-b,1]}(p) & \text{if } a=+\infty \text{ and } b<1 \ \ind_{]0,1]}(p) & \text{if } a=0 \text{ or } b=1. \ \systR \label{q:mbbefd:ab} \end{equation}
Moments
Using the definition of the exposure curve, we have $E(X) = 1/G_X'(0).$ The expectation for $\mbox{MBBEFD}(a,b)$ is $$ E(X)=\frac{\ln(\frac{a+b}{a+1})}{\ln(b)} (a+1). $$ When $a=0$ or $b=1$, the expectation is simply $E(X)=1$.
The MBBEFD distribution, second parametrization
\label{sec:mbbefd2p}
Parameter domain
For fitting purposes and for verifying parameter constraints, [@bernegger97] proposed a second parametrization $\mbox{MBBEFD}(g,b)$. Using the following parameter $g=1/p_{a,b}$, it is possible to reformulate the $\mbox{MBBEFD}(a,b)$. That is $$ g= \frac{a+b}{(a+1)b} \Leftrightarrow a=\frac{(g-1)b}{1-gb}. $$ So $g\geq 1$ guarantees that $\frac{a+b}{(a+1)b}\in [0,1]$, in addition to $b>0$. The special case $g=1$ leading to a Dirac distribution at $x=1$ corresponds to $a(1-b)=0$ in the previous parametrization. The parameter domain is $$ \widetilde{\calD}_{g,b}= \left{ (g,b)\in\R^2, b>0, g\geq 1 \right}. $$
Characterization by the exposure curve
The exposure curve is defined for $x\in[0,1]$ as \begin{equation} G_X(x)= \systL \frac{\ln(\frac{(g-1)b}{1-b}+\frac{1-gb}{1-b}b^x)}{\ln(gb)} & \text{ if } g>1, b\neq 1, b\neq 1/g \ \frac{\ln(1+(g-1)x)}{\ln( g)} & \text{ if } g>1, b= 1 \ \frac{1-b^x}{1-b} & \text{ if } g>1, bg= 1 \ x & \text{ if } g=1 \text{ or } b= 0 \systR \label{expcurve:mbbefd:gb} \end{equation} Note that the case $g>1, bg=1$ implies $g=1/b, b<1$.
Distribution, density and quantile functions
The resulting distribution function is \begin{equation} F_X(x)= \systL \left(1- \frac{1-b}{(g-1)b^{1-x}+1-gb}\right)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{ if } g>1, b\neq 1, b\neq 1/g \ \left(1- \frac{1}{1+(g-1)x}\right)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{ if } g>1, b= 1 \ (1-b^x)\ind_{[0,1[}(x) + \ind_{[1,+\infty[}(x) & \text{ if } g>1, bg= 1 \ \ind_{[1,+\infty[}(x) & \text{ if } g=1 \text{ or } b= 0 \systR \label{cdf:mbbefd:gb} \end{equation} As in the previous parametrization, there is a non-null probability at $x=1$, \begin{equation} P(X=1) = 1/g. \label{mass0:mbbefd:gb} \end{equation} The improper density function is for $x\in ]0,1[$ \begin{equation} \tilde f_X(x) = \systL -\frac{(1-b)(g-1)\ln(b)b^{1-x}}{((g-1)b^{1-x}+1-gb)^2} & \text{ if } g>1, b\neq 1, b\neq 1/g \ \frac{g-1}{(1+(g-1)x)^2} & \text{ if } g>1, b= 1 \ -\ln(b) b^x & \text{ if } g>1, bg= 1 \ 0 & \text{ if } g=1 \text{ or } b= 0 \systR \label{impd:mbbefd:gb} \end{equation} The quantile function is \begin{equation} \forall p\in [0,1],~ q_X(p) = \systL \left(1-\frac{\ln\left(\frac{gb-1}{g-1} +\frac{1-b}{(1-p)(g-1)}\right)}{\ln(b)} \right) \ind_{[0,1-1/g[}(p) + \ind_{[1-1/g,1]}(p) & \text{ if } g>1, b\neq 1, b\neq 1/g \ \frac{p}{(1-p)(g-1)}\ind_{[0,1-1/g[}(p) + \ind_{[1-1/g,1]}(p) & \text{ if } g>1, b= 1 \ \frac{\ln(1-p)}{\ln(b)}\ind_{[0,1-1/g[}(p) + \ind_{[1-1/g,1]}(p) & \text{ if } g>1, bg= 1 \ \ind_{]0,1]}(p) & \text{ if } g=1 \text{ or } b= 0 \systR \label{q:mbbefd:gb} \end{equation}
Moments
Let us compute the first two moments. $$ E(X) = 1/G_X'(0) = \systL \frac{\ln(gb)(1-b)}{\ln(b)(1-gb)} & \text{ if } g>1, b\neq 1, b\neq 1/g \ \frac{\ln(g)}{g-1} & \text{ if } g>1, b= 1 \ \frac{b-1}{\ln(b)} & \text{ if } g>1, bg= 1 \ 1 & \text{ if } g=1 \text{ or } b= 0 \systR $$
Fitting methods and pricing examples
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 \pkg{mbbefd} package.
As its name suggest, the maximum likelihood function consists in maximizing the likelihood function defined as $$ \calL(\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 $\theta$ 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)\ind_{x=1}$ and $F_X'(x)\ind_{x\neq 1}$.
For example when considering the beta distribution, the likelihood function to be maximized is $$ \calL((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 $$ \systL P(X = 1;\theta) = tl_n \ E(X; \theta) = \bar x_n\ \vdots \ E(X^k; \theta) = m_{n,k}\ \systR \text{ where } \systL tl_n = \frac{1}{n}\sum_{i=1}^n \ind_{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. \systR $$ 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 \pkg{fitDR} function of the \pkg{mbbefd} package returns an object of class \pkg{"DR"} inheriting from the \pkg{"fitdist"} class for which generic functions are available \pkg{print, plot, summary, gofstat}. There is also a bootstrap procedure implemented in the \pkg{bootDR} 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) b1 <- bootDR(f1, niter=20) summary(b1)
par(cex.main=0.8, cex.lab=0.8, cex.axis=0.8, cex=0.8) denscomp(f1, demp=TRUE, main="Histogram and theoretical densities")
par(cex.main=0.8, cex.lab=0.8, cex.axis=0.8, mar=c(2,2,2,2)) plot(b1, enhance=TRUE, main="Bootstrapped value of parameters")
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 \pkg{fitdistrplus} package, we can compare multiple fits on the same graph or with the same statistics. We just have to pass a list of \pkg{"DR"} objects.
f2 <- fitDR(x, "oibeta", method="mle") f3 <- fitDR(x, "oiunif", method="mle") gofstat(list(f1, f2, f3)) 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 \pkg{eccomp}. 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 $$ \mbox{CHF } 3,500,000 \mbox{ 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$:
$$ \mbox{CHF } 2,908,182 \mbox{ 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.
Client <- scan(textConnection( '150 75 250 200 400 325 600 500 800 700 1000 900 1250 1125 1500 1375 1750 1625 2000 1875 2500 2250 3000 2750 4000 3500 5500 4750 9000 7250 12500 10750 18000 15250 24000 21000 36000 30000 48000 42000 72000 60000 90000 81000')) MaxMPL <- Client[seq(from=1, to=length(Client),by = 2)] MeanMPLGrossLoss <- Client[seq(from=2, to=length(Client),by = 2)] GrossPremium <- scan(textConnection( '33434 14568 6324 4584 3341 1405 1169 683 613 554 700 552 1194 1490 4177 3527 3249 2712 2588 1988 657 1918' )) ExposureCurve <- c(rep(1.5,3), rep(2.0,3), rep(3.0, 4), rep(4.0, 12)) ClientData <- data.frame(MaxMPL, MeanMPLGrossLoss, GrossPremium, ExposureCurve) ClientData2 <- ClientData names(ClientData2) <- c("Max MPL '000", "Mean MPL Gross Loss '000", "Gross Premium '000", "Exposure Curve Parameter c") library(pander) panderOptions('big.mark', ',') panderOptions('table.split.cells', 8) pander(ClientData2, justify = rep('right', ncol(ClientData2))) ## Example layer D <- 1.5e3*457/550 L <- 5e3*457/550 MPLoss <- 3.5e6 retainedDed <- D/MPLoss*1000 retainedLoss <- ecMBBEFD(retainedDed, b=swissRe(4)['b'], g=swissRe(4)['g']) prem <- 1194 * (1 - retainedLoss)
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.
Losses below deductible
Losses with a maximum MPL below the deductible are not relevant and hence no reinsurance premium needs to be calculated.
Losses above deductible, but below limit
As an example we use the band that has a maximum MPL of CHF 4m, with a mean MPL gross loss of CH 3.5m and gross premium of CH 1.194m.
The client would retain $1,246/3,500=35.6\%$ of the maximum MPL. Using our exposure curve we can read of that this would reduce our expected claims burden by $35.6\%$.
Therefore the reinsurance premium for this band is 1,194k $\cdot 20.5\% =244.8k$.
Losses above limit
Here we select the band with a maximum MPL of CHF 90m and a mean gross MPL of CHF 81m. Because of the limit of CHF $4154.5$ only that amount matters. From the exposure curve we can read of again the average proportion of loss paid by the cedant $0.8\%$.
To calculate the premium we first have to derive the proportion of premium that is attached to the risk below the limit: 4,155k/8,100k $\cdot 1,918k=98.4k$.
Multiplying the premium with the reinsurer's share of the loss gives us a premium of CHF $23.4k$.
Premium Rate
For the overall portfolio we can calculate the premium rate, using the expected loss ratio of 55%:
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)
Applying the function to our data we get a premium rate of $1.47\%$. Differences to the Swiss Re paper result from rounding errors.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.