I will be back in Guanajuato, Mexico, this week, to visit Victor Rivero. And I will give a talk at the Centro de Investigacion en Matematicas (CIMAT) this Wednesday on “Multivariate Archimax Copulas“. The slides are already online.
(there is a lot of material on copulas, as requested, to provide an introduction for students not familiar with this concept).
We’ve seen Monday, in the MAT8595 course how to use the Generalized Pareto Distribution to estimate some downside risk measures, given a sample (assumed to be i.i.d., I will not mention here properties on extremes for stochastic processes) with distribution . The cumulative distribution function of the Pareto distribution is here
For some threshold , and
, we can write
From Pickands–Balkema–de Haan theorem, if is large enough, then
Given our sample , let
denote the number of observations over, threshold
. Then we can write
or equivalently
If we invert this function, we get the quantile of level ,
Actually, a threshold and then the implied number of observation exceeding that threshold, it is possible to consider a fixed number of observation, and then the associated threshold will be the associated order statistics.
The density of the Pareto distribution is here
which is here function of two paramters, and
.As discussed in the course, it is possible to use the Delta method to derive the asymptotic distribution of any quantile, and get then an approximated (asymptotic) confidence interval.
But since is usually not a parameter of interest, why not considering a reparametrization of our density, as a function of
and
(for some probability
that will be considered as fixed from now on). We can easily get (assuming that
) that
Tis expression is simple, and can be used to derive the likelihood (on the observations exceeding the threshold)
Numerically, let us write (and plot) that function. Consider some real data here
> X=as.numeric(danish) > Xs=sort(X,decreasing=TRUE) > n=length(X) > u=10 > nu=sum(X>u)
Consider, say, the 99.9% quantile,
> p=.999
The empirical quantile is here
> quantile(X,p) 99.9% 131.5519
The density and the loglikelihood functions are here
> gq=function(x,xi,q){ + ( (n/nu*(1-p) ) ^ (-xi)-1)/(xi*(q-u))* + (1+((n/nu*(1-p))^(-xi)-1)/(q-u)*x)^(-1/xi-1)} > loglik=function(param){ + xi=param[2];q=param[1] + lg=function(i) log(gq(Xs[i],xi,q)) + return(-sum(Vectorize(lg)(1:nu))) + }
We can try to plot this likelihood using
> h=201 > Q=seq(50,300,length=h) > XI=seq(.1,1,length=h) > XIQ=as.matrix(expand.grid(Q,XI)) > M=mapply(loglik,XIQ)
Unfortunately, it was not working, so I used the old style
> M=matrix(NA,h,h) > for(i in 1:h){for(j in 1:h){M[i,j]=loglik(c(Q[i],XI[j]))}}
The level curves of the log-likelihood are here
> hc=heat.colors(100) > image(Q,XI,-M,col=hc) > contour(Q,XI,-M,add=TRUE)
Again, since our interest is in the quantile, we can draw the profile likelihood and get the maximum of that function
> PL=function(Q){ + profilelikelihood=function(xi){ + loglik(c(Q,xi))} + return(optim(par=.8,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(100,500))) $minimum [1] 111.1055
and the graph is
> XQ=seq(50,300,length=101) > L=Vectorize(PL)(XQ) > plot(XQ,-L,type="l") > up=OPT$objective > abline(h=-up) > abline(h=-up-qchisq(p=.95,df=1),col="red") > I=which(-L>=-up-qchisq(p=.95,df=1)) > lines(XQ[I],rep(-up-qchisq(p=.95,df=1),length(I)), + lwd=5,col="red") > abline(v=range(XQ[I]),lty=2,col="red")
which can be seen as an alternative to
> gpd.q(tailplot(gpd(X,u)),.999) Lower CI Estimate Upper CI 64.66184 94.28956 188.91752 $objective [1] 454.6481
If we want to focus on another downside risk measure, that shouldn’t be too difficult. For instance, the expected shortfall, can be estimated as
where denotes the mean excess function, which can be writen, with a Generalized Pareto Distribution
Thus, a natural estimator for the expected shortfall is
One more time, it is possible to re-parametrize the density of the Pareto distribution, using instead of
. Here, we get
The code to get the associated log-likelihood is here
> ge=function(x,xi,es){ + (xi+(n/nu*(1-p))^(-xi)-1)/(xi*(1-xi)*(es-u))*(1+(xi+(n/nu*(1-p))^(-xi) + -1)/((es-u)*(1-xi))*x)^(-1/xi-1) + } > loglik=function(param){ + xi=param[2];es=param[1] + lg=function(i) log(ge(Xs[i],xi,es)) + return(-sum(Vectorize(lg)(1:nu))) + }
and again, we can plot it
and the profile (log) likelihood is here (for the 99.9% expected shortfall)
> PL=function(ES){ + profilelikelihood=function(xi){ + loglik(c(ES,xi))} + return(optim(par=.8,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(100,500))) $minimum [1] 143.66 $objective [1] 454.6481
which could be compared with
> gpd.sfall(tailplot(gpd(X,u)),.999) Lower CI Estimate Upper CI 96.64625 191.36972 394.87555
In the MAT8595 course, we’ve seen yesterday Hill estimator of the tail index. To be more specific, we did see see that if , with
, then Hill estimators for
are given by
for . Then we did say that
satisfies some consistency in the sense that
if
, but not too fast, i.e.
(under additional assumptions on the rate of convergence, it is possible to prove that
). Further, under additional technical conditions
In order to illustrate this point, consider the following code. First, let us consider a Pareto survival function, and the associated quantile function
> alpha=1.5
> S=function(x){ifelse(x>1,x^(-alpha),1)}
> Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root}
The code here is obviously too complicated, since this power function can easily be inverted. But later on, we will consider a more complex survival function. Here are the survival function, and the quantile function,
> u=seq(0,5,by=.01) > plot(u,Vectorize(S)(u),type="l",col="red") > u=seq(0,99/100,by=.01) > plot(u,Vectorize(Q)(u),type="l",col="blue",ylim=c(0,20))
 |
 |
Here, we need the quantile function to generate a random sample from this distribution,
> n=500 > set.seed(1) > X=Vectorize(Q)(runif(n))
Hill plot is here
> library(evir) > hill(X) > abline(h=alpha,col="blue")
We can now generate thousands of random samples, and see how those estimators behave (for some specific ‘s).
> ns=10000
> HillK=matrix(NA,ns,10)
> for(s in 1:ns){
+ X=Vectorize(Q)(runif(n))
+ H=hill(X,plot=FALSE)
+ hillk=function(k) H$y[H$x==k]
+ HillK[s,]=Vectorize(hillk)(15*(1:10))
+ }
and if we compute the average,
> plot(15*(1:10),apply(HillK,2,mean)
we do get a series of estimators that can be considered as unbiased.
So far, so good. Now, recall that being in the max-domain of attraction of the Fréchet distribution does not mean that , with
, but is means that
for some slowly varying function , not necessarily constant! In order to understand what could happen, we have to be slightly more specific. And this can be done only by looking at second order regular variation property of the survival function. Assume, here that there is some auxilary function
such that
This (positive) constant is – somehow – related to the speed of convergence of the ratio of the survival functions to the power function (see e.g. Geluk et al. (2000) for some examples).
To be more specific, assume that
then, the second order regular variation property is obtained using , and then, if
goes to infinity too fast, then the estimator will be biased. More precisely (see Chapter 6 in Embrechts et al. (1997)), if
, then, for some
,
The intuitive interpretation of this result is that if is too large, and if the underlying distribution is not exactly a Pareto distribution (and we do have this second order property), then Hill’s estimator is biased. This is what we mean when we say
(the later comes from properties of a sample mean: the more observations, the less the volatility of the mean).
Let us run some simulations to get a better understanding of what’s going on. Using the previous code, it is actually extremly simple to generate a random sample with survival function
> beta=.5
> S=function(x){+ ifelse(x>1,.5*x^(-alpha)*(1+x^(-beta)),1) }
> Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root}
If we use the code above. Here, with
> n=500 > set.seed(1) > X=Vectorize(Q)(runif(n))
the Hill plot becomes
> library(evir) > hill(X) > abline(h=alpha,col="blue")
But it’s based on one sample, only. Again, consider thousands of samples, and let us see how Hill’s estimator is behaving,
so that the (empirical) mean of those estimator is
This week, in the MAT8595 course, we will start the section on inference for extreme values. To start with something simple, we will use maximum likelihood techniques on a Generalized Pareto Distribution (we’ve seen Monday Pickands-Balkema-de Hann theorem).
In the context of parametric models, the standard technique is to consider the maximum of the likelihood (or the log-likelihod).Let denote the parameter (with
). Given some – stnardard – technical assumptions, such as
, or
on some neighbourhood of
, then
where denotes Fisher information matrix (see any textbook for mathematical statistics courses). Consider here some i.i.d. sample, from a Generalized Pareto Distribution, with parameter
, so that
If we solve (numerically) the first order condition of the maximum likelihood, we get an estimator which satisfies
The idea of this asymptotic normality is the following : if the true distribution of the sample is a GPD with parameter , then, if
is large enough, then
will have a joint normal distribution. So if we generate a lot of sample (sufficently large, say 200 observations), then the scatterplot of the estimator should the same as the scatterplot of a Gaussian distribution,
> library(evir)
> n=200
> param=matrix(NA,1000,2)
> for(s in 1:1000){
+ x=rgpd(n,xi=1/1.5,beta=1)
+ param[s,]=gpd(x,0)$par.ests
+ }
> m=apply(param,2,mean)
> S=var(param)
> library(mnormt)
> x=seq(min(param[,1])-.05,max(param[,1])+.05,length=101)
> y=seq(min(param[,2])-.05,max(param[,2])+.05,length=101)
> vx=rep(x,each=length(y))
> vy=rep(y,length(x))
> vz=dmnorm(cbind(vx,vy),m,S)
> z=matrix(vz,length(y),length(x))
> COL=rev(heat.colors(100))
> image(x,y,z,col=COL)
> points(param)
and to get a 3d representation
> x=seq(min(param[,1])-.05,max(param[,1])+.05,length=31) > y=seq(min(param[,2])-.05,max(param[,2])+.05,length=31) > vx=rep(x,each=length(y)) > vy=rep(y,length(x)) > vz=dmnorm(cbind(vx,vy),m,S) > z=matrix(vz,length(y),length(x)) > persp(x,y,t(z),shade=TRUE,col="green",theta=-30,phi=20,ticktype="detailed", + xlab="xi",ylab="sigma")
 |
 |
With 200 observations, if the true underlying distribution is a GPD, then, indeed, the joint distribution of seems to be normal. That would be interesting to generate some confidence intervals for instance, or define some tests.
To go further, see any standard textbook on statistical mathematics, e.g. Casella & Berger (2002).
Another important property is the so called delta-method (we’ve seen Monday in class that it was obtained easily using a first order Taylor expansion). The idea is that if is asymptotically normal, and if
is sufficently smooth, then
will also be asymptotically Gaussian. More precicely (see also the header of this blog)
From this property, we can get the normality of (which is another parametrization used in extreme value models), or on any quantile,
. Let us run some simulation, one more time to check that we actually have a joint normality.
> library(evir)
> n=200
> param=riskm=matrix(NA,1000,2)
> for(s in 1:1000){
+ x=rgpd(n,xi=1/1.5,beta=1)
+ param[s,]=gpd(x,0)$par.ests
+ xihat=param[s,1]
+ shat=param[s,2]
+ q=shat * (.01^(-xihat) - 1)/xihat
+ tvar=q+(shat + xihat * q)/(1 - xihat)
+ riskm[s,]=c(1/xihat,q)
+ }
> m=apply(riskm,2,mean)
> S=var(riskm)
> library(mnormt)
> x=seq(min(riskm[,1])-.05,max(riskm[,1])+.05,length=101)
> y=seq(min(riskm[,2])-.05,max(riskm[,2])+.05,length=101)
> vx=rep(x,each=length(y))
> vy=rep(y,length(x))
> vz=dmnorm(cbind(vx,vy),m,S)
> z=matrix(vz,length(y),length(x))
> image(x,y,t(z),col=COL)
> points(riskm)
As we can see bellow, with samples of size 200, we cannot use this asymptotical result: it looks like we do not have enought data. But if we run the same code with
> n=5000
 |
 |
We get the joint normality of and
. This is what we can get from this result, called delta-method in statistical textbooks. See again Casella & Berger (2002) for more details.
Another interesting tool is the concept of profile likelihood. This would be interesting here since the main interest is the tail index ,
being here some kind of auxilary parameter. See Venzon & Moolgavkar (1988) for more details. Here, we will plot
But more generally, it is possible to consider
where 
is the set of interesting parameters. Then (under standard suitable conditions) we can prove that
which can be used to derive confidence intervals. In the GPD case, for each , we have to find an optimal
. We compute the (profile) likelihood i.e.
. And we can compute the maximum of this profile likelihood. This two-stage optimization is, in general, not equivalent with the (global) maximization of the likelihood, as computed below
> n=500
> set.seed(1)
> x=rgpd(n,xi=1/1.5,beta=1)
> loglikelihood=function(xi,beta){
+ sum(log(dgpd(x,xi,mu=0,beta))) }
> XIV=(1:300)/100;L=rep(NA,300)
> for(i in 1:300){
+ XI=XIV[i]
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ L[i]=-optim(par=1,fn=profilelikelihood)$value }
> plot(XIV,L,type="l")
> XIV[which.max(L)]
[1] 0.67
> gpd(x,0)$par.ests
xi beta
0.6730145 0.9725483
We are not far away. Actually, if we want to compute the maximum of the profile likelihood (and not only compute the values of the profile likelihood on a grid, as before), we use
> PL=function(XI){
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ return(optim(par=1,fn=profilelikelihood)$value)}
> (OPT=optimize(f=PL,interval=c(0,3)))
$minimum
[1] 0.6731025
$objective
[1] 822.5574
Observe that, indeed, we are not far away from the maximum likelihood estimator of (I believe that it’s mainly a computational issue here, and theat the two are similar, here… actually, I’d be glad to hear about cases where maximum of the profile likelihood is not the same as the maximum of the likelihood). The interesting point is that we can use this technique to compute a confidence interval, and even visualize it on a graph
> up=OPT$objective > abline(h=-up) > abline(h=-up-qchisq(p=.95,df=1),col="red") > I=which(L>=-up-qchisq(p=.95,df=1)) > lines(XIV[I],rep(-up-qchisq(p=.95,df=1),length(I)), + lwd=5,col="red") > abline(v=range(XIV[I]),lty=2,col="red")
The vertical lines are the lower and the upper bound of a 95% confidence interval for parameter .
To go further, see Murphy, S.A & van der Vaart, A.W. (2000). On Profile Likelihood.
A reinsurance case study for tomorrow’s class. The goal will be to price some nonproportional reinsurance contract, for business interruption claims. Consider the following dataset,
> library(gdata) > db=read.xls( + "https://perso.univ-rennes1.fr/arthur.charpentier/SIN_1985_2000-PE.xls", + sheet=1) Content type 'application/vnd.ms-excel' length 183808 bytes (179 Kb) open URL ================================================== downloaded 179 Kb
As for any (standard) insurance contract, there are two parts in the pricing
Here, we do not have covariates (but it might be possible to use some, like the kind of industry, the location, etc).
Let us start with the expected number of claims, per year. Here is the daily frequency,
The data are rather old… but somehow, it is a good thing since after ten years, we can expect that most of the claims have been settled (we’ll discuss claims dynamic starting next week). To plot the graph above, we use
> date=db$DSUR > D=as.Date(as.character(date),format="%Y%m%d") > vD=seq(min(D),max(D),by=1) > sD=table(D) > d1=as.Date(names(sD)) > d2=vD[-which(vD%in%d1)] > vecteur.date=c(d1,d2) > vecteur.cpte=c(as.numeric(sD),rep(0,length(d2))) > base=data.frame(date=vecteur.date,cpte=vecteur.cpte) > plot(vecteur.date,vecteur.cpte,type="h",xlim=as.Date(as.character( + c(19850101,20111231)),format="%Y%m%d"))
Then, we can get a prediction of the daily number of business interruption claims, e.g. for any day in 2010 (assume that we had to price a reinsurance contract a few years ago), using a (standard) Poisson regression
> regdate=glm(cpte~date,data=base,family=poisson(link="log")) > nd2010=data.frame(date=seq(as.Date(as.character(20100101),format="%Y%m%d"), + as.Date(as.character(20101231),format="%Y%m%d"),by=1)) > pred2010 =predict(regdate,newdata=nd2010,type="response") > sum(pred2010) [1] 159.4757
Observe that using old data has drawbacks, since we got much more uncertainty if we use a regression on time (to include some possible trend)
Say we have something like 160 claims over a given year, on average.
> plot(D,db$COUTSIN,type="h")
Let us now focus on the cost of those claims. We have 2,400 claims in our dataset, to fit a model (or at least estimate how much a reinsurance contract might cost us). Assume that we would like to purchase a reinsurance contract for our very large claims. Like the two largest per year. Over 16 years, the decutible should be close to the cost of the 32nd largest claim, which was close to 15 million.
> quantile(db$COUTSIN,1-32/2400)/1e6 98.66667% 15.34579 > abline(h=quantile(db$COUTSIN,1-32/2400),col="blue")
So consider some reinsurance contract with a deductible of 15 million. Unfortunately, we cannot find unlimited covers. So let us assume that a reinsurance company agrees for such a deductible, but with a limited cover of 35 million. The average cost (for the reinsurance company) is where
A first idea is to look at the first cost, i.e. the empirical average of that indemnity, on our portfolio. The indemnity function is
> indemn=function(x) pmin((x-15)*(x>15),50-15)
we can check on a few losses that it is actually what we wish to compute
> indemn(5) [1] 0 > indemn(20) [1] 5 > indemn(50) [1] 35
Now, if the compute the average repayment by the reinsurance company, over 16 years, we get
> mean(indemn(db$COUTSIN/1e6)) [1] 0.1624292
So, per claim, the reinsurance company will pay, on average 162,430. With 160 claims per year, the pure premium should be close to 26 million
> mean(indemn(db$COUTSIN/1e6))*160 [1] 25.98867
(again, for a 35 million cover, for some claims that should occur, on average, twice a year). As we will see, a standard model in reinsurance is the Pareto distribution (or to be more specific, a Generalized Pareto one),
There are three parameters here
The strategy is to consider a threshold below our deductible, e.g. 12 million. Then, given that the loss exceed 12 million, we can fit a Generalized Pareto distribution,
> gpd.PL <- gpd(db$COUTSIN,12e6)$par.ests
> gpd.PL
xi beta
7.004147e-01 4.400115e+06
and compute
> E <- function(yinf,ysup,xi,beta,threshold){
+ as.numeric(integrate(function(x) (x-yinf)*dgpd(x,xi,mu=threshold,beta),
+ lower=yinf,upper=ysup)$value+
+ (1-pgpd(ysup,xi,mu=threshold,beta))*(ysup-yinf))
+ }
Here, given that a claim exceeds 12 million, the average repayment is close to 6 million
> E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 6058125
Now, we have to take into account the probability to reach 12 million, which is here
> mean(db$COUTSIN>12e6) [1] 0.02639296
So, if we summarize, we have on average 160 claims per year,
> p [1] 159.4757
Only 2.6% will exceed 12 million
> mean(db$COUTSIN>12e6) [1] 0.02639296
So, the yearly frequency of claism larger than 12 million is 4.2 claims
> p*mean(db$COUTSIN>12e6) [1] 4.209036
And for a claim that exceed 12 million, the average repayment is
> E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 6058125
So, the pure premium should be close to
> p*mean(db$COUTSIN>12e6)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 25498867
which (hopefully) is close to the empirical value we got. Actually, it is also possible to look at the impact of the threshold parameter, since it is clearly and intermediate value that could be changed. I mean, why 12 and not 10? Consider
> esp=function(threshold=12e6,p=sum(pred2010)){
+ (gpd.PL <- gpd(db$COUTSIN,threshold)$par.ests)
+ return(p*mean(db$COUTSIN>threshold)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],threshold))
+ }
We can plot the pure premium as a function of that threshold,
> seuils=seq(1e6,15e6,by=1e6) > plot(seuils,Vectorize(esp)(seuils),type="b",col="red")
which is between 24 and 26 for large thresholds. Again, that is only the first step, and we can price a higher reinsurance layer, like a reinsurance contract with a deductible of 50 million (we have our previous reinsurance contract for claims below that threshold), and a cover of 50 million, for instance. For those high layers, it become interesting to have a parametric model, which should be more robust than the empirical average.
In order to illustrate the next section of the non-life insurance course, consider the following example1, inspired from http://sciencepolicy.colorado.edu/…. This is the so-called “Normalized Hurricane Damages in the United States” dataset, for the period 1900-2005, from Pielke et al. (2008). The dataset is available in xls format, so we have to spend some time to import it,
> library(gdata)
> db=read.xls(
+ "http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls",
+ sheet=1)
trying URL 'http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls'
Content type 'application/vnd.ms-excel' length 119296 bytes (116 Kb)
opened URL
==================================================
downloaded 116 Kb
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
LANGUAGE = "fr_CA:fr",
LC_ALL = (unset),
LANG = "fr_CA.UTF-8"
are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
The problem with excel spreadsheets is that some columns might have pre-specified format (here, losses are with a format 000,000,000 for instance)
> tail(db)
Year Hurricane.Description State Category Base.Economic.Damage
202 2005 Cindy LA 1 320,000,000
203 2005 Dennis FL 3 2,230,000,000
204 2005 Katrina LA,MS 3 81,000,000,000
205 2005 Ophelia NC 1 1,600,000,000
206 2005 Rita TX 3 10,000,000,000
207 2005 Wilma FL 3 20,600,000,000
Normalized.PL05 Normalized.CL05 X X.1
202 320,000,000 320,000,000 NA NA
203 2,230,000,000 2,230,000,000 NA NA
204 81,000,000,000 81,000,000,000 NA NA
205 1,600,000,000 1,600,000,000 NA NA
206 10,000,000,000 10,000,000,000 NA NA
207 20,600,000,000 20,600,000,000 NA NA
To get data in a format we can play with, consider the following function,
> stupidcomma = function(x){
+ x=as.character(x)
+ for(i in 1:10){x=sub(",","",as.character(x))}
+ return(as.numeric(x))}
and let’s convert those values into numbers,
> base=db[,1:4] > base$Base.Economic.Damage=Vectorize(stupidcomma)(db$Base.Economic.Damage) > base$Normalized.PL05=Vectorize(stupidcomma)(db$Normalized.PL05) > base$Normalized.CL05=Vectorize(stupidcomma)(db$Normalized.CL05)
Here is the dataset we will use, from now on,
> tail(base)
Year Hurricane.Description State Category Base.Economic.Damage
202 2005 Cindy LA 1 3.20e+08
203 2005 Dennis FL 3 2.23e+09
204 2005 Katrina LA,MS 3 8.10e+10
205 2005 Ophelia NC 1 1.60e+09
206 2005 Rita TX 3 1.00e+10
207 2005 Wilma FL 3 2.06e+10
Normalized.PL05 Normalized.CL05
202 3.20e+08 3.20e+08
203 2.23e+09 2.23e+09
204 8.10e+10 8.10e+10
205 1.60e+09 1.60e+09
206 1.00e+10 1.00e+10
207 2.06e+10 2.06e+10
We can visualize the normalized costs of hurricanes, from 1900 till 2005, with the 207 hurricanes (here the x-axis is not time, it is simply the index of the loss)
> plot(base$Normalized.PL05/1e9,type="h",ylim=c(0,155))
As usual, there are two components when computing the pure premium of an insurance contract. The number of claims (or here hurricanes) and the individual losses of each claim. We’ve seen – above – individual losses, let us focus now on the annual frequency.
> TB <- table(base$Year) > years <- as.numeric(names(TB)) > counts <- as.numeric(TB) > years0=(1900:2005)[which(!(1900:2005)%in%years)] > db <- data.frame(years=c(years,years0), + counts=c(counts,rep(0,length(years0)))) > db[88:93,] years counts 88 2003 3 89 2004 6 90 2005 6 91 1902 0 92 1905 0 93 1907 0
On average, we experience about 2 (major) hurricanes per year,
> mean(db$counts) [1] 1.95283
In predictive modeling (here, we wish to price a reinsurance contract for, say, 2014), we need probably to take into account some possible trend in the hurricane occurrence frequency. We can consider either a linear trend,
> reg0 <- glm(counts~years,data=db,family=poisson(link="identity"), + start=lm(counts~years,data=db)$coefficients)
or an exponential one,
> reg1 <- glm(counts~years,data=db,family=poisson(link="log"))
We can plot those three predictions, and get a prediction for the number of (major) hurricanes in 2014,
> plot(years,counts,type='h',ylim=c(0,6),xlim=c(1900,2020))
> cpred1=predict(reg1,newdata=data.frame(years=1890:2030),type="response")
> lines(1890:2030,cpred1,col="blue")
> cpred0=predict(reg0,newdata=data.frame(years=1890:2030),type="response")
> lines(1890:2030,cpred0,col="red")
> abline(h=mean(db$counts),col="black")
> (predictions=cbind(constant=mean(db$counts),linear=
+ cpred0[126],exponential=cpred1[126]))
constant linear exponential
126 1.95283 3.573999 4.379822
> points(rep((1890:2030)[126],3),prediction,col=c("black","red","blue"),pch=19)
Observe that changing the model will change the pure premium: with a flat prediction, we expect less than 2 (major) hurricanes, but with the exponential trend, we expect more than 4…
This is for the expected frequency. Now, we should find a suitable model to compute the pure premium of a reinsurance treaty, with a (high) deductible, and a limited (but large) cover. As we will seen in class next week, the appropriate model is a Pareto distribution (see Hagstrœm (1925), Huyghues-Beaufond (1991) or a survey – in French – published a few years ago).
We can use Hill’s plot to estimate the tail index,
> library(evir) > hill(base$Normalized.PL05)
Clearly, costs of major hurricanes are heavy tailed.
Now, consider an insurance company, in the U.S., with 5% market share (just to illustrate). We will consider there \tilde Y_i= Y_i/20. The losses are given below. Consider a reinsurance treaty, with a deductible of 2 (billion) and a limited cover of 4 (billion),
For our Pareto model, consider only losses above 500 millions,
> threshold=.5
> (gpd.PL <- gpd(base$Normalized.PL05/1e9/20,threshold)$par.ests)
xi beta
0.4424669 0.6705315
Keep in mind the 1 hurricane out of 8 reaches that level
> mean(base$Normalized.CL05/1e9/20>.5) [1] 0.1256039
Given that the loss exceeds 500 millions, we can now compute the expected value of the reinsurance contact,
To compute it we can use
> E <- function(yinf,ysup,xi,beta){
+ as.numeric(integrate(function(x) (x-yinf)*dgpd(x,xi,mu=threshold,beta),
+ lower=yinf,upper=ysup)$value+
+ (1-pgpd(ysup,xi,mu=threshold,beta))*(ysup-yinf))
+ }
[Nov 5th] there is a typo in the previous function, since the threshold should be used, here, as a parameter in the function, if you want to play with that function an see the impact of the threshold (see a more recent post on the same topic, but a different dataset)… but here, we do not change the threshold, so it is not a big deal.
Now, it is probably time to bring all the pieces together. We might expect a bit less than 2 (major) hurricanes per year,
> predictions[1] [1] 1.95283
and each hurricane has 12.5% chances to cost more than 500 million for our insurance company,
> mean(base$Normalized.PL05/1e9/20>.5) [1] 0.1256039
and given that an hurricane exceeds 500 million loss, then the expected repayment by the reinsurance company is (in millions)
> E(2,6,gpd.PL[1],gpd.PL[2])*1e3 [1] 330.9865
So the pure premium of the reinsurance contract is simply
> predictions[1]*mean(base$Normalized.PL05/1e9/20>.5)* + E(2,6,gpd.PL[1],gpd.PL[2])*1e3 [1] 81.18538
for a cover of 4 billion, in excess of 2.
1.This example will be found in the Reinsurance and Extremal Events chapter in the forthcoming Computational Actuarial Science with R, by Eric Gilleland and Mathieu Ribatet.
Last week, we had a discussion with some colleagues about the fact that – in order to prepare for the SOA exams – we did not have time (so far) to mention results on extreme values in our actuarial program. I did gave an introduction in my nonlife actuarial models class, but it was only an introduction, in three hours, in order to illustrate reinsurance pricing. And I told my students that if they wanted to know more about extreme values, they should start a master program in actuarial science and finance, since I will give a course on extremes (and copulas) next winter.
But actually, extreme values are everywhere ! For instance, there is a Prudential TV commercial where has people place large, round stickers on a number line to represent the age of the oldest person they know. This forms some kind of histogram. The message is to have Prudential prepare you to have adequate money for all these years. And actually, anyone can add his or her own sticker at the Prudential website.
Patrick Honner, on his blog (http://mrhonner.com/…), did mention this interesting representation. But this idea is not new, as mentioned in a post, published three years ago. In 1932, Emil Gumbel gave a talk in France on the “âge limite“. And as he wrote it “on peut donc supposer que la distribution de l’âge limite – c’est à dire la probabilité que cet âge ait une valeur donnée – soit Gaussienne“. In 1932 (not aware of Fisher and Tippett work, he thought that the limiting distribution for a maximum would be Gaussian). But a few years after, he read about Fisher’s work, and observed also that “la distribution d’une valeur extrêmes peut être représentée pour un nombre suffisant d’observations par la formule doublement exponentielle, pourvu que la distribution initiale se comporte asymptotiquement comme une exponentielle. La formule devient rigoureuse si la distribution initiale est exponentielle“, as he wrote in 1935. And in 1937, he wrote a paper on “les centennaires” that can also be related to the work of Bortkiewicz on rare events. One should also mention one of the most important paper in extreme value theory, published in 1974 by Balkema and de Haan, on Residual Life Time at Great Age.
Because in this experiment, the question is “How Old is the Oldest Person You Know?“, so it is the distribution of a maximum. And from Fisher-Tippett theorem, if we assume that the age is bounded (and that there exists some finite upper limit), then the limiting distribution for the maxima (or to be more rigorous, a affine transformation of the maxima) should be Weibull distribution. And this is what it looks like
> plot(-x,dweibull(x,2.25,4),type="l",lwd=2)
As an actuary, the only thing I know about demography, is the distribution of the age of death. For instance, consider the following French life table
> alive <- read.table( + "https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv", + sep=";",header=TRUE)$Lx > nb= -diff(alive) > ages=0:110 > plot(ages,nb,type="h")
This is the distribution of the age of the death in a given population. Which is not the same as the distribution mentioned above! What we look for is the following: given that someone is alive, what could be the distribution of his-her age ? Actually, if we assume that the yearly number of birth is constant with time (as well as death probability), then we can compute easily to number of people of age : we take everyone born (exactly)
years ago, and remove all those who died at at
,
, etc. So the function should be
> probadeath=nb/sum(nb) > nbx=function(x) 1-sum(probadeath[1:(x+1)]) > surv=Vectorize(nbx)(ages) > distrage=surv/sum(surv)
which looks like
But this assumption of constant number of birth is not that relevent. And actually, what we need is the distribution of the age within a population… This is a population pyramid, actually. The French one can be downloaded from http://www.insee.fr/fr/ppp/bases-de-donnees/….
> population <- read.table("popinsee2007.csv",sep=";",header=TRUE)$POPTOT07
> ages=0:107
> plot(ages,population/sum(population),type="h")
(the red line being the one obtained previously, using some natality assumptions). Now, let us use this population to generate acquaintances.
> agemax=function(nsim=1000,size=20){
+ agemax=rep(NA,nsim)
+ for(i in 1:nsim){
+ X=sample(ages,prob=population/sum(population),size=size,replace=TRUE)
+ agemax[i]=max(X)}
+ return(agemax)}
Here, we assume that everyone knows 20 other people, randomly chosen in the entire population, then we return the age of the oldest. And we do that for 1,000 people. Here is the distribution, we obtain
> XS=agemax(10000,20) > plot(table(XS)/length(XS),type="h",xlim=c(0,108))
where the red line is a Weibull distribution (a transformed one, actually, since in extremely value theory, the distance to the upper bound of the distribution has a Weibull density),
> library(MASS) > fit=fitdistr(108-XS,dweibull,list(shape=1,scale=1)) > lines(ages,dweibull(108-ages,fit$estimate[1],fit$estimate[2]),col="red")
Which is quite close to the distribution obtained in the commercial, don’t you think ? But still, it should be possible to be more accurate, since people should think of their parents, or grandparents. So I guess it could be possible to build a more accurate algorithm, to get something closer to the distribution obtained on the Prudential website. But first, let us wait to have more stickers, more observations… and then I’ll be back to play with it !
A few months ago, I did publish a post (in French) following my reading of Leonard Mlodinow’s the Drunkard’s Walk. More precisely, I mentioned a paragraph that I found extremely informative
But it looks like those gangsters were not only stealing money. They were also stealing ideas, here from a British statistician, manely Leonard Henry Caleb Tippett. Leonard Tippett is famous in Extreme Value Theory for his theorem (the so-called Fisher-Tippett theorem, which gives the possible limiting distributions for a normalized version of the maximum from an i.i.d. sequence, see old posts). According to Martin Gardner, Leonard Tippett suggested to use middle numbers (not the last ones) of larger ones to generate (pseudo) random sequences, or more precisely, in 1927, “published a table of 41,600 random numbers, obtained by taking the middle digits of the area of parishes in England”
I could not get a copy of the book Random Sampling Numbers by Leonard Tippett (I could only find reviews, e.g. Nair (1938)) but I do believe that this technique should work to generate sequences that do look like sequences of random numbers. Note that several techniques were mentioned in previous posts (in French) published a few years ago.
Now, I should also take some time to apologize because, sometimes, I am the one playing the gangster: I do steal a lot of illustrations on the internet. And I would like to apologize to the authors. On my previous blog, I did try – once – to add a short line at the end of a post, explaining where the illustration was coming from (trying to give credit to the illustrator). Less than 10 days after adding this short line, I received an email from a ‘publisher’, telling me that there were rights attached to the picture, and that I had 24 hours to remove it (if not, their lawyers will see what to do). Of course, I did remove the picture, and the mention. Now, I use pictures, and no mention. And I feel guilty. So I wanted to apologize for stealing others’ work. I am still discussing to hire an illustrator, to illustrate my blog. Work in progress….
During the course, we have seen that it is natural to assume that not only the individual claims frequency can be explained by some covariates, but individual costs too. Of course, appropriate families should be considered to model the distribution of the cost , given some covariates
.Here is the dataset we’ll use,
> sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt", + header=TRUE,sep=";") > sinistres=sinistre[sinistre$garantie=="1RC",] > sinistres=sinistres[sinistres$cout>0,] > contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt", + header=TRUE,sep=";") > couts=merge(sinistres,contrat) > tail(couts) nocontrat no garantie cout exposition zone puissance agevehicule 1919 6104006 11933 1RC 5376.04 0.37 E 6 1 1920 6107355 12349 1RC 51.63 0.74 E 4 1 1921 6108364 13229 1RC 1320.00 0.74 B 9 1 1922 6109171 11567 1RC 1320.00 0.74 B 13 1 1923 6111208 14161 1RC 970.20 0.49 E 10 5 1924 6111650 14476 1RC 1940.40 0.48 E 4 0 ageconducteur bonus marque carburant densite region 1919 32 57 12 E 93 10 1920 45 57 12 E 72 10 1921 32 100 12 E 83 0 1922 56 50 12 E 93 13 1923 30 90 12 E 53 2 1924 69 50 12 E 93 13
Here, each line is a claim. Usual families to model the cost are the Gamma distribution, or the inverse Gaussian. Or the lognormal distribution (which is not in the exponential family, but one can assume that the logarithm of the cost can be modeled with a Gaussian distribution). Consider here only one covariate, e.g. the age of the car, and two different models: a Gamma one, and a lognormal one.
> age=0:20 > reggamma.sp <- glm(cout~agevehicule,family=Gamma(link="log"), + data=couts) > Pgamma <- predict(reggamma.sp,newdata=data.frame(agevehicule=age),type="response")
For the Gamma regression, it is a simple GLM, so it is not difficult. For a lognormal distribution, one should remember that the expected value of a lognormal distribution is not the exponential of the underlying Gaussian distribution. A correction should be made, here to get an unbiased estimator for the average cost,
> reglm.sp <- lm(log(cout)~agevehicule,data=baseCOUT) > sigma <- summary(reglm.sp)$sigma > mu <- predict(reglm.sp,newdata=data.frame(agevehicule=age)) > Pln <- exp(mu+sigma^2/2)
We can plot those two predictions on a single graph,
> plot(age,Pgamma,xlab="",ylab="",col="red",type="b",pch=4) > lines(age,Pln,col="blue",type="b")
Here it is,
Observe that it is also possible to use splines, since there might be no reason for the age to appear here in a multiplicative way,
Here, the two models are rather close. Nevertheless, one should remember that the Gamma model can be extremely sensitive to large claims (I mean here really large claims). On the other hand, with the log-transformation for the lognormal model, it seams that this model is less sensitive to large events. Actually, if I use the complete dataset, the regressions are the following,
i.e. with a lognormal distribution, the average cost is decreasing with the age of the car, while it is increasing with a Gamma model. The main reason here is that there is one large (not to say huge) claim in the dataset,
> couts[which.max(couts$cout),] cout exposition zone puissance agevehicule ageconducteur 7842 4024601 0.22 B 9 13 19 marque carburant densite region 7842 2 E 93 24
One young driver got a $ 4 million claim, with a 13 year old car. This is an outliers for the Gamma regression, that clearly influences the estimation (the second largest if only one third of this one). Since there is a clear influence of large claims on the estimation of the average cost, a natural idea might be to remove those large claims. Or perhaps to see them as different from normal claims: normal claims can be explained by some covariates, but perhaps that those large claims should be shared not only within its own class, but within all the insured on the portfolio. To formalize this idea, observe that we can write
where the blue part is associated to normal-sized claims, while large ones correspond to the red part. It is then possible to run three regressions: one on normal sized claims, one on large claims, and one on the indicator of having a large claims, given that a claim occurred. The code here is something like that: a large claim – here – is above $ 10,000 (one has a fix it)
> s= 10000 > couts$normal=(couts$cout<=s) > mean(couts$normal) [1] 0.9818087
which represent 2% of the claims in our dataset.We can run 3 sets of regressions, with smoothed regression on the age of the car. The first one to model large claims individual costs,
> indice = which(couts$cout>s) > mean(couts$cout[indice]) [1] 34471.59 > library(splines) > regB=glm(cout~bs(agevehicule),data=couts, + subset=indice,family=Gamma(link="log")) > ypB=predict(regB,newdata=data.frame(agevehicule=age),type="response") > ypB2=mean(couts$cout[indice])
the second one to model normal claims individual costs,
> indice = which(couts$cout<=s) > mean(couts$cout[indice]) [1] 1335.878 > regA=glm(cout~bs(agevehicule),data=couts, + subset=indice,family=Gamma(link="log")) > ypA=predict(regA,newdata=data.frame(agevehicule=age),type="response") > ypA2=mean(couts$cout[indice])
And finally, a third one, on the probability of having a normal sized claim, given that a claim occurred
> regC=glm(normal~bs(agevehicule),data=couts,family=binomial) > ypC=predict(regC,newdata=data.frame(agevehicule=age),type="response") > regC2=glm(normal~1,data=couts,family=binomial) > ypC2=predict(regC2,newdata=data.frame(agevehicule=age),type="response")
Note that we to have, each time something that can be interpreted either as , or
– i.e. no covariate is considered on the later. On the graph below, we did plot
where Gamma regressions – with splines – are considered for the average costs, while logistic regressions – again with splines – are considered to model probabilities.
(but careful with splines: on borders, since we do not have a lot of observations, the behavior can be… odd. And adjustments should be made to obtain an adequate level of premium). If it is legitimate to assume that normal-sized claims can be explained by some covariates, perhaps large claims (or extremely large ones) are just purely random, i.e. not function of any covariate, at all. I.e.
To go one step further, it might also be possible to assume that not only the size of the claim (given that it is a large one) is not a function of any covariate, but perhaps neither is the probability of having an extremely large claim, too
From the first part, we’ve seen that the distribution considered had an impact on the prediction, and in the second, we’ve seen that the definition of large claims (and how to deal with them) also has an impact. So clearly, actuaries have some leverage when working on ratemaking…
In insurance, the law of large numbers (named loi des grands nombres initially by Siméon Poisson, see e.g. http://en.wikipedia.org/…) is usually mentioned to legitimate large portfolios, because of pooling and diversification: the larger the pool, the more ‘predictable’ the losses will be (in a given period). Of course, under standard statistical assumption, namely finite expected value, and independence (see http://freakonometrics.blog.free.fr/…. for a discussion, in French). Since in insurance, catastrophes are usually rare – and extremely costly – and actuaries might be interested to model occurrence of that small number of events (see e.g. Aldous’ book on that specific topic, that can be downloaded from http://stat.berkeley.edu/…). The theorem behind is sometimes called the law of small numbers (from the book published by Ladislaus Bortkiewicz, but we’ll get back to that story later on, see also Whitaker (1914) http://biomet.oxfordjournals.org/… or the book recently published by Michael Falk, Jürg Hüsler and Rolf-Dieter Reiss).
The so-called Poisson distribution (see http://en.wikipedia.org/…) was introduced by Siméon Poisson in 1837 (in Recherches sur la Probabilité des Jugements en Matière Criminelle et en Matière Civile, Précédées des Règles Générales du Calcul des Probabilités, see http://gallica.bnf.fr/…). But it had been defined more than a century before, by Abraham De Moivre, in 1711, in De Mensura Sortis seu; de Probabilitate Eventuum in Ludis a Casu Fortuito Pendentibus (see e.g. the review in http://www.jstor.org/…). Let denote a counting random variable, then it said to be Poisson distributed if there is
such that
De Moivre obtained that distribution from an approximation of the binomial distribution. Recall that the binomial distribution is a standard distribution in actuarial science, for instance to model the number of deaths among insured. If individual death probabilities are identical, say
, and if deaths are independent events, then
And if and
, then
Again, this is an asymptotic theorem, which is valid when we have a lot of observations (
), but also, the probability of occurrence should be extremely small (since
), which is why to use the term small numbers. Siméon Poisson was not interested by mathematical approximations: his main point was to get a distribution with nice goodness of fit properties for the data he was working on. He wanted to get a better understanding of cours d’assises (jury panel, might be a valid translation of the French term). A jury consists of 12 jurors who voted to determine whether a defendant was guilty. When guilt was predominant, with at least 8 votes against 4, the defendant was convicted (which was 47% of criminal cases). 5 with 7 votes against, the opinion of professional judges was requested (11% of criminal trials again). Using these statistics we can demonstrate that a defendant brought before an assize court is guilty of the order of 68%, and the probability that a juror is not wrong by voting (condemning an innocent or releasing a culprit) was about 54%. He sought to calculate the probability that a defendant is wrongfully convicted, and gets 2%. And 28% of exonerated defendants are in fact guilty. Siméon Poisson introduced this law to get probabilities easily. But the law he considered is central in probability….
The heuristic of the main theorem, related to the Poisson distribution is the following: let denote i.i.d random variables taking values in
(in a general setting, one component can be the time, the other one an upper region of interest, where some stochastic process might be). Let
. If
as
(or
to be a little bit more specific about the assumptions), let
denote the (random variable characterizing) count of events
, then
can be approximated by a Poisson distribution with parameter
.
The heuristic is that if we consider a large number of observations, and if we count how many are in a given (small) region, then the number of such observations is Poisson distributed.
n=1000 X=runif(n)*10-1.5 Y=runif(n)*10-1.5 plot(X,Y,axis=FALSE,cex=.6) u=seq(-1,1,by=.01) v=sqrt(1-u^2) polygon(c(u,rev(u)),c(v,rev(-v)),col="yellow",border=NA) I=(X^2+Y^2)<1 points(X[I],Y[I],cex=.6,pch=19,col="red")
If we run some simulations,
The parameter of the Poisson distribution is the area of the yellow disk, over the area of the square, i.e.
> (lambda=10*pi) [1] 31.41593 > lines(0:60-.5,dpois(0:60,lambda),type="b",col="red")
To get an interpretation related to insurance modeling, let denote an upper layer in a reinsurance contract, i.e.
for some deductible
. Let
‘s denote individual losses. Then the number of claims that hit this upper layer can be modeled with a Poisson distribution. More precisely, if deductible
becomes extremely large (and
), we obtain the point-over-threshold model in extreme value theory (see e.g. http://brale.math.hr/~iugrina/… or http://fire.nist.gov/bfrlpubs/…): if
has a Poisson distribution and, conditionally on
,
are independent identically distributed generalized Pareto random variables, then
has the generalized extreme value distribution. Thus, exceedances models (for rare events) are closely related to Poisson processes.
As mentioned above, the Poisson distribution appears when events occur somehow randomly and independently, over time. It is then natural to study the time between two occurences (or two claims, in an insurance context).
It is neither Siméon Poisson nor De Moivre, but Ladislaus Von Bortkiewicz who first mentioned the Poisson distribution as the law of small numbers. In 1898 (see https://archive.org/…), he studied the number number of soldiers killed by being kicked by a horse, from 1875 till 1894, in 200 corps (more precisely 10 corps over 20 ans).
He did obtain the following distribution (here, the parameter of the Poisson distribution is 0.61, i.e. the average number of death per year)
| number of death per year |
Empirical counts |
Poisson distribution |
| 0 | 109 | 108.67 |
| 1 | 65 | 66.21 |
| 2 | 22 | 20.22 |
| 3 | 3 | 4.11 |
| 4 | 1 | 0.63 |
| 5 and more | 0 | 0.08 |
It is possible to find a lot of cases where the Poisson distribution fits extremely well. For instance, if we consider the number of hurricanes, that landed in Florida after 1850,
| number of hurricanes per year |
empirical frequency |
Poisson frequency |
| 0 | 30 | 27.16 |
| 1 | 48 | 47.99 |
| 2 | 37 | 42.41 |
| 3 | 29 | 24.98 |
| 4 | 8 | 11.03 |
| 5 | 3 | 3.90 |
| 6 | 3 | 1.15 |
| 7 | 1 | 0.29 |
| 8 and more | 0 | 0.08 |
The return period was introduced by Emil Gumbel, in hydrology, to link probabilities and durations (see e.g. http://freakonometrics.blog.free.fr/…). A decennial event has an occurence probability of 1/10. 10 is then the average waiting time before occurence. This does not mean that the event will not occur before 10 years, or has to occur before 10 years. Consider a return period (in years), then the yearly probability of non-occurrence is
.
And the probability of non-occurence over years is then
. It is standard to summarize this property with the following table,
|
return period |
||||||
|
Number of years ( |
10 | 20 | 50 | 100 | 200 | |
| 10 | 65.1% | 40.1% | 18.3% | 9.6% | 4.9% | |
| 20 | 87.8% | 64.2% | 33.2% | 18.2% | 9.5% | |
| 50 | 99.5% | 92.3% | 63.6% | 39.5% | 22.5% | |
| 100 | 99.9% | 99.4% | 86.7% | 63.4% | 39.5% | |
| 200 | 99.9% | 99.9% | 98.2% | 86.6% | 63.3% | |
The diagonal in the table above is extremely interesting. It looks like there is some kind of convergence towards a limiting value (here 63.2%). Indeed, the number of events observed over n years have a Binomial distribution, with probability , which will converge towards the Poisson distribution with parameter 1. The probability of not having a catastrophe is then
, which is equal to 0.632.
The Poisson distribution keeps appearing when computing probabilies of rare events. For instance, the probability to have at least one incident in a nuclear plant in France, over a 50 year period. Assume that the annual probability of an incident in a reactor is small, e.g. 0.05%. Assume further that reactors are independent among them, and in time. The probability to have an incident over 80 reactors in 50 years is (exactly)
Of course, a linear approximation is not correct (even if it was mentioned in some French newspaper, as explained in an old post http://freakonometrics.blog.free.fr/…)
On the other hand
> p=0.00005 > 1-(1-p)^(50*80) [1] 0.1812733 > 1-exp(-50*80*p) [1] 0.1812692
which is the probability that is null when
has a Poisson distribution with parameter
. We clearly see here an application of De Moivre’s approximation in risk management.
Another way of looking at this problem is based on the following idea: given the fact that in 45 years of observations on 450 reactors worldwide (roughly), three major accidents were observed including Three Mile Island (1979) and Fukushima (2011), i.e. the average time between accidents can be estimated at 16 years. For a single reactor, we can assume that the average time to wait before an incident is 450 times 16 years, i.e 7200 years. Or the probability to have one incident, over one year, for one reactor is 1 over 7200 (this is the idea behind the return period concept). If we assume that the arrival of accidents occurs randomly and independently of each other (as defined above) then the number of major accidents observed over a period of 50 years in France follows a Poisson distribution with parameter 50 / (7200/80). Also, the probability of having no major accident over 50 years, with 80 reactors can be estimated by
i.e.
> 1-exp(-50*80/7200) [1] 0.4262466
(keeping in mind all the uncertainty around the estimated waiting time before a major accident to a single reactor!).
Just a short post to get back on results mentioned at the end of the course. Since copulas are obtained using (univariate) quantile functions in the joint cumulative distribution function, they are – somehow – related to the marginal distribution fitted. In order to illustrate this point, consider an i.i.d. sample 
from a Student-t distribution,
library(mnormt) r=.5 n=200 X=rmt(n,mean=c(0,0),S=matrix(c(1,r,r,1),2,2),df=4)
Thus, the true copula is Student-t. Here, with 4 degrees of freedom. Note that we can easily get the (true) value of the copula, on the diagonal
dg=function(t) pmt(qt(t,df=4),mean=c(0,0), S=matrix(c(1,r,r,1),2,2),df=4) DG=Vectorize(dg)
Four strategies are considered here to define pseudo-copula base variates,
Now that we have a sample with margins in the unit square, let us construct the empirical copula,
Let us now compare those four approaches.
X0=cbind((X[,1]-mean(X[,1])/sd(X[,1])), (X[,2]-mean(X[,2])/sd(X[,2]))) Y=pnorm(X0)
Then, the following code is used to compute the value of the empirical copula, on the diagonal,
diagonale=function(t,Z) mean((Z[,1]<=t)&(Z[,2]<=t)) diagY=function(t) diagonale(t,Y) DiagY=Vectorize(diagY) u=seq(0,1,by=.005) dY=DiagY(u)
On the graph below, 1,000 samples of size 200 have been generated. All trajectories are the estimation of the copula on the diagonal. The black plain line is the true value of the copula
Obviously, it is not good at all. Mainly because the distribution of 
can’t be a copula, since margins are not even uniform on the unit interval.
U=pt(X,df=4)
This time, the fit is much better.
F1=fitdistr(X0[,1],dt,list(df=5),lower = 0.001) F2=fitdistr(X0[,2],dt,list(df=5),lower = 0.001) V=cbind(pt(X0[,1],df=F1$estimate),pt(X0[,2],df=F2$estimate))
Here, it is also very good. Even better than before, when the true distribution is considered.
(it is like using Lillie test for goodness of fit, versus Kolmogorov-Smirnov, seehere for instance, in French).
R=cbind(rank(X[,1])/(n+1),rank(X[,2])/(n+1))
Here it is even better then the previous one
If we compare Box-plots of the value of the copula at point (.2,.2), we obtain the following, with on top ranks, then fitting with the good family, then using the true distribution, and finally, using a non-proper distribution.
Just to illustrate one more time a result mentioned in a previous post, “in statistics, having too much information might not be a good thing“.
Since several students got the intuition that natural catastrophes might be non-insurable (underlying distributions with infinite mean), I will post some comments on testing procedures for extreme value models.
A natural idea is to use a likelihood ratio test (for composite hypotheses). Let 
denote the parameter (of our parametric model, e.g. the tail index), and we would like to know whether is smaller or larger than 
(where in the context of finite versus infinite mean 
). I.e. either belongs to the set 
or to its complementary 
. Consider the maximum likelihood estimator 
, i.e.
Let 
and 
denote the constrained maximum likelihood estimators on 
and 
respectively,
Either 
and 
(on the left), or 
and 
(on the right)
So likelihood ratios
are either equal to
or
If we use the code mentioned in the post on profile likelihood, it is easy to derive that ratio. The following graph is the evolution of that ratio, based on a GPD assumption, for different thresholds,
> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> library(evir)
> X=base1$Loss.in.DKM
> U=seq(2,10,by=.2)
> LR=P=ES=SES=rep(NA,length(U))
> for(j in 1:length(U)){
+ u=U[j]
+ Y=X[X>u]-u
+ loglikelihood=function(xi,beta){
+ sum(log(dgpd(Y,xi,mu=0,beta))) }
+ XIV=(1:300)/100;L=rep(NA,300)
+ for(i in 1:300){
+ XI=XIV[i]
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ L[i]=-optim(par=1,fn=profilelikelihood)$value }
+ plot(XIV,L,type="l")
+ PL=function(XI){
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ return(optim(par=1,fn=profilelikelihood)$value)}
+ (L0=(OPT=optimize(f=PL,interval=c(0,10)))$objective)
+ profilelikelihood=function(beta){
+ -loglikelihood(1,beta) }
+ (L1=optim(par=1,fn=profilelikelihood)$value)
+ LR[j]=L1-L0
+ P[j]=1-pchisq(L1-L0,df=1)
+ G=gpd(X,u)
+ ES[j]=G$par.ests[1]
+ SES[j]=G$par.ses[1]
+ }
>
> plot(U,LR,type="b",ylim=range(c(0,LR)))
> abline(h=qchisq(.95,1),lty=2)
with on top the values of the ratio (the dotted line is the quantile of a chi-square distribution with one degree of freedom) and below the associated p-value
> plot(U,P,type="b",ylim=range(c(0,P))) > abline(h=.05,lty=2)
In order to compare, it is also possible to look at confidence interval for the tail index of the GPD fit,
> plot(U,ES,type="b",ylim=c(0,1)) > lines(U,ES+1.96*SES,type="h",col="red") > abline(h=1,lty=2)
To go further, see Falk (1995), Dietrich, de Haan & Hüsler (2002), Hüsler & Li (2006) with the following table, or Neves & Fraga Alves (2008). See also here or there (for the latex based version) for an old paper I wrote on that topic.
This week, we conclude the part on extremes with an application of extreme value theory to risk measures. We have seen last week that, if we assume that above a threshold 
, a Generalized Pareto Distribution will fit nicely, then we can use it to derive an estimator of the quantile function (for percentages such that the quantile is larger than the threshold)
It the threshold is 
, i.e. we keep the 
largest observations to fit a GPD, then this estimator can be written
The code we wrote last week was the following (here based on log-returns of the SP500 index, and we focus on large losses, i.e. large values of the opposite of log returns, plotted below)
> library(tseries)
> X=get.hist.quote("^GSPC")
> T=time(X)
> D=as.POSIXlt(T)
> Y=X$Close
> R=diff(log(Y))
> D=D[-1]
> X=-R
> plot(D,X)
> library(evir)
> GPD=gpd(X,quantile(X,.975))
> xi=GPD$par.ests[1]
> beta=GPD$par.ests[2]
> u=GPD$threshold
> QpGPD=function(p){
+ u+beta/xi*((100/2.5*(1-p))^(-xi)-1)
+ }
> QpGPD(1-1/250)
97.5%
0.04557386
> QpGPD(1-1/2500)
97.5%
0.08925095
This is similar with the following outputs, with the return period of a yearly event (one observation out of 250 trading days)
> gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/250, ci.type = + "likelihood", ci.p = 0.95,like.num = 50) Lower CI Estimate Upper CI 0.04172534 0.04557655 0.05086785
or the decennial one
> gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/2500, ci.type = + "likelihood", ci.p = 0.95,like.num = 50) Lower CI Estimate Upper CI 0.07165395 0.08925558 0.13636620
Note that it is also possible to derive an estimator for another population risk measure (the quantile is simply the so-called Value-at-Risk), the expected shortfall (or Tail Value-at-Risk), i.e.
The idea is to write that expression
so that we recognize the mean excess function (discussed earlier). Thus, assuming again that above (and therefore above that high quantile) a GPD will fit, we can write
or equivalently
If we substitute estimators to unknown quantities on that expression, we get
The code is here
> EpGPD=function(p){
+ u-beta/xi+beta/xi/(1-xi)*(100/2.5*(1-p))^(-xi)
+ }
> EpGPD(1-1/250)
97.5%
0.06426508
> EpGPD(1-1/2500)
97.5%
0.1215077
An alternative is to use Hill’s approach (used to derive Hill’s estimator). Assume here that 
, where 
is a slowly varying function. Then, for all 
,
Since is a slowly varying function, it seem natural to assume that this ratio is almost 1 (which is true asymptotically). Thus
i.e. if we invert that function, we derive an estimator for the quantile function
which can also be written
(which is close to the relation we derived using a GPD model). Here the code is
> k=trunc(length(X)*.025)
> Xs=rev(sort(as.numeric(X)))
> xiHill=mean(log(Xs[1:k]))-log(Xs[k+1])
> u=Xs[k+1]
> QpHill=function(p){
+ u+u*((100/2.5*(1-p))^(-xiHill)-1)
+ }
with the following Hill plot
For yearly and decennial events, we have here
> QpHill(1-1/250) [1] 0.04580548 > QpHill(1-1/2500) [1] 0.1010204
Those quantities seem consistent since they are quite close, but they are different compared with empirical quantiles,
> quantile(X,1-1/250) 99.6% 0.04743929 > quantile(X,1-1/2500) 99.96% 0.09054039
Note that it is also possible to use some functions to derive estimators of those quantities,
> riskmeasures(gpd(X,quantile(X,.975)),1-1/250) p quantile sfall [1,] 0.996 0.04557655 0.06426859 > riskmeasures(gpd(X,quantile(X,.975)),1-1/2500) p quantile sfall [1,] 0.9996 0.08925558 0.1215137
(in this application, we have assumed that log-returns were independent and identically distributed… which might be a rather strong assumption).
These data were collected at Copenhagen Reinsurance and comprise 2167 fire losses over the period 1980 to 1990, They have been adjusted for inflation to reflect 1985 values and are expressed in millions of Danish Kron. Note that it is possible to work with the same data as above but the total claim has been divided into a building loss, a loss of contents and a loss of profits.
> base1=read.table( + "http://freakonometrics.free.fr/danish-univariate.txt", + header=TRUE) > base2=read.table( + "http://freakonometrics.free.fr/danish-multivariate.txt", + header=TRUE)
Consider here the first dataset (we deal – so far – with univariate extremes),
> X=base1$Loss.in.DKM > D=as.Date(as.character(base1$Date),"%m/%d/%Y") > plot(D,X,type="h")
The graph is the following,
A natural idea is then to plot
i.e.
> Xs=sort(X) > logXs=rev(log(Xs)) > n=length(X) > plot(log(Xs),log((n:1)/(n+1)))
Points are on a straight line here. The slope can be obtained using a linear regression,
> B=data.frame(X=log(Xs),Y=log((n:1)/(n+1))) > reg=lm(Y~X,data=B) > summary(reg) Call: lm(formula = Y ~ X, data = B) Residuals: Min 1Q Median 3Q Max -0.59999 -0.00777 0.00878 0.02461 0.20309 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.089442 0.001572 56.88 <2e-16 *** X -1.382181 0.001477 -935.55 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.04928 on 2165 degrees of freedom Multiple R-squared: 0.9975, Adjusted R-squared: 0.9975 F-statistic: 8.753e+05 on 1 and 2165 DF, p-value: < 2.2e-16 > reg=lm(Y~X,data=B[(n-500):n,]) > summary(reg) Call: lm(formula = Y ~ X, data = B[(n - 500):n, ]) Residuals: Min 1Q Median 3Q Max -0.48502 -0.02148 -0.00900 0.01626 0.35798 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.186188 0.010033 18.56 <2e-16 *** X -1.432767 0.005105 -280.68 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.07751 on 499 degrees of freedom Multiple R-squared: 0.9937, Adjusted R-squared: 0.9937 F-statistic: 7.878e+04 on 1 and 499 DF, p-value: < 2.2e-16 > reg=lm(Y~X,data=B[(n-100):n,]) > summary(reg) Call: lm(formula = Y ~ X, data = B[(n - 100):n, ]) Residuals: Min 1Q Median 3Q Max -0.33396 -0.03743 0.02279 0.04754 0.62946 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.67377 0.06777 9.942 <2e-16 *** X -1.58536 0.02240 -70.772 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1299 on 99 degrees of freedom Multiple R-squared: 0.9806, Adjusted R-squared: 0.9804 F-statistic: 5009 on 1 and 99 DF, p-value: < 2.2e-16
The slope here is somehow related to the tail index of the distribution. Consider some heavy tailed distribution, i.e. 
, so that 
, where 
is some slowly varying function. Equivalently, the exists a slowly varying function 
such that 
. Then
i.e. since a natural estimator for 
is the order statistic 
, the slope of the straight line is the opposite of tail index 
. The estimator of the slope is (considering only the 
largest observations)
Hill‘s estimator is based on the assumption that the denominator above is almost 1 (which means that 
, as 
), i.e.
Note that, if 
, but not two fast, i.e. as , then 
(one can even get 
with stronger convergence assumptions). Further
Based on that (asymptotic) distribution, it is possible to get a (asymptotic) confidence interval for
> xi=1/(1:n)*cumsum(logXs)-logXs > xise=1.96/sqrt(1:n)*xi > plot(1:n,xi,type="l",ylim=range(c(xi+xise,xi-xise)), + xlab="",ylab="",) > polygon(c(1:n,n:1),c(xi+xise,rev(xi-xise)), + border=NA,col="lightblue") > lines(1:n,xi+xise,col="red",lwd=1.5) > lines(1:n,xi-xise,col="red",lwd=1.5) > lines(1:n,xi,lwd=1.5) > abline(h=0,col="grey")
It is also possible to work with 
, then 
. And similarly 
as (and again 
with additional assumptions on the rate of convergence), and
(obtained using the delta-method). Again, we can use that result to derive (asymptotic) confidence intervals
> alpha=1/xi > alphase=1.96/sqrt(1:n)/xi > YL=c(0,3) > plot(1:n,alpha,type="l",ylim=YL,xlab="",ylab="",) > polygon(c(1:n,n:1),c(alpha+alphase,rev(alpha-alphase)), + border=NA,col="lightblue") > lines(1:n,alpha+alphase,col="red",lwd=1.5) > lines(1:n,alpha-alphase,col="red",lwd=1.5) > lines(1:n,alpha,lwd=1.5) > abline(h=0,col="grey")
The Deckers-Einmahl-de Haan estimator is
where for
Then (given again conditions on the speed of convergence i.e. , with as ),
Finally, Pickands‘ estimator
it is possible to prove that, as ,
Here the code is
> Xs=rev(sort(X)) > xi=1/log(2)*log( (Xs[seq(1,length=trunc(n/4),by=1)]- + Xs[seq(2,length=trunc(n/4),by=2)])/ + (Xs[seq(2,length=trunc(n/4),by=2)]-Xs[seq(4, + length=trunc(n/4),by=4)]) ) > xise=1.96/sqrt(seq(1,length=trunc(n/4),by=1))* +sqrt( xi^2*(2^(xi+1)+1)/((2*(2^xi-1)*log(2))^2)) > plot(seq(1,length=trunc(n/4),by=1),xi,type="l", + ylim=c(0,3),xlab="",ylab="",) > polygon(c(seq(1,length=trunc(n/4),by=1),rev(seq(1, + length=trunc(n/4),by=1))),c(xi+xise,rev(xi-xise)), + border=NA,col="lightblue") > lines(seq(1,length=trunc(n/4),by=1), + xi+xise,col="red",lwd=1.5) > lines(seq(1,length=trunc(n/4),by=1), + xi-xise,col="red",lwd=1.5) > lines(seq(1,length=trunc(n/4),by=1),xi,lwd=1.5) > abline(h=0,col="grey")
It is also possible to use maximum likelihood techniques to fit a GPD distribution over a high threshold.
> library(evd) > library(evir) > gpd(X,5) $n [1] 2167 $threshold [1] 5 $p.less.thresh [1] 0.8827873 $n.exceed [1] 254 $method [1] "ml" $par.ests xi beta 0.6320499 3.8074817 $par.ses xi beta 0.1117143 0.4637270 $varcov [,1] [,2] [1,] 0.01248007 -0.03203283 [2,] -0.03203283 0.21504269 $information [1] "observed" $converged [1] 0 $nllh.final [1] 754.1115 attr(,"class") [1] "gpd"
or equivalently (or almost)
> gpd.fit(X,5) $threshold [1] 5 $nexc [1] 254 $conv [1] 0 $nllh [1] 754.1115 $mle [1] 3.8078632 0.6315749 $rate [1] 0.1172127 $se [1] 0.4636270 0.1116136
The interest of the latest function is that it is possible to visualize the profile likelihood of the tail index,
> gpd.profxi(gpd.fit(X,5),xlow=0,xup=3)
or
> gpd.profxi(gpd.fit(X,20),xlow=0,xup=3)
Hence, it is possible to plot the maximum likelihood estimator of the tail index, as a function of the threshold (including a confidence interval),
> GPDE=Vectorize(function(u){gpd(X,u)$par.ests[1]}) > GPDS=Vectorize(function(u){ + gpd(X,u)$par.ses[1]}) > u=c(seq(2,10,by=.5),seq(11,25)) > XI=GPDE(u) > XIS=GPDS(u) > plot(u,XI,ylim=c(0,2)) > segments(u,XI-1.96*XIS,u,XI+ + 1.96*XIS,lwd=2,col="red")
Finally, it is possible to use block-maxima techniques.
> gev.fit(X) $conv [1] 0 $nllh [1] 3392.418 $mle [1] 1.4833484 0.5930190 0.9168128 $se [1] 0.01507776 0.01866719 0.03035380
The estimator of the tail index is here the last coefficient, on the right.
Since it is rather difficult to install a package in class rooms, here is the source of rcodes used here (to fit a GPD for exceedances)
> source("http://freakonometrics.blog.free.fr/public/code/gpd.R")
Next time, we will discuss how to use those estimators.