Tail index estimation

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",
+ "http://freakonometrics.free.fr/danish-multivariate.txt",

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.

Quel lancer de dé faire à Serpents et Échelles ?

Hier soir, j’évoquais l’utilisation des chaînes de Markov au jeu Serpents et Échelles. Comme me le faisait remarquer Jean-Philippe, de manière assez étrange, les enfants sont toujours beaucoup plus content après avoir fait un 6 qu’après avoir fait un 1. Mais est-ce réellement la valeur optimale que le dé doit prendre ? Ça dépend bien sûr de la position. Par exemple, au premier lancer, on peut se demander ce que deviendrait le nombre de tours moyen nécessaires pour finir le jeu, conditionnellement aux 6 lancers possibles. On peut calculer, toujours conditionnellement aux 6 valeurs possibles de dés (et des positions où on va se retrouver), l’espérance du nombre de tours pour finir la partie,

```esperance=function(h0){
INITIAL = as.numeric(which(M[h0+1,]>0))-1
ESPERANCE=rep(NA,length(INITIAL))
names(ESPERANCE)=INITIAL
for(k in 1:length(INITIAL)){
initial=rep(0,n+1); initial[INITIAL[k]]=1
distrib=initial%*%M
game=rep(NA,1000)
for(h in 1:length(game)){
game[h]=distrib[n+1]
distrib=distrib%*%M}
ESPERANCE[k]=sum(1-game)}
return(ESPERANCE)}```

(à partir du code mis en ligne hier), e.g. pour le premier lancer,

```> esperance(0)
1        2        3        5        6       14
32.16499 31.99947 31.82348 31.47954 31.36441 29.83020```

où la valeur indiquée est la position où l’on pourrait se retrouver (la dernière correspond à un 4). On note que le meilleur premier lancer possible est celui qui nous amène sur la première échelle, i.e. la 4ème case,

Si on regarde case par case, on a les valeurs des “meilleurs” lancers de dés suivant (sans forcément l’unicité) en fonction de sa position sur le plateau

(avec en rouge les serpents et en bleu les échelles). Amusant non ?

Minuit, l’heure du crime ?

Il y a quelques jours, j’étais invité par Rémi pour me rendre au poste de police pour discuter de la possibilité de consulter quelques données. Et Rémi m’a donné un historique de près d’un million d’appels téléphonique d’urgence, sur Montréal, passés au cours des six dernières années. Histoire de prendre un peu en main les données, j’ai voulu regarder (inspiré par une remarque de Rémi) si minuit était vraiment une heure de crime1

Ce dont je dispose ce sont les heures d’appels d’urgence (le fameux 911), avec un code (formellement, le code indiqué par la personne lors de l’appel, pas celui établi par la patrouille une fois sur place). Globalement, la distribution des appels (par seconde) au cours de la journée ressemble à ça,

avec les six années en noir, 2005 en bleu et 2010 en rouge (histoire de visualiser une stabilité temporelle).
On peut aussi se focaliser sur les causes données lors de l’appel, afin de voir quel criminel commet un crime à minuit. Par exemple, pour les alarmes consécutives à un cambriolage, on obtient la répartition suivante,

Manifestement, minuit est un creux en matière de cambriolage. N’en déplaise à Arsène Lupin. Les hold-ups, eux aussi surviennent pendant la journée,

Les intrusions-effractions aussi

Sinon, de manière générale, les infractions au code criminel surviennent en fin de journée (mais baissent en fin de soirée)

et les conflits, suivent la même distribution

Autrement dit, les conflits surviennent quand on rentre du travail… mais pas à minuit. En revanche, si on commence à regarder les appels pour signaler descoups de feux, ils surviennent un peu après les conflits… mais un peu avant minuit.

Ah, mais finalement je crois que j’ai trouvé le crime qui survient à minuit,

Et dans ma base, cela correspond au code 063, ce qui signifie trouble à la paix. Autrement dit, si un appel arrive au 911 à minuit, c’est rarement pour signaler un (vrai) crime (au sens de ceux qu’on trouve dans les polars), mais plus pour signaler un problème de beuverie (que l’on peut aussi trouver chez Ian Rankin d’ailleurs). Damned ! Moi qui pensait résoudre des crimes avec mes séries de chiffres…
.
1 Je n’ai pas réussi à trouver d’où pouvait sortir cette idée reçue. En lisant le dernier Ian Rankin (ou plutôt la dernière enquête de John Rebus), je me suis fait la réflexion que dans beaucoup d’enquêtes, les policiers arrivaient vers minuit. J’ai eu l’impression que c’était aussi souvent le cas dans bon nombre d’épisodes de séries télé “policières” (CSI, Castle, Bones, ). En fait si quelqu’un sait où je pourrais trouver une base avec des dates de décès des personnages dans les séries policières ou des polars, je suis preneur.

Wald, score et rapport de vraisemblance

Vendredi, nous devrions aborder un peu la problématique des tests asymptotiques. Avant, rappelons que l’estimateur du maximum de vraisemblance de  vérifie la propriété asymptotique suivante.

Or on a une propriété intéressante sur les convergences de vecteurs Gaussiens (pour rester assez général). Soit  une suite de vecteurs aléatoires telle que  où . Alors

où  (pour un paramètre univarié, ça sera alors ). Mais ce résultat n’est intéressant que si  est connue. En fait, si on a la suite des variances vérifie , alors

Traduit sur notre estimateur du maximum de vraisemblance, cela signifie que

ou encore

Trois tests peut alors être mis en œuvre (et tant qu’à faire, on peut essayer de l’illustrer sur un exemple: encore et toujours du pile/face).

```> set.seed(1)
> n=20
> X=sample(0:1,size=n,replace=TRUE)
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
> LL=sapply(p,logL)
> plot(p,LL,type="l",col="red",lwd=2)
> p0=.5
> points(p0,logL(p0),pch=3,cex=1.5,lwd=2)
> abline(v=p0,lty=2)```

On a alors 20 tirages de pile ou face, on a obtenu 11 piles, et on se demande si la pièce est “juste“.Mais avant toute chose, commençons par calculer l’estimateur du maximum de vraisemblance, ainsi que le score et l’information de Fisher. Dans ce modèle de variables de Bernoulli, on connaît des formes explicites de ces quantités. Mais ici, on va utiliser la méthode la plus générale qui soit, i.e. on va maximiser la log-vraisemblance, puis on va calculer le score et l’information de Fisher à la valeur de  que l’on souhaite tester.

```> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=p0,method="BFGS")\$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> nx=sum(X==1)
> f = expression(nx*log(p)+(n-nx)*log(1-p))
> Df = D(f, "p")
> Df2 = D(Df, "p")
> p=p0
> score=eval(Df)
> (IF=-eval(Df2))
[1] 80
> 1/(p0*(1-p0)/n)
[1] 80```

On note que l’information de Fisher calculée par double différenciation de la log-vraisemblance est identique à la formule analytique. Pour rappels, on a

Trois tests (équivalent asymptotiquement comme on le verra en cours) sont alors possibles.

Tout d’abord le test de Wald propose de travailler sur la différence entre l’estimateur du maximum de vraisemblance, et la valeur que l’on cherche à tester (comme le montre le dessin ci-dessous). Cette différence doit être “petite” si  est la vraie valeur.On peut alors utiliser la statistique suivante

Asymptotiquement, en utilisant le résultat évoqué au début, cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

