Maximum likelihood estimates for multivariate distributions

September 22, 2012
By

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Consider our loss-ALAE dataset, and – as in Frees & Valdez (1998) – let us fit a parametric model, in order to price a reinsurance treaty. The dataset is the following,

```> library(evd)
> data(lossalae)
> Z=lossalae
> X=Z[,1];Y=Z[,2]```

The first step can be to estimate marginal distributions, independently. Here, we consider lognormal distributions for both components,

```> Fempx=function(x) mean(X<=x)
> Fx=Vectorize(Fempx)
> u=exp(seq(2,15,by=.05))
> plot(u,Fx(u),log="x",type="l",
+ xlab="loss (log scale)")
> Lx=function(px) -sum(log(Vectorize(dlnorm)(
+ X,px[1],px[2])))
> opx=optim(c(1,5),fn=Lx)
> opx\$par
[1] 9.373679 1.637499
> lines(u,Vectorize(plnorm)(u,opx\$par[1],
+ opx\$par[2]),col="red")```

The fit here is quite good,

For the second component, we do the same,

```> Fempy=function(x) mean(Y<=x)
> Fy=Vectorize(Fempy)
> u=exp(seq(2,15,by=.05))
> plot(u,Fy(u),log="x",type="l",
+ xlab="ALAE (log scale)")
> Ly=function(px) -sum(log(Vectorize(dlnorm)(
+ Y,px[1],px[2])))
> opy=optim(c(1.5,10),fn=Ly)
> opy\$par
[1] 8.522452 1.429645
> lines(u,Vectorize(plnorm)(u,opy\$par[1],
+ opy\$par[2]),col="blue")```

It is not as good as the fit obtained on losses, but it is not that bad,

Now, consider a multivariate model, with Gumbel copula. We’ve seen before that it worked well. But this time, consider the maximum likelihood estimator globally.

```> Cop=function(u,v,a) exp(-((-log(u))^a+
+ (-log(v))^a)^(1/a))
> phi=function(t,a) (-log(t))^a
> cop=function(u,v,a) Cop(u,v,a)*(phi(u,a)+
+ phi(v,a))^(1/a-2)*(
+ a-1+(phi(u,a)+phi(v,a))^(1/a))*(phi(u,a-1)*
+ phi(v,a-1))/(u*v)
> L=function(p) {-sum(log(Vectorize(dlnorm)(
+ X,p[1],p[2])))-
+ sum(log(Vectorize(dlnorm)(Y,p[3],p[4])))-
+ sum(log(Vectorize(cop)(plnorm(X,p[1],p[2]),
+ plnorm(Y,p[3],p[4]),p[5])))}
> opz=optim(c(1.5,10,1.5,10,1.5),fn=L)
> opz\$par
[1] 9.377219 1.671410 8.524221 1.428552 1.468238```

Marginal parameters are (slightly) different from the one obtained independently,

```> c(opx\$par,opy\$par)
[1] 9.373679 1.637499 8.522452 1.429645
> opz\$par[1:4]
[1] 9.377219 1.671410 8.524221 1.428552```

And the parameter of Gumbel copula is close to the one obtained with heuristic methods in class.

Now that we have a model, let us play with it, to price a reinsurance treaty. But first, let us see how to generate Gumbel copula… One idea can be to use the frailty approach, based on a stable frailty. And we can use Chambers et al (1976) to generate a stable distribution. So here is the algorithm to generate samples from Gumbel copula

```> alpha=opz\$par[5]
> invphi=function(t,a) exp(-t^(1/a))
> n=500
> x=matrix(rexp(2*n),n,2)
> angle=runif(n,0,pi)
> E=rexp(n)
> beta=1/alpha
> stable=sin((1-beta)*angle)^((1-beta)/beta)*
+ (sin(beta*angle))/(sin(angle))^(1/beta)/
+ (E^(alpha-1))
> U=invphi(x/stable,alpha)
> plot(U)```

Here, we consider only 500 simulations,

Based on that copula simulation, we can then use marginal transformations to generate a pair, losses and allocated expenses,

```> Xloss=qlnorm(U[,1],opz\$par[1],opz\$par[2])
> Xalae=qlnorm(U[,2],opz\$par[3],opz\$par[4])```

In standard reinsurance treaties – see e.g. Clarke (1996) – allocated expenses are splited prorata capita between the insurance company, and the reinsurer. If $X$ denotes losses, and $Y$ the allocated expenses, a standard excess treaty can be has payoff

$\begin{cases} 0\cdot \boldsymbol{1}(X\leq R)\\ \displaystyle{\left(X-R+ \frac{X-R}{X}\cdot Y\right) \right )}\cdot\boldsymbol{1}(R

where $L$ denotes the (upper) limit, and $R$ the insurer’s retention. Using monte carlo simulation, it is then possible to estimate the pure premium of such a reinsurance treaty.

```> L=100000
> R=50000
> Z=((Xloss-R)+(Xloss-R)/Xloss*Xalae)*
+ (R<=Xloss)*(Xloss<L)+
+ ((L-R)+(L-R)/R*Xalae)*(L<=Xloss)
> mean(Z)
[1] 12596.45```

Now, play with it… it is possible to find a better fit, I guess…

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.