Maximum likelihood estimates for multivariate distributions

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  denotes losses, and  the allocated expenses, a standard excess treaty can be has payoff

where  denotes the (upper) limit, and  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…