```> alpha=0.05
> (T=(pml-p0)^2*IF)
[1] 0.1999970
> T<qchisq(1-alpha,df=1)
[1] TRUE```

Ici, on accepte l’hypothèse que la pièce n’est pas pipée.

Ensuite le test du rapport de vraisemblance propose de travailler sur les valeurs de la log-vraisemblance. Là encore, cette différence doit être “petite” si  est la vraie valeur. Graphiquement, on mesure une distance normalisée

La statistique est celle d’un test bilatéral,

On pose alors

(le 2 sera détaillé un peu en cours, c’est mis là de manière à avoir trois tests équivalents). Asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

```> (T=2*(logL(pml)-logL(p0)))
[1] 0.2003347
> T<qchisq(1-alpha,df=1)
[1] TRUE```

Là encore, on accepte l’hypothèse que la pièce n’est pas pipée.

Enfin, le test du score propose de travailler sur la pente de la log-vraisemblance en  . Cette pente doit être “petite” si  est la vraie valeur.

La pente correspond au score,

La statistique de test est alors ici

Et là encore, asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

```> (T=slope^2/IF)
[1] 0.2
> T<qchisq(1-alpha,df=1)
[1] TRUE```

Et le test accepte ici encore l’hypothèse que la pièce n’est pas pipée.

Extracting information from a keyboard…

Yesterday, Baptiste published a post on “ethno-photography” (here). As he mentioned it, in Paris 8, they experience a real absence of serious cleaning of office equipment. He then shows the keyboard of the only computer they can use in the sociology department (for forty researchers),

Apart from the fact that everyone in France should be ashamed to see how much is spent in universities (which is the first information we have from that picture), we should also be able to guess in which langage people work in this department.
I considered three books (two in French, one in English) and I would like to see the frequency of each letter,

• Mauss, manuel d’éthnographie (here), 1926
• Durkheim, Livre II: Les croyances élémentaires in Les formes élémentaires de la vie religieuse (here), 1912
• Ferri, Criminal Sociology (here), 1896

Those three books are in rich text format, I just changed it to get text files… Then, it is easy to count appearance of letters. E.g. for Mauss,

```> library(corpora)
> textfile=scan("MAUSS-manuel.txt",
+ what="char", sep="\n")
> textfile<-tolower(textfile)
> M=NA
> for(i in 1:length(textfile)){
+ line=textfile[i]
+ M=c(M,strsplit(as.character(line),"")[[1]])
+ }
> T=table(M)
> T
M
'     -           \t     !     "     %     &     (     )     ,     .     /
5308  1049 86589    44     3     3     2     2   370   391  6609  4909    12
:     ;     ?     @     ]     _     ~     ’     =     «     »     ¬     °
819  1178   113     1     1     4     1    39     1   108   107   823     3
…     0     1     2     3     4     5     6     7     8     9     a     à
1    69   213    83    73    34    48    33    28    64   151 30559  1651
â     ä     b     c     ç     d     e     é     è     ê     ë     f     g
224     3  3562 14678   110 17713 63955 10354  1798  1000     5  4555  4911
h     i     î     ï     j     k     l     m     n     ñ     o     º     ô
4359 30851   226    47  1147   247 24792 12844 32525     6 25562     2   151
ö     œ     p     q     r     s     t     u     ù     û     ü     v     w
12    52 12696  4667 28237 37630 32945 25001   211    40     9  4787   164
x     y     z
1996  1222   343```

Then, we can summarize in to see proportion of standard 26 letters, and we have, for Mauss,

and for Durkeim,

If we compare the two, we have almost the same proportions,

If we look at our book in English now, we have

i.e., if we compare with Mauss for instance

So we have much more E in French than in English, but still, people writing in English use a lot the E. So looking at the E should not give us any clue…. But we can see that in English, the H is as common as the L, or the C. Not in French, where L is much more frequent than the H. But on the picture, the C is more clear than the H. We can also look at the U, which is common in French, not in English… Here, on the keyboard, it is perfectly clear… so I guess people use it frequently.
So I would say that they write more in French than the write in English, on that computer.
Actually, the same idea has been used a long time ago on calculators to see that Benford’s law works: some numbers are really used (as well as the legend pretends that some pages in logarithm books were never used….), see here orthere. So Baptiste, if one day the keyboard is cleaned up, please send me another picture after a few weeks to see if things have changed….
An for those who cannot imagine how it is to work in some universities in France, just look at his blog (here). Pictures are unbelievable….Good luck Baptiste….

Qui a eu l’idée d’élargir les trottoirs en France ?

J’étais tombé l’autre jour sur un joli dessin expliquant pourquoi les policiers et les syndicats n’estimaient jamais le même nombre de manifestants, avec une explication scientifique, et vaguement rationnelle, ici.

C’est  effectivement joli… et on retrouve un rapport de 3 dont pas mal de monde semble parler, sauf que si c’était vrai, le rapport entre les estimations de la police et des syndicats devrait être stable d’une manifestation à l’autre, non ?
Si on compare le 19 mars dernier (ici) et le 16 octobre (), on arrive aux chiffres suivants,

Les traits pointillés sont la droite de régression. La courbe en trait plein et la régression (passant par l’origine car on cherche un rapport de proportionnalité) en prenant en compte un effet taille, i.e. en introduisant des poids différents dans ma régression. J’ai pris la taille des villes, autrement dit, plus la ville est grosse, plus je donne un poids important au rapport… L’idée étant d’intégrer un effet médiatique, même si la pente n’est plus le ratio moyen…  Avec cette dernière méthode, j’obtiens des pentes de 4,06 pour mars, et 6,12 pour octobre… ce qui est à un peu différents des rapports nationaux mentionnés sur le site du monde, i.e. (3 millions selon les syndicats et 1,2 selon la police, en mars, contre  2,5 millions contre 825 milliers l’autre jour). Notons qu’on peut vérifier rapidement que ce n’est pas un problème statistique, les valeurs sont sensiblement différentes, si l’on compare les distributions des estimateurs..

Bref, quelque chose semble avoir changé entre ces deux dates…. la taille des trottoirs ?

Margin of error, and comparing proportions in the same sample

Irecently tried to answer a simple question, asked by @adelaigue. Actually, I thought that the answer would be obvious… but it is a little bit more compexe than what I thought. In a recent survey about elections in Brazil, it was mentionned in a French newspapper that “Mme Rousseff, 62 ans, de 46,8% des intentions de vote et José Serra, 68 ans, de 42,7%” (i.e. proportions obtained from the survey). It is also mentioned that “la marge d’erreur du sondage est de 2,2% ” i.e. the margin of error is 2.2%, which means (for the journalist) that there is a “grande probabilité que les 2 candidats soient à égalité” (there is a “large probability” to have equal proportions).
Usually, in sampling theory, we look at the margin of error of a single proportion. The idea is that the variance of $\widehat{p}$, obtained from a sample of size  is

thus, the standard error is

The standard 95% confidence interval, derived from a Gaussian approximation of the binomial distribution is

The largest value is obtained when p is 1/2, and then we have a worst case confidence interval (an upper bound) which is

So with a margin of error  means that . Hence, with a 5% margin of error, it means that n=400. While 2.2% means that n=2000:
> 1/.022^2
[1] 2066.116
Classically, we compare proportions between two samples: surveys at two different dates, surveys in different regions, surveys paid by two different newpapers, etc. But here, we wish to compare proportions within the same sample. This has been consider in an “old” paper published in 1993 in the American Statistician,

It contains nice figures to illustrate the difference between the standard approach,

