---
title: "Exposure rating, destruction rate models and the mbbefd package"
author: "Christophe Dutang, Markus Gesmann, Giorgio Spedicato"
date: "`r Sys.Date()`"
output: 
  pdf_document:
    number_sections: true
bibliography: mbbefd.bib
vignette: |
  %\VignetteIndexEntry{Exposure rating, destruction rate models and the mbbefd package} 
  %\VignetteEngine{knitr::rmarkdown} 
  \usepackage[utf8]{inputenc}
header-includes:
- \newcommand{\ind}{1\!\!1}
- \newcommand{\R}{\mathbb R}
- \newcommand{\N}{\mathbb N}
- \newcommand{\calD}{\mathcal D}
- \newcommand{\calL}{\mathcal L}
- \newcommand{\systL}{\left\{\begin{array}{ll}}
- \newcommand{\systR}{\end{array}\right.}
- \newcommand{\matL}{\left(\begin{matrix}}
- \newcommand{\matR}{\end{matrix}\right)}
- \newcommand{\detL}{\left|\begin{matrix}}
- \newcommand{\detR}{\end{matrix}\right|}
- \newcommand{\pkg}{\textsf}

---


```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
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"  
```


<!------------------ SECTION 1 ------------------> 
# 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 `r citep(bib['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 `r citep(bib['bernegger97'])`.
Finally, we provide an example of pricing of XL contract.



## Loss Distribution

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

To keep it simple I 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 via:

$$
\begin{aligned}
\sigma^2&=\log(1+CV^2) \\
\mu &=\log(M) - \sigma^2/2 
\end{aligned}
$$

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

```{r, fig.height=2.5, fig.width=6, echo=FALSE, warning=FALSE}
library(lattice)
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, t="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) 
is then:

$$
E[X] = \int_{-\infty}^\infty x\, f(x) dx
$$

Assuming losses are not negative ($P(X<0)=0$) this simplifies to:

$$ 
E[X] = \int_{0}^\infty x\, f(x) dx
$$

## Limited Expected Loss

Let's further assume that losses are limited by $\alpha$ in my 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 I 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 $\int u(x) v'(x) \, dx = u(x) v(x) - \int v(x) \, u'(x)  dx$ helps me with the first part of the integral:

$$
\int_0^\alpha x \, f(x) dx = \alpha F(\alpha) - \int_0^\alpha F(x) dx
$$

Therefore:

$$
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.


```{r LEV, fig.height=2.5, fig.width=5, echo=FALSE}
alpha <- 100
xyplot(Value ~ Loss | Type, 
       data=subset(dat, !Type %in% "PDF: f(x)"), 
       ylab="", 
       panel=function(x,y,...){
         x1 <- c(x[x<=alpha])
         if(panel.number()<2){
           panel.polygon(c(x1, rev(x1)),
                          rev(c(rep(1, length(x[x<=alpha])),
                           rev(y[x<=alpha]))),
                         col="skyblue", border=NA)
         }else{
           panel.polygon(c(x1, rev(x1)),
                         c(rep(0, length(x[x<=alpha])),
                           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)
       },
       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))
```

## Loss Cost of a Layer 

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


```{r LEV2, fig.height=2.5, fig.width=5, echo=FALSE}
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]$.

```{r}
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):

```{r}
(ILF <- integrate(S, 0, 100)$value / integrate(S, 0, 80)$value)
```

Therefore we would expect the LEV to increase by `r round((ILF-1)*100,1)`% 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:

```{r ExposureCurve, dev.args=list(pointsize=8), fig.height=3, fig.width=3}
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.

<!------------------ SECTION 2 ------------------> 
# 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.`r citep(bib['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.`r citep(bib['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.`r citep(bib['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, 
`r citep(bib['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)\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
$$



<!------------------ SECTION 3 ------------------> 
# 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 \pkg{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 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).
```{r simufit, echo=TRUE, fig.height=5, fig.width=5, warning=FALSE}
library(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(mfrow=c(1, 2))
denscomp(f1, demp=TRUE)
plot(b1, enhance=TRUE)
```
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.

```{r simufit2, echo=TRUE, fig.height=5, fig.width=10, warning=FALSE}
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
```{r simufit3, echo=TRUE, fig.height=5, fig.width=5, warning=FALSE}
eccomp(list(f1, f2, f3), leg=c("mbbefd", "oibeta", "oiunif"), do.points=FALSE)
```



## Computation of rate on line (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=`r round(457/550,2)`$:
$$
\mbox{CHF } `r format(3500000*457/550, big.mark=",")` 
\mbox{ xs CHF } `r format(1500000*457/550, big.mark=",")`
$$
 
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.

```{r SwissReExample, message=FALSE, echo=FALSE, warning=FALSE}
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 have to consider three scenarios:

### 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 I 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=`r round(retainedDed*100,1)`% of the 
maximum MPL. Using our exposure curve we can read of that this would reduce our 
expected claims burden by `r round(retainedLoss*100,1)`%.

Therefore the reinsurance premium for this band is 1,194k $\cdot$ `r round((100-retainedLoss*100),1)`% =`r round(prem,1)`k. 

### Losses above limit

Here I select the band with a maximum MPL of CHF 90m and a mean gross
MPL of  CHF 81m. Because of the limit of CHF `r round(L, 1)` only that amount matters. From the exposure curve we can read
of again the average proportion of loss paid by the cedant
`r round(ecMBBEFD(D/L, b=swissRe(4)['b'], g=swissRe(4)['g']),1)`%. 

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=`r round(4155/81000*1918,1)`k.

Multiplying the premium with the reinsurer's share of the loss gives
us a premium of CHF `r round(4155/81000*1918*(1-ecMBBEFD(D/L, b=swissRe(4)['b'], g=swissRe(4)['g'])),1)`k.

### Rate on Line

For the overall portfolio we can calculate the rate on line, which is the loss cost
of the layer compared to the total loss cost, using the expected loss ratio of 55%:

```{r SwissRe2}
XL_ROL <- function(Deductible, Limit, MaxMPL, 
                   MeanMPL, GrossPremium, C, ULR){
  DedPerMPL <- Deductible / ifelse(MaxMPL < Deductible, Deductible, 
                                   ifelse(MaxMPL<Limit, MaxMPL, 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
}

rol <- XL_ROL(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 rate on line of `r round(rol*100,2)`%.
Differences to the Swiss Re paper are a result of rounding errors.


# References