and the one we would like to study here.

This point is mentioned in the book by Kish, survey sampling (thanks Benoit for the reference),

Let and denote empirical frequencies we have obtained from the sample, based on  observations. Then since

and

we have

Thus, a natural margin of error on the difference between the two proportion is here

which is here 4 points
> n=2000
> p1=46.8/100
> p2=42.7/100
> 1.96*sqrt((p1+p2)-(p1-p2)^2)/sqrt(n)
[1] 0.04142327
Which is exactly the difference we have here ! Hence, the probability of reaching such a value is quite small (2%)
> s=sqrt(p1*(1-p1)/n+p2*(1-p2)/n+2*p1*p2/n)
> (p1-p2)/s
[1] 1.939972
> 1-pnorm(p1-p2,mean=0,sd=sqrt((p1+p2)-(p1-p2)^2)/sqrt(n))
[1] 0.02619152

Actually, we can compare the three margin of errors we have so far,

• the upper bound
• the “average one”

where

• the more accurate one we just obtained,

where .
> p=seq(0,.5,by=.01)
> ic1=rep(1.96/sqrt(4*n),length(p))
> ic2=1.96*sqrt(p*(1-p))/sqrt(n)
> delta=.01
> ic31=1.96*sqrt(2*p-delta^2)/sqrt(n)
> delta=.2
> ic32=1.96*sqrt(2*p-delta^2)/sqrt(n)
> plot(p,ic32,type=”l”,col=”blue”)
> lines(p,ic31,col=”red”)
> lines(p,ic2)
> lines(p,ic1,lty=2)
So on the graph below, the dotted line is the standard upper bound, the plain line in black being a more accurate one when the probability is  (the x-axis). The red line is the true margin of error with a large difference between candidates (20 points) and the blue line with a small difference (1 point).

Remark: an alternative is to consider a chi-square test, comparering two multinomial distributions, with probabilities  and  where is the average proportion, i.e. 44.75%. Then

i.e.  =3.71
> p=(p1+p2)/2
> (x2=n*((p1-p)^2/p+(p2-p)^2/p))
[1] 3.756425
> 1-pchisq(x2,df=1)
[1] 0.05260495
Under the null hypothesis, should have a chi-square distribution, with one degree of freedom (since the average is fixed here). Here the probability to reach that level is around 5% (which can be compared with the 2% we add before).

So finally, I would think that here, stating that there is a “large probability” is not correct…

Too large datasets for regression ? What about subsampling….

recently, a classmate working in an insurance company told me he had too large datasets to run simple regressions (GLM, which involves optimization issues), and that they were thinking of a reward for the one who will write the best R-code (at least the fastest). My first idea was to use subsampling techniques, saying that 10 regressions on 100,000 observations can take less time than a regression on 1,000,000 observations. And perhaps provide also better results…

• Time to run a regression, as a function of the number of observations

Here, I generate a dataset as follows

and we fit

where is a spline function (just to make it as general as possible, since in insurance ratemaking, we include continuous variates that do not influence claims frequency linearly in the score). Yes, there might be also useless variables, including one of them which is strongly correlated with one that has an impact in the regression. The code to generate the dataset is simply

```> n=10000
> X1=rexp(n)
> X2=sample(c("A","B","C"),size=n,replace=TRUE)
> X3=runif(n)
> Z=rmnorm(n,c(0,0),matrix(c(1,0.8,.8,1),2,2))
> X4=Z[,1]
> X5=Z[,2]
> X6=X1^2
> E=runif(n)
> lambda=.2*X5-4*dbeta(X3,2,5)+X1+
+1*(X2=="A")-2*(X2=="B")-5*(X2=="C")
> Y=rpois(n,exp(lambda))
> base=data.frame(Y,X1,X2,X3,X4,X5,X6,E)```

We would like the study the time it takes to run a regression, as a function of the size (i.e. the number of lines ) of the dataset.

```> system.time( glm(Y~bs(X1)+X2+X3+X4+
+ X5+X6+offset(log(E)),family=poisson,
+ data=base) )
utilisateur     système      écoulé
0.25        0.00        0.25```

Here, the time I look at is the last one. But so far, it was rather simple, but it is not the best model I can get. Let us use a stepwise (backward) variable selection,

```> system.time( step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=poisson,
+ data=base)) )
Start:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))
Step:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + offset(log(E))
Df Deviance    AIC
<none>        2236.0 2882.1
- X5      1   2240.1 2884.2
- X4      1   2244.1 2888.2
- X3      1   4783.2 5427.3
- X2      2   5311.4 5953.5
- bs(X1)  3   6273.7 6913.8
utilisateur     système      écoulé
1.82        0.03        1.86```

Finally, from the first regression, we have points in black (based on 200 simulated datasets), and with a stepwise procedure, we have the points in red.

i.e. it might look linear (proportional), but if it was linear, then on a log-log scale, we should have also straigh lines, with slope 1,

Actually, it looks like a convex function.

The interpretation of that convexity might lead to misinterpretation. On the graph below on the left, on a dataset two times bigger than the previous one (black point) will be less than two times longer to run, while on the right, it will be more than two timess longer,

Convexity can simply be interpreted as “too large datasets take time, and too small too…”. Which is a first step: it should be interesting, in some cases, to run several regressions on smaller datasets….

• Running 100 regressions on 100 lines, or running 1 regression on 10,000 lines ?

Here, we have datasets with =200,000 lines. The questions is how long will it take if we subdived into subsamples (of equal size), and run regressions ?

```> nk=trunc(n/k)rep(1:k,each=nk); nt=nk*k
> base=data.frame(Y[1:nt],X1[1:nt],
+ X2[1:nt],X3[1:nt],X4[1:nt],X5[1:nt],
+ X6[1:nt],E[1:nt],classe)
> system.time( for(j in 1:k){
+  glm(Y~bs(X1)+X2+X3+X4+X5+
+ X6+offset(log(E)),family=poisson
+ ,data=base,subset=classe==j) })
utilisateur     système      écoulé
1.31        0.00        1.31
> system.time( for(j in 1:k){
+      step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=
+ poisson,data=base,subset=classe==j)) })
Start:  AIC=183.97
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))```

[…]

```  Df Deviance    AIC
<none>        117.15 213.04
- X2      2   250.15 342.04
- X3      1   251.00 344.89
- X4      1   420.63 514.53
- bs(X1)  3   626.84 716.74
utilisateur     système      écoulé
11.97        0.03       12.31```

On the graph below, we have the time (y-axis, here on a log scale) it took to run regression on samples of size , as function of (x-axis), including the time it took to run the regression on a dataset of size which is the concentration of dots on the left (i.e. =1), both on the 6 regressors – in black – and with a strepwise procedure – in red. One has to keep in mind that I did not remove the printing option in the stepwise procedure, so it might be difficult to compare the two clouds (black vs. red). Nevertheless, we clearly see that if we run regression on samples of size , when is not too large, i.e. less than 10 or 15, it is not longer than the regression on =200,000 lines.

So here we see that running 100 regressions on 2,000 lines is longer than running 1 regression on 200,000 lines… But maybe we are not comparing things that are actually comparable: what if it takes a bit longer, but we strongely improve the quality of our estimators ?

• What about the quality of the output ?

Here, we consider only one dataset, with =100,000 lines (just to make it run a bit faster). And =20 subsets. Recall that the generated dataset is from

and we fit

Here, we plot here and a confidence interval, defined as

The lightblue segment is the initial estimator, while the blue one is obtained from the stepwise procedure. The grey area represent the estimation on the overall sample, while the segments on the right are the estimators (each on samples of size ).

We can see that we have much more volatility on those estimators, but the average (horizontal doted lines) are not so bad… The true value (i.e. the one used to generate the dataset is the dotter black horizontal line).
And if we repeat that on 1,000 simulated dataset, we obtaind the following distribution for (blue line), so we have an unbiased estimator of our parameter (the verticular line being here the true value), here including a stepwise procedure,

But if we add the the red curve is the average of the the previous one being now the clear blue line in the back, we see that taking average of estimators on subsamples is not bad at all, on the contrary,

and for those who think that the stepwise procedure is a mistake, here is what we get without it,

So what we can see is that running 20 regressions can take (a little) more time (from what we’ve seen earlier) than running only one on the whole dataset…. but it provides better estimates. So the tradeoff is not that simple, and maybe running several regressions on huge datasets can be a proper alternative.

Tails of copulas, une lecture graphique

Suite à une formation que je faisais en fin de semaine à Brest (les slides sont ici et ), je voulais revenir sur les histoires de tails of copulas, pour reprendre le titre de l’article (ici) de Gary Venter (et qui correspond à des choses que j’avais pu présenter il y a quelques années à Berlin, les slides étant en ligne ici).

• Quantifier la dépendance de queue

L’idée est de noter qu’il est noter qu’il existe deux manières de quantifier la dépendance de queue. La première est liée à l’approche de Joe (1990, ici, ou 1997 pour le livre), qui a introduit un (strong) tail dependence index. Par exemple pour la queue inférieure,

soit

La seconde est liée à une idée que l’on retrouve dans les travaux de Janet Heffernan, Stuart Coles ou Jonathan Tawn. L’intuition est la suivante (on peut la retrouver en ligne ici). Si  et  ont la même loi et que l’on suppose les variables indépendantes, alors

En revanche, si les variables sont comonotones (c’est à dire égales comme on suppose les lois identiques),

Aussi, on peut supposer qu’il existe un indice tel que

Le soucis est que le cas d’indépendance correspond à =2, alors que le cas de dépendance forte correspond au cas =1. Il est alors usuel de faire une transformation affine pour se ramener sur [0,1], et que la force de la dépendance soit croissante avec , e.g.

Posons alors

qui pourra être interprété comme un (weak) tail dependence index.
Bref, ces deux mesures donnent de l’information sur le comportement dans les queues de distribution.

• Les fonctions de concentration dans les queues

L’idée est de noter qu’il est possible d’étudier ces fonctions afin de mieux comprendre le comportement dans les queues. En s’inspirant de Gary Venter, on peut définir

pour étudier le comportement dans la queue inférieure, et

pour la queue supérieure,où  est la copule de survie associée à , au sens où

et

Cet outil permettra de modéliser la dépendance forte. On peut également poser, afin d’étudier la dépendance faible,

ou

• Application statistique

L’idée est de noter qu’il est facile d’estimer ces fonctions. Ces outils peuvent être utiles pour mieux comprendre le comportement dans les queues.
Par exemple pour une copule Gaussienne de corrélation 0,5, on a la forme théorique suivante pour les fonctions de concentration (au sens fort)

Statistiquement, il est possible d’estimer ces quantités en comptant simplement le nombre d’observations dans le coin inférieur gauche, ou le coin supérieur droit.  Si on dispose d’un échantillon, on peut alors regarder ce que donnent les versions

et

Pour un échantillon de taille n=500, on obtient les intervalles de confiance à 90% de la forme suivante,

Le code R ressemble à ça

```> library(evd); data(lossalae)
> cor(lossalae,method="spearman")
Loss     ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000```

avec le code suivant pour la version empirique,

```> z=seq(0,.5,by=.001)
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+  Lemp[i]=sum((U<=z[i])&(V<=z[i]))/sum(U<=z[i])
+  Remp[i]=sum((U>=1-z[i])&(V>=1-z[i]))/sum(U<=z[i])
+ }```

et pour la version théorique,

```> Lg=(pcopula(copclayton,cbind(z,z)))/(z)
> Rg=((1-2*(1-z)+pcopula(copclayton,cbind(1-z,1-z))))/(z)
> plot(c(1-z,z),c(Lg,Rg))```

De plus, on a des fonctions similaires pour la dépendance au sens faible, avec le code suivant pour la version théorique,

```> Lg=log(pcopula(cop,cbind(z,z)))/log(z)
> Rg=log((1-2*(1-z)+pcopula(cop,cbind(1-z,1-z))))/log(z)
> Lg=1/Lg*2-1
> Rg=1/Rg*2-1```

et celui là pour la version empirique

```> z=seq(0,.5,by=.001)
> v <- lossalae
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+  Lemp[i]=log(mean((U<=z[i])&(V<=z[i])))/log(mean(U<=z[i]))
+  Remp[i]=log(mean((U>=1-z[i])&(V>=1-z[i])))/log(mean(U<=z[i]))
+ }
> Lemp=1/Lemp*2-1
> Remp=1/Remp*2-1```

Bref, on peut utiliser ces fonctions sur des vrais échantillons. Considérons l’exemple classique loss-alae (où l’on couple les frais dans des sinistres assurés, et les frais payés par l’assureur). On souhaite ajuster une copule, sans trop savoir laquelle. On peut commencer par étudier la dépendance forte, et comparer avec une copule Gaussienne. La copule Gaussienne de référence possède ici le même rho de Spearman que l’échantillon dont on dispose,

```> cor(lossalae,method="spearman")
Loss     ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000
> library(copula)
> paramgauss=.47
> paramclayton=.9
> paramgumbel=1.45
> copgauss=normalCopula(paramgauss)
> copclayton=claytonCopula(paramclayton, dim = 2)
> copgumbel=gumbelCopula(paramgumbel, dim = 2)```

On obtient ici

La courbe verte est l’intervalle de confiance (ponctuel) à 95% pour une copule Gaussienne et un échantillon de même taille. On voit qu’on modélise mal la structure de dépendance. Avec une copule duale de Clayton, on obtient

et enfin pour une copule de Gumbel,

Bref, la copule de Gumbel semble réellement bien adaptée… Si on creuse en étudiant la dépendance au sens faible, on peut valider là aussi ce modèle. En effet, si la référence est la copule Gaussienne,

ou pour une copule de Clayton,

alors qu’une copule de Gumbel donnerait

Le bagplot, comme extension du boxplot

Lundi 15, Michel Grun Rehomme avait présenté les travaux qu’il avait obtenu avec Noureddine Ben Lagha sur des problèmes d’écrêtement de coûts de sinistres, que l’on peut rapprocher de la détection de points atypiques, ou outliers. Ce matin, Alexander McNeil nous a rappelé que les bagplot pouvaient être très intéressants en gestion des risques (mentionné ici). Je ne peux pas m’empêcher de revenir un petit peu sur cette notion, en m’écartant des applications financières proposées par Alexander pour proposer des applications en démographie et en climatologie (inspiré par des travaux de Rob Hyndman).

• Le boxplot et la régression quantile

Un outil classique (en dimension 1) est le boxplot, proposé par John Tukey en 1977. Par exemple, sur les coûts de sinistres en assurance auto, on obtient

On peut aussi faire un boxplot par classes (ou tranches) d’âge,

qui correspond, en passant à la limite (i.e. en considérant l’âge comme une variable continue), à un ensemble de régressions quantiles.

On peut aussi regarder l’impact de l’ancienneté du véhicule (à la souscription du contrat) sur le coût individuel des sinistres,

avec en rouge la régression (linéaire) du log du coût sur l’âge du véhicule. Là aussi, on peut faire des régressions quantiles,

Le code pour obtenir les régressions quantiles lissées est assez simple, une fois installé les library(quantreg) et library(splines),

```> library(quantreg);  library(splines)
> fit = rq(log(cout)~bs(agecond,df=8),data=BASECOUT,tau=.5)
> xu=seq(10,110)
> yu=predict(fit,newdata=data.frame(agecond=xu))
> lines(xu,yu,col="red",lwd=2)```

Cette notion a été généralisée en multivariée sous le nom de bagplot. Par exemple, pour déterminer la “médiane” en dimension 2, on va recherche un ensemble (appelé bag) contenant 50% des observations.

• La profondeur

Cela se fait à l’aide du concept de profondeur statistique, introduite par John Tukey en 1975. En dimension 1, si

désigne la fonction de répartition, la “profondeur de demi-espace” s’écrit

qui est maximal pour la médiane (on obtient alors ). Plus généralement, pour un quantile d’ordre ,

Cette notion se généralise en dimension supérieure en posant

On notera que Fraiman et Muñiz ont proposé une extension en 2001. S’il n’y a pas d’ex-aequos, en dimension 1,

qui appartient à , et dont la médiane est la plus profonde. Liu a proposé une profondeur dite simpliciale en 1984,

D’un point de vue statistique, la fonction de répartition est inconnue, mais on dispose de . On note alors l’équivalent empirique de la profondeur. En dimension 2, ça sera  où . La méthode de construction est détaillée ci-dessous,

• Construction d’un bagplot

Pour cela, on note  le plus petit nombre d’observations contenu dans un demi plan passant par . On note  l’ensemble des points  tels que  contienne au moins  points. On notera qu’on peut montrer que , et que ces régions sont convexes. Rousseeuw et Ruts ont proposé un algorithme permettant de construire un tel ensemble de région en 1998. Sur l’exemple suivant, on considère 20 points en dimension 2,

Notons qu’en changeant de méthode, on obtient des résultats qui peuvent être très changeant.

```> library(mnormt)
> Z=rmnorm(20,mean=c(0,0),varcov=matrix(c(1,.6,.6,1),2,2))
> source("http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/bagplot/bagplot.R")
> bagplot(Z,factor=1,dkmethod=1)
> bagplot(Z,factor=1,dkmethod=2)```

La méthodologie a été proposée par Rousseeuw, Ruts et Tukey en 1999 (ici).

• Copules Gaussiennes, versus copules de Clayton

On peut essayer de voir si les bags changent de forme lorsque la copule sous-jacente change. L’exemple suivant montre le bagplot d’un vecteur gaussien,

alors que l’exemple suivant montre le bagplot pour une copule duale de Clayton, avec des marges Gaussiennes.

Dans le premier cas, on retrouve que les bagplots ressemblent des ellipses (ce qui correspond à ce que j’avais raconté ici sur les ellipses de confiance). On retrouve cette ellipse sur le graphique dessous, avec 5000 points (simulés) au lieu de 200,

pour la copule Gaussienne, et

pour la copule de Clayton. La petite animation ci-dessous montre l’évolution des bag en fonction de la corrélation du vecteur Gaussien,

• Illustration sur les données Loss-ALAE

Considérons la base (célèbre) des coûts de sinistres, et des frais associés.

`>  library(evd); data(lossalae)`

La figure ci-contre montre l’échantillon de 1500 points, avec le logarithme des coûts, et le logarithme des frais. Classiquement, on s’intéresse à la loi conditionnelle des frais sachant le coût du sinistre. On peut alors tracer des boxplot conditionnels (comme pour le coût des sinistres conditionnels à l’âge du véhicule),

Mais l’idée des bagplot est sensemblement différentes. On s’intéresse véritablement au couple de variables.

• Du bagplot fonctionnel

Rob Hyndman a également utilisé ces méthode de bagplot quand la taille devient infiniment grande, i.e. dans un cadre fonctionnel. Par exemple, on peut utiliser cette technique pour détecter des années atypiques sur des évolutions de température, ou dans un ajustement de modèle de type Lee-Carter en mortalité. Ces deux cas sont abordées par Rob Hyndman dans la library(rainbow). Pour faire simple et rapide, on peut définir une profondeur fonctionnelle,

Pour commencer, on peut regarder les température moyenne mensuelles de la température à la surface de la mer, entre 1951 et 2007, obtenue ici.

On note que 4 années se détachent des autres, 1982/1983 et 1997/1998, qui sont particulièrement élevées, mais pas seulement. Certaines semblent décallées de quelques mois. Le code est ici simplement

```> library(rainbow)
> fboxplot(data = ElNino, plot.type = "functional", type = "bag")
> fboxplot(data = ElNino, plot.type = "bivariate", type = "bag")```

Dans le cas des données démographiques, il faut bricoler un peu. Repartons des données français que nous avons étudiées plus d’une fois,

```> DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> DEATH\$Age=as.numeric(as.character(DEATH\$Age))
> DEATH\$Age[is.na(DEATH\$Age)]=110
> EXPOSURE\$Age=as.numeric(as.character(EXPOSURE\$Age))
> EXPOSURE\$Age[is.na(EXPOSURE\$Age)]=110
> YEAR=unique(DEATH\$Year);nC=length(YEAR)
> AGE =unique(DEATH\$Age);nL=length(AGE)
> rownames(MUH)=AGE
> colnames(MUH)=YEAR
> rownames(MUF)=AGE
> colnames(MUF)=YEAR
> MUH=MUH[1:90,]
> MUF=MUF[1:90,]
> MUHF=fts(x = AGE[1:90], y = log(MUH), xname = "Age",
+          yname = "Log Mortality Rate")
> MUFF=fts(x = AGE[1:90], y = log(MUF), xname = "Age",
+          yname = "Log Mortality Rate")
> fboxplot(data = MUHF, plot.type = "functional", type = "bag")
> fboxplot(data = MUHF, plot.type = "bivariate", type = "bag")```

avec un code similaire pour les femmes. Le graphique suivant montre l’évolution des taux de mortalité des hommes en fonctions de leur âge, pour différentes années,

On observe que certaines années se détachent assez nettement, en particulier sur les âges allant de 15 à 40 ans. Comme le montre le graphique à gauche, ce sont les années de guerre qui ressortent fortement (surtout la première guerre mondiale). Pour les femmes, l’évolution des taux de mortalités en fonction des âges semble en revanche beaucoup plus régulière, puisqu’aucune année atypique ne semble ressortir,

Ceci se confirme sur le graphique de droite, où, contrairement aux hommes, aucune année ne semble réellement se détacher. Notons qu’au lieu de chercher des années atypiques, on peut chercher faire une lecture duale, en cherchant les âges qui se ont eut un changement en terme d’évolution du taux de mortalité sur un siècle différent. Par exemple pour les femmes, les taux de mortalité ont l’évolution suivante au cours du siècle,

ce qui permet de repérer quelques âges un peu atypiques, i.e. les bébés, qui ont eu un gain en terme de taux de mortalité considérable, relativement aux autres âges. En revanche, pour les hommes, il n’y a pas particulièrement d’âge qui soit statistiquement différent des autres,

• Applications en terme de stress scenarios

Nous avons vu qu’il y avait presque toujours (pour les hommes comme pour les femmes) des outliers, i.e. des courbes atypiques. Il serait donc normal de s’attendre à en trouver dans les scénarios de prédiction. Les bagplot sont alors un outils intéressant. Si on regarde les prédictions obtenues par la méthode de Lee & Carter par exemple, ou plutôt des scénarios générés par la méthode de Lee & Carter,

```> LCH <- fdm(BASEH,order=10)
> LCHf=forecast(LCH,100,method="ets")
> LCHg=simulate(LCHf,1)
> LCHg1=as.matrix(LCHg[,,1])
> MUHf=LCHg1[1:90,]
> MUHF=fts(x = AGE[1:90], y = log(MUHf), xname = "Age",
+         yname = "Log Mortality Rate")
> library(rainbow)
> fboxplot(data = MUHF, plot.type = "functional", type = "bag")
> fboxplot(data = MUHF, plot.type = "bivariate", type = "bag")```

On peut alors identifier clairement des années atypiques – pour un scenario donné. Bref, on a un outils pour vérifier que l’on génère effectivement des points atypiques, un peu extrêmes, mais surtout, on peut les identifier. On peut alors creuser un peu pour mieux comprendre ce qui se passe….

Pour aller plus loin sur ces histoires de bagplot et aux régressions quantiles, je renvoie ici ou (Thierry Roncalli avait remis la régression quantile au goût du jour il y a quelques années, à l’époque où je programmais encore en Gauss, je renvoie aussi ici pour des applications financières).

• Les quantiles bivariés

En fait, le point de départ du billet était très lié d’une notion qui méritera un autre billet: comment définir des quantiles (et plus généralement des mesures de risques) dans un cadre multivarié….

Incertitude(s) dans les triangles de paiements

Un petit billet aujourd’hui sur les méthodes de provisionnements, et sur les calculs d’incertitude à 1 an, ou à l’ultime, de montants de provisions pour sinistres restant à payer…. Ce sujet avait été évoqué ici (surtout dans les commentaires à l’époque).

• Mack et incertitude à long terme

Le mse dans le cadre du modèle de Mack est proposé dans la library(ChainLadder), via la fonction MackChainLadder(). Sur le triangle utilisé dans les exemples et formations en provisionnement, on obtient les mse estimés suivants

```> source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")

Latest Dev.To.Date Ultimate    IBNR Mack.S.E CV(IBNR)
1  4,456       1.000    4,456     0.0    0.000      NaN
2  4,730       0.995    4,752    22.4    0.639   0.0285
3  5,420       0.993    5,456    35.8    2.503   0.0699
4  6,020       0.989    6,086    66.1    5.046   0.0764
5  6,794       0.978    6,947   153.1   31.332   0.2047
6  5,217       0.708    7,367 2,149.7   68.449   0.0318

Totals
Latest:    32,637.00
Ultimate:  35,063.99
IBNR:       2,426.99
Mack S.E.:     79.30
CV(IBNR):       0.03```

Mais il s’agit fondamentalement d’une formule sur l’estimation ultime de la charge. Pour une année de survenance donnée, la formule dérivée par Thomas Mack ressemble à ça

(honteusement emprunté dans des slides de Mario Wüthrich, par paresse, ici), où pour rappel,

• Muerz & Wüthrich et incertitude à court terme

Pour faire simple, il existe une formule fermée, un peu complexe, et une simplification, reposant sur l’idée que, si on note

alors, Michael Merz et Mario Wüthrich ont montré que

Le code sur le triangle mentionné ici donne le résultat suivant (ce code est parfois appelé dans divers documents MackMerzWuthrich.R)

```> source("https://perso.univ-rennes1.fr/arthur.charpentier/merz-wuthrich-triangle.R")
> MSEP_Mack_MW(PAID,0)
MSEP Mack MSEP observable approchée MSEP observable exacte
1  0.0000000                  0.000000               0.000000
2  0.6393379                  1.424131               1.315292
3  2.5025153                  2.543508               2.543508
4  5.0459004                  4.476698               4.476698
5 31.3319292                 30.915407              30.915407
6 68.4489667                 60.832875              60.832898
7 79.2954414                 72.574735              72.572700```

La dernière ligne (notée ici 7 correspond au total – toutes années de survenance confondues). Ici, on regarde l’erreur à un an. Et Mario Wüthrich concluait – lors d’une jolie école d’été à Lyon il y a quelques années – en notant le point suivant

Sur le triangle étudié ici, on était à 91%. Mais supposer que l’on soit inférieur à 1 est un point qu’il me semble intéressant de discuter…

• incertitude plus grande à court terme ?

Dans un billet datant un peu (ici), on laissait entendre qu’il était sûrement possible que l’incertitude à court terme excède l’incertitude à long terme.
Arnaud Lacoume a trouvé un vrai triangle dans lequel l’incertitude de court terme excède l’incertitude à long terme (détaillé ici dans une feuille xls). Pour information, il s’agit du risque GAV, garanties accidents de la vie, bref, rien de très exotique…

`> Triangle = read.table("https://perso.univ-rennes1.fr/arthur.charpentier/GAV-triangle.csv",sep=";")`

Si on utilise le code gentiment offert par Arnaud, on note qu’effectivement, l’incertitude à 1 an est grande, très grande…

```> MSEP_Mack_MW(Triangle,0)
MSEP Mack MSEP observable approchée MSEP observable exacte
1        0.00                       0.0                    0.0
2    12459.74                  129692.2               152605.9
3   209431.14                  214136.5               214419.6
4   258003.38                  198072.3               198773.0
5  3055297.40                 3048489.5              3065525.1
6 58429393.29                57056117.3             67375794.0
7 58669646.13                57301552.4             67586106.6```

On note ici que le MSE exact à un an (car Merz et Wüthrich proposent aussi une formule approchée, basée sur une approximation qui contraint l’estimation à un an à être inférieure à l’estimation à horizon ultime) est – si on regarde toutes années de survenance confondues – 15% plus élevé que le MSE à ultime. Avec un résultat du même ordre pour la dernière année (ce n’est donc pas un problème d’agrégation). On peut jouer à simuler des triangles, via un processus de déroulé de paiements.

Regardons par exemple le code suivant pour simuler des triangles. L’idée est ici que l’on n’a pas d’incertude sur la valeur à l’ultime,

```> M=matrix(NA,10,10)
> for(i in 1:10){
+ X=rexp(10)
+ Y=cumsum(X)/sum(X)
+ M[i,1:(11-i)]=Y[1:(11-i)]*100}
> Triangle = M```

On tire des variables positives, car on veut un modèle à incréments positifs. Par exemple, une des simulations a donné le triangle suivant

```> Triangle
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 35.4 36.3 38.2 47.6 50.7 66.1 71.3 72.1 98.7   100
[2,] 15.5 20.2 22.3 23.7 63.9 68.2 75.9 79.3 95.1    NA
[3,]  4.1 30.3 32.3 33.3 49.4 53.1 62.8 77.1   NA    NA
[4,]  1.1 19.4 38.9 41.9 48.0 66.5 68.3   NA   NA    NA
[5,]  0.3  7.2 10.8 17.8 28.1 30.1   NA   NA   NA    NA
[6,] 15.6 32.2 46.0 50.6 68.1   NA   NA   NA   NA    NA
[7,]  4.0  9.9 10.6 15.8   NA   NA   NA   NA   NA    NA
[8,]  6.1 24.0 29.4   NA   NA   NA   NA   NA   NA    NA
[9,]  0.0  0.5   NA   NA   NA   NA   NA   NA   NA    NA
[10,]  2.3   NA   NA   NA   NA   NA   NA   NA   NA    NA```

En faisant le calcul sur les mesures d’incertitude à la Merz & Würthrich, on obtient une très grande incertitude à court terme,

```> MSEP_Mack_MW(Triangle,0)
MSEP Mack MSEP observable approchée MSEP observable exacte
1   0.000000                  0.000000               0.000000
2   6.820602                 12.776729             118.050196
3  13.377536                 13.119557              13.126625
4  17.415678                 13.638880              13.651131
5  11.797752                  5.574916               5.580840
6  26.767233                 18.407916              18.429069
7  27.184076                 23.725348              23.773581
8  44.100071                 16.410374              16.444281
9   7.227724                  3.450343               3.460982
10 54.012270                 49.443633              49.593284
11 99.184167                 85.042393             145.008282```

Ici je suis 50% plus élevé ! Bon, mais c’est beau de montrer pourquoi une idée partagée par beaucoup de monde n’est pas toujours vraie…mais il reste à comprendre plus en détails le pourquoi … à suivre donc !

En attendant, je mets ici le code utilisé pour quantifier cette incertitude à un

```### MERZ & WUTHRICH

MWChainLadder <- function(Trian, bool = 0, warning = FALSE){

##  variables
I <- nrow(Trian)
J <- ncol(Trian)
diag <- matrix(NA,I,1)
diag_inv <- matrix(NA,I,1)
S_I <- matrix(NA,1,J)
S_II <- matrix(NA,1,J)
Phi <- matrix(NA,I,1)
Psi <- matrix(NA,I,1)
Delta <- matrix(NA,I,1)
Gamma <- matrix(NA,I,1)
Lambda <- matrix(NA,I,1)
Upsilon <- matrix(NA,I,1)
cov_obs <- matrix(0,I,J)
cov_reel <- matrix(0,I,J)
msep_obs <- matrix(NA,1,J+1)
msep_reel <- matrix(NA,1,J+1)
Delta[1] <- Phi[1] <- Psi[1] <- Upsilon[1] <- Lambda[1] <- 0

Mack\$sigma[J-1] <- sqrt(min(Mack\$sigma[J-2]^4/
Mack\$sigma[J-3]^2,min(Mack\$sigma[J-2]^2,Mack\$sigma[J-3]^2)))
#plot(Mack)

for (i in 1:I){
diag[i] = Trian[i,I-i+1]
diag_inv[i] = Trian[I-i+1,i]
}

for (j in 1:J){
S_I[j] <- sum(Trian[1:(I-j),j])
S_II[j] <- sum(Trian[1:(I-j+1),j])
}
S_I[I] <- 0

Delta[2] <- Mack\$sigma[I-1]^2/(S_I[I-1]*(Mack\$f[I-1])^2)
Phi[2] <- 0
Psi[2] <- Mack\$sigma[I-1]^2/(diag[2]*(Mack\$f[I-1])^2)
Upsilon[2] <- Mack\$sigma[I-1]^2/(S_II[I-1]*(Mack\$f[I-1])^2)
Lambda[2] <- diag[2]*Mack\$sigma[I-1]^2/((Mack\$f[I-1]^2)*S_II[I-1]*S_I[I-1])
for (i in 3:I){
Delta[i] <- Mack\$sigma[I-i+1]^2/(S_I[I-i+1]*(Mack\$f[I-i+1])^2) + sum(
(diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
Mack\$sigma[(I-i+2):(J-1)]^2/(S_I[(I-i+2):(J-1)]*(Mack\$f[(I-i+2):(J-1)])^2))
Phi[i] <- sum((diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
Mack\$sigma[(I-i+2):(J-1)]^2/(diag_inv[(I-i+2):(J-1)]*(Mack\$f[(I-i+2):(J-1)])^2))
Psi[i] <- Mack\$sigma[I-i+1]^2/(diag[i]*(Mack\$f[I-i+1])^2)
Upsilon[i] <- Phi[i] + Mack\$sigma[I-i+1]^2/(S_II[I-i+1]*(Mack\$f[I-i+1])^2)

Lambda[i] <- diag[i]*Mack\$sigma[I-i+1]^2/((Mack\$f[I-i+1]^2)*S_II[I-i+1]*S_I[I-i+1]) + sum(
(diag_inv[(I-i+2):(J-1)]/S_II[(I-i+2):(J-1)])^2*
Mack\$sigma[(I-i+2):(J-1)]^2/(S_I[(I-i+2):(J-1)]*(Mack\$f[(I-i+2):(J-1)])^2))
}
Gamma = Phi + Psi

#MSEP par acc. year
for (i in 1:I){
msep_obs[i] = Mack\$FullTriangle[i,J]^2 * (Gamma[i] + Delta[i])
msep_reel[i] = Mack\$FullTriangle[i,J]^2 * (Phi[i] + Delta[i])
}

#Covariance
for (i in 2:(I-1)){
for (k in (i+1):I){
cov_obs[i,k] <- Mack\$FullTriangle[i,J]*Mack\$FullTriangle[k,J]*(Upsilon[i] + Lambda[i])
cov_reel[i,k] <- Mack\$FullTriangle[i,J]*Mack\$FullTriangle[k,J]*(Phi[i] + Lambda[i])
}
}

#MSEP aggregated
msep_obs[I+1] = sum(msep_obs[2:I]) + 2*sum(cov_obs)
msep_reel[I+1] = sum(msep_reel[2:I]) + 2*sum(cov_reel)

##################################################################################################################
## EXACTES
##################################################################################################################
#DÈclaration des variables
facteur <- matrix(1,I,1)
Phi_exact <- matrix(NA,I,1)
Psi_exact <- matrix(NA,I,1)
Gamma_exact <- matrix(NA,I,1)
Upsilon_exact <- matrix(NA,I,1)
cov_obs_exact <- matrix(0,I,J)
cov_reel_exact <- matrix(0,I,J)
msep_obs_exact <- matrix(NA,1,J+1)
msep_reel_exact <- matrix(NA,1,J+1)
Phi_exact[1] <- Upsilon_exact[1]<- Gamma_exact[1] <- 0

# i=2
Phi_exact[2] <- Psi_exact[2] <- 0
Gamma_exact[2] <- Mack\$sigma[I-1]^2/(diag[I-1]*(Mack\$f[I-1])^2)
Upsilon_exact[2] <- Mack\$sigma[I-1]^2/(S_II[I-1]*(Mack\$f[I-1])^2)

# loop for prior
for (i in 3:I){
facteur[i] <-prod(1 + diag_inv[(I-i+2):(J-1)]*(Mack\$sigma[(I-i+2):(J-1)]^2)/((S_II[(I-i+2):(J-1)]*
Mack\$f[(I-i+2):(J-1)])^2))
Phi_exact[i] <- (1 + Mack\$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack\$f[I-i+1])^2))* (facteur[i]-1)
Psi_exact[i] <- (1 + Mack\$sigma[I-i+1]^2/(S_II[i]*(Mack\$f[I-i+1])^2)) * Phi[i] /
(1 + Mack\$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack\$f[I-i+1])^2))
Upsilon_exact[i] <- (1 + Mack\$sigma[I-i+1]^2/(S_II[I-i+1]*(Mack\$f[I-i+1])^2))* facteur[i]-1
Gamma_exact[i] <- (1 + Mack\$sigma[I-i+1]^2/(diag_inv[I-i+1]*(Mack\$f[I-i+1])^2))* facteur[i]-1
}

#MSEP per accident year
for (i in 1:I){
msep_obs_exact[i] = Mack\$FullTriangle[i,J]^2 * (Gamma_exact[i] + Delta[i])
msep_reel_exact[i] = Mack\$FullTriangle[i,J]^2 * (Phi_exact[i] + Delta[i])
}

#Covariance
for (i in 2:(I-1)){
for (k in (i+1):I){
cov_obs_exact[i,k] <- Mack\$FullTriangle[i,J]*Mack\$FullTriangle[k,J]*(Upsilon_exact[i] + Lambda[i])
cov_reel_exact[i,k] <- Mack\$FullTriangle[i,J]*Mack\$FullTriangle[k,J]*(Psi_exact[i] + Lambda[i])
}
}

#MSEP ag.
msep_obs_exact[I+1] <- sum(msep_obs_exact[2:I]) + 2*sum(cov_obs_exact)
msep_reel_exact[I+1] <- sum(msep_reel_exact[2:I]) + 2*sum(cov_reel_exact)

##################################################################################################################
##Outputs Fonction
##################################################################################################################
msep_Mack <- array(0,c(1,I+1))
msep_Mack[1:I] <- Mack\$Mack.S.E[,I]
msep_Mack[I+1] <- Mack\$Total.Mack.S.E
Vari <- array(0,c(1,I+1))
for (i in 1:I){
Vari[i] <- Mack\$FullTriangle[i,J]^2 * Psi[i]
}
Vari[I+1] <- sum(Vari[1:I])

# complete vision
if(bool==1){
result <- cbind(t(msep_Mack), t(sqrt(Vari)), t(sqrt(msep_obs)), t(sqrt(msep_obs_exact)), t(sqrt(msep_reel)), t(sqrt(msep_reel_exact)))
result <- as.data.frame(result)
l <- list("Variance CDR", "MSEP Mack","MSEP obs. approx.","MSEP obs. exact", "MSEP real approx.","MSEP real exact")
names(result) <- l
return(result)
}

#si l'on souhaite uniquement l'aspect solvabilité
if(bool==0){
result <- cbind(t(msep_Mack), t(sqrt(msep_obs)), t(sqrt(msep_obs_exact)))
result <- as.data.frame(result)
l <- list("MSEP Mack","MSEP obs. approx","MSEP obs. exact")
names(result) <- l
return(result)
}
}```

Delta method avec R, ou transformation des paramètres

Un petit billet aujourd’hui sur une question qui m’a été posée (et qui m’a permis de découvrir des nouvelles fonctions sous R). On suppose que l’on a un modèle économétrique de la forme

que l’on peut encore écrire

Supposons que ce soit la formulation en terme de ces derniers coefficients qui nous intéresse, et que l’on cherche un intervalle de confiance.
La Delta Method nous dit que si on dispose d’un estimateur asymptotiquement gaussien (mais c’est le cas pour la plupart des estimateurs que l’on construit en pratique), tout transformation bijective de cet estimateur reste asymptotiquement Gaussien. Bref, on va pouvoir utiliser ces résultats pour en déduire un intervalle de confiance (asymptotique) pour le nouveau paramètre. En effet, l’estimateur par moindre carrés vérifie la normalité asymptotique (sous quelques hypothèses)

ce qui donne avec les notations précédentes

Autrement dit, on peut facilement en déduire un intervalle de confiance (asymptotique). Pour cela, on a trois méthodes (au moins) avec R. Pour la suite, je vais considerer la régression sur la base de données cars (comme toujours), et on s’intéresse à l’estimation de

(ne me demandez pas pourquoi, il me fallait un exemple numérique, j’ai pris un truc nonlinéaire pour qu’au moins ça puisse avoir un vague intérêt).

• En faisant un peu de calcul formel

Pour les calculs symboliques, il y a la fonction D dans le package de base,

```> f <- expression(sin(cos(x + y^2)))
> ( Dfx <- D(f, "x") )
-(cos(cos(x + y^2)) * sin(x + y^2))```

ce qui pourra alors être utilisé pour faire les calculs de la variance asymptotique.

```> x=c(0,1,2); y=1
> eval(D.sc)
[1] -0.7216061 -0.8316919 -0.0774320```

On peut également utiliser cette fonction pour faire plein de jolis dessins,

```> x <-seq(-5, 5, length=100)
> fx=x^2
> Df <-D(expression(x^2), "x")
> Dfx<-eval(Df)
> plot(x, fx, type = "l", las = 1,ylim=c(-2,2))
> for (i in 1:length(x)) {
+ abline(a = fx[i] - Dfx[i] * x[i], b = Dfx[i],
+ col = rainbow(100)[i])
+ }```

Mais bon, c’était juste pour rigoler un peu… On peut utiliser cette fonction pour calculer l’écart type de l’estimateur de la grandeur qui nous intéresse.

```> reg1 = lm(dist~speed, data = cars)
> estmean = coef(reg1)
> estvar = summary(reg1)\$cov.unscaled * summary(reg1)\$sigma^2
> f = function(x,y){
+ -y/(x+y)
+ }
> X=f(estmean[1],estmean[2])
> g=expression(-y/(x+y))
> x<-as.numeric(estmean[1]); y<-as.numeric(estmean[2])
> V=eval(c(eval(D(g, "x")),eval(D(g, "y"))))
> sqrt(t(V) %*% estvar %*% V)
[,1]
[1,] 0.1063255```

Voilà pour la première méthode. Petite remarque: pour faire le calcul symbolique de la dérivée, j’utilise expression et non pas function.

• En utilisant la fonction deltamethod de la librarie msm

Ici on commence comme avant, sauf que la phase de calcul est immédiat…

```> reg2 <- lm(dist~speed, data = cars)
> estmean = coef(reg2)
> estvar = summary(reg2)\$cov.unscaled * summary(reg2)\$sigma^2
> deltamethod (~ -x2 / (x1 + x2), estmean, estvar)
[1] 0.1063255```

et d’ailleurs, ça peut s’appliquer pour un paramètre en dimension 2

```> deltamethod (list(~ -x2 / (x1 + x2), ~ x1+x2), estmean, estvar)
[1] 0.1063255 6.3664368```

En revanche, on obtient des estimateur des écart-types, et pas de la matrice de variance-covariance des estimateurs. Autrement dit on ne peut pas en déduire l’ellipse de confiance, seulement des intervalles de confiance marginaux (en utilisant la normalité asymptotique).

• En utilisant la fonction delta.method de librarie alr3

Pour finir, une dernière fonction. La seule chose qui change c’est la synthaxe. Mais la sortie est plus élégante,

```> reg3 = lm(dist~speed, data = cars)
> delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255```

En fait, l’écart-type indiqué ici est précisément celui qui est utilisé pour écrire l’intervalle de confiance. Dans cet exemple, l’intervalle de confiance approché (car on suppose que la normalité – en théorie seulement asymptotique – est vérifiée) et estimé (car les paramètres sont inconnus et estimés) est de la forme suivant, pour un seuil de confiance à 95%

```> DM=delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255
> DM\$estimate+qt(.975,df=nrow(cars)-1)*DM\$se
[1] 0.5018277
> DM\$estimate+qt(.025,df=nrow(cars)-1)*DM\$se
[1] 0.07448936```

C’est à dire qu’on a un intervalle de confiance de la forme [0.074;0.502].