Chain Ladder, avec R

Un billet rapide pour mettre en ligne des parties du code tapé en cours, mercredi dernier. On avait commencé par convertir la feuille du classeur excel en un fichier texte, pour faciliter la lecture,

> setwd("C:\\Users\\savsalledecours\\Desktop")
> triangle=read.table("exACT2040.csv",header=TRUE,sep=";")	
> triangle
  ANNEE   X0   X1   X2   X3   X4   X5
1  2000 3209 4372 4411 4428 4435 4456
2  2001 3367 4659 4696 4720 4730   NA
3  2002 3871 5345 5398 5420   NA   NA
4  2003 4239 5917 6020   NA   NA   NA
5  2004 4929 6794   NA   NA   NA   NA
6  2005 5217   NA   NA   NA   NA   NA

L’idée – quand on importe un triangle – est de récupérer une base sous la forme précédente, avec des valeurs manquantes dans la partie inférieure du triangle (on verra l’intérêt quand on fait une régression). On avait ensuite calculé les facteurs de transition, et en même temps complété le triangle,

> T=triangle[,2:7]
> rownames(T)=triangle$ANNEE
> T2=T
> n=ncol(T)
> L=rep(NA,n-1)
> for(j in 1:(n-1)){
+ L[j]=sum(T[1:(n-j),j+1])/sum(T[1:(n-j),j])
+ T2[(n-j+1):n,j+1]=L[j]*T2[(n-j+1):n,j]
+ }

Les facteurs de transition sont ici,

> L
[1] 1.380933 1.011433 1.004343 1.001858 1.004735

et le triangle complété

> T2
       X0       X1       X2       X3       X4       X5
2000 3209 4372.000 4411.000 4428.000 4435.000 4456.000
2001 3367 4659.000 4696.000 4720.000 4730.000 4752.397
2002 3871 5345.000 5398.000 5420.000 5430.072 5455.784
2003 4239 5917.000 6020.000 6046.147 6057.383 6086.065
2004 4929 6794.000 6871.672 6901.518 6914.344 6947.084
2005 5217 7204.327 7286.691 7318.339 7331.939 7366.656

Le montant de provision est ici en faisant la différence entre la charge ultime (dans la dernière colonne) et les derniers paiements observés (sur la seconde diagonale)

> CU=T2[,n]
> Pat=diag(as.matrix(T2[,n:1]))
> Ri=CU-Pat
> R=sum(Ri)

soit, numériquement

> R
[1] 2426.985

On avait alors vu que l’on pouvait calculer un tail factor, en supposant une décroissance exponentielle des facteurs de transition, et on rajoutait alors une colonne correspondant au montant ultime, par année d’accident,

> logL=log(L-1)
> t=1:5
> b=data.frame(logL,t)
> reg=lm(logL~t,data=b)
> logLp=predict(reg,newdata=data.frame(t=6:100))
> Lp=exp(logLp)+1
> Linf=prod(Lp)
> T3=T2
> T3$Xinf=T3$X5*Linf

On a ici

> T3
       X0       X1       X2       X3       X4       X5     Xinf
2000 3209 4372.000 4411.000 4428.000 4435.000 4456.000 4459.149
2001 3367 4659.000 4696.000 4720.000 4730.000 4752.397 4755.755
2002 3871 5345.000 5398.000 5420.000 5430.072 5455.784 5459.639
2003 4239 5917.000 6020.000 6046.147 6057.383 6086.065 6090.366
2004 4929 6794.000 6871.672 6901.518 6914.344 6947.084 6951.993
2005 5217 7204.327 7286.691 7318.339 7331.939 7366.656 7371.862

(je laisse reprendre le code pour calculer le montant de provisions). Enfin, on avait montré comment utiliser une régression pondérée, pour calculer les facteurs de transition,

> T4=as.matrix(T$X0,n,1)
> for(j in 1:(n-1)){
+ Y=T[,j+1]
+ X=T[,j]
+ base=data.frame(X,Y)
+ reg=lm(Y~0+X,weights=1/X)
+ T4=cbind(T4,
+ predict(reg,
+ newdata=data.frame(X=T4[,j]
+ )))
+ }

Ce qui donnait la même projection que la méthode Chain Ladder

> T4
  [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
1 3209 4431.414 4482.076 4501.543 4509.909 4531.263
2 3367 4649.601 4702.758 4723.184 4731.961 4754.367
3 3871 5345.591 5406.705 5430.188 5440.279 5466.039
4 4239 5853.775 5920.698 5946.414 5957.464 5985.673
5 4929 6806.619 6884.435 6914.337 6927.186 6959.986
6 5217 7204.327 7286.691 7318.339 7331.939 7366.656

La suite mercredi prochain, même si on risque d’aller très vite sur la méthode de Mack (et les calculs d’erreur quadratique moyenne pour arriver à la régression de Poisson). A suivre donc…

Easter

This morning, there was an interesting post entitled “why does Easter move around so much?” online on http://economist.com/blogs/economist-explains/…

In my time series classes, I keep saying that sometimes, series can exhibit seasonlity, but the seasonal effect can be quite irregular. It is the cas for river levels, where snowmelt can have a huge impact, and it is irregular. Similarly, chocolate sales (even monthly, or quarterly) depends on Easter. Because it can be either in March, or in April, the seasonal pattern is not as regular as flower sales for instance (Valentine beeing always on February 14th, as far as I remember). If we look at the word eggs on http://google.com/trends/q=eggs…, we do observe a cycle related to Easter.

The title of the article published by http://economist.com/blogs/economist-explains/… claims that there is a lot of variability on Eater’s day. Let us check ! The answer to the question “When is Easter ?” can be the following (if we want a short answer): Easter Sunday is the first Sunday after the first full moon after vernal equinox. For more details, see e.g. http://ortelius.de/east. The algorithm used to compute the date of Easter can is online, on http://smart.net/~mmontes/….

> century = year/100
> G = year % 19
> K = (century - 17)/25
> I = (century - century/4 - (century - K)/3 + 19*G + 15) % 30
> I = I - (I/28)*(1 - (I/28)*(29/(I + 1))*((21 - G)/11))
> J = (year + year/4 + I + 2 - century + century/4) % 7
> L = I - J
> EasterMonth = 3 + (L + 40)/44
> EasterDay = L + 28 - 31*(EasterMonth/4)

Actually, this algorithm can be found in some R packages. Here we use the date of Easter from AD 1000 and AD 3000,

> library(timeDate)
> E=Easter(1000:3000)
> D=as.Date(E)
> table(months(D))/2001

    april     march 
0.7651174 0.2348826

(April being before March, in the alphabetical order) If we look at the distribution of the date, it is the following, the starting point being March 1st,

> J=as.numeric(D-as.Date(paste("01/03/",1000:3000,sep=""),"%d/%m/%Y"))
> hist(J,breaks=seq(20,55),col="light green")

And if we look at the autocorrelation function, we can observe that indeed, after 19 years, there is a strong correlation (that could be seen in the algorithm given previously),

> plot(acf(J))

But in order to get a better understanding of the dynamics, we can also look at transiftion matrices. Define

> Q=quantile(J,seq(0,1,by=.25))
> Q[1]=Q[1]-1
> C=cut(J,Q)

Then, the one year transition matrix is (in %)

> k=1; n=length(C)
> B=data.frame(X1=(C[1:(n-k)]),X2=(C[(k+1):n]))
> (T=table(B$X1,B$X2))

          (20,31] (31,39] (39,46] (46,55]
  (20,31]       0       0     265     277
  (31,39]     316       0      13     182
  (39,46]     224     264       0       0
  (46,55]       1     247     211       0
> P=T/apply(T,1,sum)
> round(P*1000)/10

          (20,31] (31,39] (39,46] (46,55]
  (20,31]     0.0     0.0    48.9    51.1
  (31,39]    61.8     0.0     2.5    35.6
  (39,46]    45.9    54.1     0.0     0.0
  (46,55]     0.2    53.8    46.0     0.0

I.e. if  Easter was early in the year (say in March, in the first quartile), then very likeliy, the year after, it will be late in the year (with 50% chance in the third quartile, and 50% chance in the fourth one).

Le libre, et l’universitaire

On a pu lire ces derniers jours plusieurs articles sur le libre accès aux publications scientifique. Je pense à la tribune “qui a peur de l’open access publiée sur http://lemonde.fr/sciences/… . Etant hébergé par http://hypotheses.org/ qui milite en faveur du libre accès, et d’une plus vaste diffusion des travaux et de la culture scientifique, je n’y suis pas insensible. Sur les publications, la conférence des présidents d’université a aussi publié un document intéressant, intitulé “bien universel, par essence, l’article scientifique, n’est pas un bien comme les autres” en ligne sur http://cpu.fr/

En fait, le but de mon billet aujourd’hui était non pas de discuter sur le libre accès aux publications scientifiques, mais à revenir sur l’utilisation des logiciels libres dans le monde universitaire, et plus spécifiquement sur la place du logiciel libre dans l’enseignement.

En lisant Scott Wilson, qui a publié “open source in higher education: how far have we come?”,  sur http://guardian.co.uk/higher-education-network/… je l’ai trouvé un peu trop optimiste. Dans son article  il note que “universities are ahead of the curve in adopting open source, we should now lead the public sector in exploring its full potential“. Jusqu’à aujourd’hui (ou disons récemment pour être honnête), j’étais assez étranger au débat. Certes, j’utilise des logiciels libres, mais juste parce que je les trouve meilleurs. Ce n’était pas du militantisme, juste du pragmatisme. Et depuis quelque temps, je prends conscience que j’ai (ou que nous avons, en tant qu’universitaires) une responsabilité collective. Par exemple dans mes enseignements, j’utilise R (autant le nommer, ce n’est un secret pour personne que la moitié des billets sur mon blog contiennent du R). Mais quand j’étais étudiant, j’ai appris à manipuler des logiciels de statistiques payant (appelons les S pour utiliser une notation générique qui devrait m’éviter des procès pour diffamation). J’ai aussi utilisé S lorsque j’ai travaillé en tant qu’actuaire, car c’était installé sur nos machines, et parce que j’avais appris ce logiciel. Sans parler de premières amours, disons qu’on garde une forme d’attachement, ou de sécurité, à travailler avec les logiciels que l’on a vus (voire appris) durant nos études. On a peur de changer de logiciel, on aime la stabilité
(car admettons le, on est tous aussi nuls les uns que les autres avec l’informatique). C’est un peu ce que l’on peut lire dans un article paru il y a quelques jours dans http://ledevoir.com/politique/quebec/… où l’on apprend que “Québec adopte un autre décret pour la mise à jour de 76 000 postes de la fonction publique, sans appel d’offres, par des logiciels vendus par […]. La facture totale est de 30 millions, pour les licences seulement. Cette mesure d’exception vise à « assurer la sécurité de nos postes »” (c’est moi qui met en gras). Ce n’est pas de sécurité informatique face à du piratage dont on parle ici, mais de sécurité du personnel qui utilise lesdits logiciels (le confort en quelque sorte). On comprends ainsi que S (et je mets dans le même panier tous les logiciels commerciaux, ou presque) offre des licences gratuites aux étudiants. On retrouve là la logique des fabricants de cigarettes il y a quelques décennies (suis-je optimiste ?) lorsqu’ils essayaient d’attirer les adolescents vers leurs cigarettes, tenant compte ainsi de l’addiction qui se crée ?

Comme je le disais au début de mon paragraphe précédent c’est par pragmatisme que, pour ma recherche, je me suis tourné vers les logiciels libres. Mais pour mon enseignement, je n’avais jamais vraiment réfléchi de manière profonde (parce que je ne fais pas de cours d’informatique: si on utilise un langage c’est parce que l’outil informatique est indispensable pour faire de la modélisation). Par simplicité (ou par fainéantise), j’enseignais des langages que je connaissais. J’utilise R dans mes cours, mais parmi mes proches voisins de bureau, Jean-Philippe utilise S (pour faire des régressions avancées en assurance non-vie) et Mathieu utilise un tableur (que l’on pourrait appeler T, pour…. faire comprendre qu’un tableur suffit pour faire de l’actuariat ? que l’actuariat est plus proche de la comptabilité que de la modélisation statistique ?). Récemment, plusieurs petites phrases m’ont poussé à m’interroger. Je pense que ça a commencé quand un(e) étudiant(e) a écrit dans l’évaluation du cours de modèle de prévision que je donnais voilà quelque temps, que mon cours se résumait à”faire de la saisi de données dans un programme gratuit sur internet” (je passe le fait que dans les questions posées aux étudiants qui remplissent l’évaluation, il est demandé ce que les étudiants pensent des présentations faites avec le logiciel P, comme si faire des présentation en LaTeX n’était pas admissible). Il y a eu ce commentaire l’autre jour sur mon blog, où quelqu’un me demandait le code en S pour refaire ce que je faisais en R. Il y a eu ce questionnement l’autre jour : “est-ce gênant si les étudiants du baccalauréat en actuariat sortent de la formation sans connaître S ?” (que j’avais entendu à Rennes 1 lors des discussions du contenu des cours du Master de statistique et d’économie “il est bon que nos étudiants aient 20 heures de cours à S, ça leur permet de mieux se placer sur le marché du travail“).

Le modèle économique est assez incroyable quand on y pense : avec des fonds publics, l’université offre un service avant-vente à des entreprises qui vendent des logiciels commerciaux. Alors qu’il existe du libre qui fait mieux. Et qui continuera à faire mieux si la communauté d’utilisateurs est active. Je pense que les logiciels commerciaux n’ont pas leur place dans les formations universitaires, et qu’il serait temps que les universitaires prennent leur part de responsabilité dans ce qui s’est mis en place Je ne souhaite pas que les universitaires prennent leur part du profit même si ça serait la moindre des choses (il suffit de survoler la page comparison of statistical packages sur wikipedia pour apprendre que rendre accro un étudiant à un logiciel qu’il utilisera pendant 10 ans, c’est assurer entre 10,000 $ et 60,000 $ de licence (pour un utilisateur), sans parler des formations, des livres, etc). Je pense qu’il y va de notre responsabilité en tant qu’universitaire de promouvoir le logiciel libre. Cela permet déjà de comprendre ce qu’on utilise (je faisais un billet la semaine passée, suite à des soucis rencontrés sur le tableur – commercial – T), mais en plus, comme le notait Scott Wilson, “open source not only promotes creativity but helps make technology more democratic allowing a community to work together to solve common problems.” Et je pense qu’il a raison. Maintenant, je suis nul pour faire des tribunes, et je voulais juste que mes collègues prennent conscience de certains aspects d’un débat qui ne fera que gonfler dans les mois qui viennent…

Somewhere else, part 41

A couple of very interesting posts this week (and also a more funny one),

and a lot of posts and articles worth reading,


“Shaking Tokyo” vostfr (court 29 mn) par ancoliemagic

Sinon, toujours quelques articles intéressants, en français cette fois

Did I miss something ?

Benford law and lognormal distributions

Benford’s law is nowadays extremely popular (see e.g. http://en.wikipedia.org/…). It is usually claimed that, for a given set data set, changing units does not affect the distribution of the first digit. Thus, it should be related to scale invariant distributions. Heuristically, scale (or unit) invariance means that the density of the measure https://latex.codecogs.com/gif.latex?%20X (or probability function) https://latex.codecogs.com/gif.latex?f(x) should be proportional to https://latex.codecogs.com/gif.latex?f(kx). Thus, because densities integrate to 1, the proportionality coefficient has to be https://latex.codecogs.com/gif.latex?k^{-1}, and therefore, https://latex.codecogs.com/gif.latex?f should satisfy the following functional equation, https://latex.codecogs.com/gif.latex?%20kf(kx)=f(x), for all https://latex.codecogs.com/gif.latex?%20x in https://latex.codecogs.com/gif.latex?%20(1,\infty) and https://latex.codecogs.com/gif.latex?%20k in https://latex.codecogs.com/gif.latex?%20(0,\infty). The solution of this functional equation is https://latex.codecogs.com/gif.latex?%20f(x)=x^{-1}, I guess this can be proved easily solving ordinary differential equation

https://latex.codecogs.com/gif.latex?%20\frac{d}{dk}%20(kf(kx))=0

Now if https://latex.codecogs.com/gif.latex?%20D denotes the first digit of https://latex.codecogs.com/gif.latex?%20X, in base 10, then

https://latex.codecogs.com/gif.latex?%20\mathbb{P}(D=d)=\frac{\displaystyle{\int_d^{d+1}%20f(x)dx}}{{\displaystyle{\int_1^{10}%20f(x)dx}}}=\cdots=\frac{\displaystyle{\log\left(1+\frac{1}{d}\right)}}{\log(10)}Which is the so-called Benford’s law. So, this distribution looks like that

> (benford=log(1+1/(1:9))/log(10))
[1] 0.30103000 0.17609126 0.12493874 0.09691001 0.07918125 
[6] 0.06694679 0.05799195 0.05115252 0.04575749
> names(benford)=1:9
> sum(benford)
[1] 1
> barplot(benford,col="white",ylim=c(-.045,.3))
> abline(h=0)

To compute the empirical distribution from a sample, use the following function

> firstdigit=function(x){
+ if(x>=1){x=as.numeric(substr(as.character(x),1,1)); zero=FALSE}
+ if(x<1){zero=TRUE}
+ while(zero==TRUE){
+ x=x*10; zero=FALSE
+ if(trunc(x)==0){zero=TRUE}
+ }
+ return(trunc(x))
+ }

and then

> Xd=sapply(X,firstdigit)
> table(Xd)/1000

In Benford’s Law: An Empirical Investigation and a Novel Explanation, we can read

It is not a mathematical article, so do not expect any formal proof in this paper. At least, we can run monte carlo simulation, and see what’s going on if we generate samples from a lognormal distribution with variance https://latex.codecogs.com/gif.latex?%20\sigma^2. For instance, with a unit variance,

> set.seed(1)
> s=1
> X=rlnorm(n=1000,0,s)
> Xd=sapply(X,firstdigit)
> table(Xd)/1000
Xd
    1     2     3     4     5     6     7     8     9 
0.288 0.172 0.121 0.086 0.075 0.072 0.073 0.053 0.060 
> T=rbind(benford,-table(Xd)/1000)
> barplot(T,col=c("red","white"),ylim=c(-.045,.3))
> abline(h=0)

Clearly, it not far away from Benford’s law. Perhaps a more formal test can be considered, for instance Pearson’s https://latex.codecogs.com/gif.latex?%20\chi^2 (goodness of fit) test.

> chisq.test(T,p=benford)

	Chi-squared test for given probabilities

data:  T 
X-squared = 10.9976, df = 8, p-value = 0.2018

So yes, Benford’s law is admissible ! Now, if we consider the case where https://latex.codecogs.com/gif.latex?%20\sigma is smaller (say 0.9), it is a rather different story,

compared with the case where https://latex.codecogs.com/gif.latex?%20\sigma is larger (say 1.1)

It is possible to generate several samples (always the same size, here 1,000 observations), just change the variance parameter https://latex.codecogs.com/gif.latex?%20\sigma and compute the https://latex.codecogs.com/gif.latex?%20p-value of the test. There might be one tricky part: when generating samples from lognormal distributions with small variance, it might be possible that some digits do not appear at all. On that case, there is a problem with the test. So we just use here

> T=table(Xd)
> T=T[as.character(1:9)]
> T[is.na(T)]=0
> PVAL[i]=chisq.test(T,p=benford)$p.value

Boxplots of the https://latex.codecogs.com/gif.latex?%20p-value of the test are the following,

When https://latex.codecogs.com/gif.latex?%20\sigma is too small, it is clearly not Benford’s distribution: for half (or more) of our samples, the https://latex.codecogs.com/gif.latex?%20p-value is lower than 5%. On the other hand, when https://latex.codecogs.com/gif.latex?%20\sigma is large (enough), Benford’s distribution is the distribution of the first digit of lognormal samples, since 95% of our samples have  https://latex.codecogs.com/gif.latex?%20p-values higher than 5% (and the distribution of the https://latex.codecogs.com/gif.latex?%20p-value is almost uniform on the unit interval). Here is the proportion of samples where the https://latex.codecogs.com/gif.latex?%20p-value was lower than 5% (on 5,000 generations each time)

Note that it is also possible to compute the https://latex.codecogs.com/gif.latex?%20p-value of Komogorov-Smirnov test, testing if the https://latex.codecogs.com/gif.latex?%20p-value has a uniform distribution,

> ks.test(PVAL[,s], "punif")$p.value

Indeed, if https://latex.codecogs.com/gif.latex?%20\sigma is larger than 1.15 (around that value), it looks like Benford’s law is a suitable distribution for the first digit.

Rationality, and MS Excel (and other calculators)

This morning, Mathieu had a nice experience in his course on computational method in actuarial science. But let us start with some mathematical formal definitions.

First, recall that https://latex.codecogs.com/gif.latex?y^x is – somehow – a standard expression. No one should be surprised to see such an expression. Generally (as explained in http://en.wikipedia.org/… ), this function is defined only when https://latex.codecogs.com/gif.latex?y\in\mathbb{R}_+. The idea is that the definition of https://latex.codecogs.com/gif.latex?y^x is that

https://latex.codecogs.com/gif.latex?y^x%20=%20\exp\left(x\log[y]\right)

And it is a definition. Such a function exists only if https://latex.codecogs.com/gif.latex?y\in\mathbb{R}_+ (maybe excluding https://latex.codecogs.com/gif.latex?0). This would be a standard definition in real-analysis.

Now, this ‘power’ function appears also in complex analysis, when dealing with unit roots. From instance, if  https://latex.codecogs.com/gif.latex?z=y^{\frac{1}{k}}e^{i%20\frac{2n\pi}{k}}, where https://latex.codecogs.com/gif.latex?y\in\mathbb{R}_+ and https://latex.codecogs.com/gif.latex?k\in\mathbb{N}_\star, for some https://latex.codecogs.com/gif.latex?n\in\mathbb{N}, then https://latex.codecogs.com/gif.latex?z^k=y. Thus, in complex-analysis it might be more complex to define properly https://latex.codecogs.com/gif.latex?y^x since it might not be unique. But we can relate (sometimes, when https://latex.codecogs.com/gif.latex?x is the inverse of an integer, or maybe a rational number ?) with roots of polynomial functions. So far, nothing new…

Let us get back to Mathieu’s problem. Actually, in his course, he wanted to compute https://latex.codecogs.com/gif.latex?(-8)^{\frac{1}{3}}. With a French version of Excel, entering

you do get https://latex.codecogs.com/gif.latex?-2. If you look at the ‘help’ window, you have some more details

It looks like this hat function can be used to define objects such as https://latex.codecogs.com/gif.latex?y^x. But with

you get

(meaning that this is a problem…). It is also possible to use the power (puissance in French) function of Excel,

Here, you also get

The weird part here is that, in the ‘help’ window, you can read that this power function can be used with any number in https://latex.codecogs.com/gif.latex?\mathbb{R}.

Another point… what about  ? Somehow, it is just the square of the previous one (with the fraction)… Here, typing

you get

(similarly with the power function). So clearly, it is not that simple to use this power function. Now, if you use Google (which is now my new online calculator when I am in class, when I cannot use R), if the power is a fraction (or to be more specific the inverse of an integer), then it works as Excel

 

you get

 But if you type (which should be close, from a continuity property of the power function)

 

you get

and similarly

On Wolfram Mathworld, enter

Mathematica does recognize that we try to deal with unit roots: the result is here

with – as expected – a numerical approximation

With Matlab, Mathieu did obtain the same as Mathematica (its decimal approximation). And to conclude, with R, Mathieu did obtain

> (-8)^(1/3)
[1] NaN
> (-8)^(.333333333333333)
[1] NaN

So for R, you cannot use this hat function on negative numbers.

Now, how can we interpret those outputs ?

1) My understanding is that clearly, with MS Excel, https://latex.codecogs.com/gif.latex?x^{ab}\neq%20\left(x^a\right)^bsince

https://latex.codecogs.com/gif.latex?(-8)^{\frac{2}{3}}\neq%20\left((-8)^{\frac{1}{3}}\right)^2

which is problematic. For instance, in insurance, with monthly discounts, we do have functions like https://latex.codecogs.com/gif.latex?u^{\frac{k}{12}}. What if

https://latex.codecogs.com/gif.latex?u^{\frac{k}{12}}\neq%20\left(u^{\frac{1}{12}}\right)^k

2) The problem comes – probably (MS Excel is not an open software, so it might be hard to check) –  from the fact that https://latex.codecogs.com/gif.latex?y^{\frac{1}{n}} is interpreted as an inverse of a (possibly) bijective function. To be more specific, https://latex.codecogs.com/gif.latex?x=y^{\frac{1}{n}} means that https://latex.codecogs.com/gif.latex?x^n=y. When https://latex.codecogs.com/gif.latex?n is an odd integer, then (in real-analysis) there is a unique inverse, and thus, https://latex.codecogs.com/gif.latex?y^{\frac{1}{n}} is uniquely defined, since https://latex.codecogs.com/gif.latex?x\mapsto%20x^n is a bijective https://latex.codecogs.com/gif.latex?\mathbb{R}\rightarrow\mathbb{R} function. This is what MS Excel (and Google) is doing: https://latex.codecogs.com/gif.latex?x\mapsto%20x^3 is a bijective https://latex.codecogs.com/gif.latex?\mathbb{R}\rightarrow\mathbb{R} function, so https://latex.codecogs.com/gif.latex?(-8)^{\frac{1}{3}} means that we need to find the unique (real) value https://latex.codecogs.com/gif.latex?x such that https://latex.codecogs.com/gif.latex?x^3=-8. Thus, somehow, it makes sense to return https://latex.codecogs.com/gif.latex?-2.

3) There is still a problem with Google, and Mathematica. That is fine to return unit roots in https://latex.codecogs.com/gif.latex?\mathbb{C}. But how comes there is only one value ? I mean, yes https://latex.codecogs.com/gif.latex?1+\sqrt{3}%20\%20i is a possible answer, since

https://latex.codecogs.com/gif.latex?(1+\sqrt{3}%20\%20i)^3=-8

but one can also observe that , and similarly, https://latex.codecogs.com/gif.latex?(-2)^3=-8 and

https://latex.codecogs.com/gif.latex?(1-\sqrt{3}%20\%20i)^3=-8

One can check with

With R, since we do not deal with power function here, but with roots, if we want to find https://latex.codecogs.com/gif.latex?x such that https://latex.codecogs.com/gif.latex?x^3=-8, the function is

> polyroot(c(8,0,0,1))
[1]  1+1.732051i -2+0.000000i  1-1.732051i

Which is different… Weird isn’t it ?

Somewhere else, part 40

Once again, a lot of interesting posts, here and there (but, as usual, outside this blog)

 

Total fatalities, in the U.S. per day (above) or per month (below)

et toujours quelques articles en français, dont un billet passionnant

  • “Il ne faut pas prendre des gens pour des cons mais ne jamais oublier qu’ils en sont” http://sametmax.com… très instructif (car oui, je suis un con)

et une information relevé par tomroud (qui nous confirme une fois de plus qu’il suit vraiment tout)

enfin, pour les amateurs

Did I miss something ?

Somewhere else, part 39

A very interesting post discovered a few days ago,

and as usual, several posts and articles, here and there,

avec comme souvent quelques articles en français

Did I miss something ?

Solvabilité et provisionnement

Mercredi, nous allons aborder en cours les aspects de solvabilité des compagnies d’assurance IARD. Plus particulièrement, nous parlerons des provisions pour sinistres à payer, ou “provision for claims outstanding (PCO)” selon la terminologie anglaise, i.e. “the estimated total cost of ultimate settlement of all claims incurred before the date of record, whether reported or not, less any amounts already paid out in respect thereof.” Je renvoie à la lecture de Le contrôle de la solvabilité des compagnies d’assurance  en ligne sur le site de l’OCDE, pour une vision globale des approches de ces provisions. La SOA avait publié un rapport en 2009, Comparison of Incurred But Not Reported IBNR Methods que j’encourage à lire.

Nous aborderons mercredi les triangles. Parmi les triangles que nous manipulerons

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

qui contient plusieurs fichiers, dont

> PAID
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372 4411 4428 4435 4456
[2,] 3367 4659 4696 4720 4730   NA
[3,] 3871 5345 5338 5420   NA   NA
[4,] 4239 5917 6020   NA   NA   NA
[5,] 4929 6794   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

ainsi que le triangle évoqué sur http://rworkingparty.wikidot.com/

> OthLiabData = read.csv("http://www.casact.org/research/reserve_data/othliab_pos.csv",header=TRUE, sep=",")
> library(ChainLadder)
> OL = SumData=ddply(OthLiabData,.(AccidentYear,DevelopmentYear,DevelopmentLag),summarise,IncurLoss=sum(IncurLoss_h1-BulkLoss_h1),
+ CumPaidLoss=sum(CumPaidLoss_h1), EarnedPremDIR=sum(EarnedPremDIR_h1))
> LossTri = as.triangle(OL, origin="AccidentYear",
+ dev = "DevelopmentLag", value="IncurLoss")
> Year = as.triangle(OL, origin="AccidentYear",
+ dev = "DevelopmentLag", value="DevelopmentYear")
> TRIANGLE=LossTri
> TRIANGLE[Year>1997]=NA
> TRIANGLE
      dev
origin      1      2      3      4      5      6      7      8      9     10
  1988 128747 195938 241180 283447 297402 308815 314126 317027 319135 319559
  1989 135147 208767 270979 304488 330066 339871 344742 347800 353245     NA
  1990 152400 238665 297495 348826 359413 364865 372436 372163     NA     NA
  1991 151812 266245 357430 400405 423172 442329 460713     NA     NA     NA
  1992 163737 269170 347469 381251 424810 451221     NA     NA     NA     NA
  1993 187756 358573 431410 476674 504667     NA     NA     NA     NA     NA
  1994 210590 351270 486947 581599     NA     NA     NA     NA     NA     NA
  1995 213141 351363 444272     NA     NA     NA     NA     NA     NA     NA
  1996 237162 378987     NA     NA     NA     NA     NA     NA     NA     NA
  1997 220509     NA     NA     NA     NA     NA     NA     NA     NA     NA

Examen intra, éléments de correction

L’énoncé de l’examen intra est en pdf ici et comme annoncé par courriel, la correction de l’intra est dans le pdf en ligne. Comme personne ne semble en désaccord avec les réponses proposées, les notes seront mises en ligne très bientôt. Concertant les questions 18 et 19 quelques compléments d’explications (que je n’avais pas tapé dans le pdf). On avait vu que l’estimateur du maximum de vraisemblance pour une régression de Poisson était asymptotiquement Gaussien,

https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\beta}}_{P}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{P}))

(asymptotiquement) avec

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}

Quand on a une régression de type binomiale négative, si on note de manière très générale https://latex.codecogs.com/gif.latex?\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i) (on avait vu en cours qu’il existait plusieurs spécifications possibles pour cette variance conditionnelle). Dans ce cas,

https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\beta}}_{BN}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{BN}))

avec

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}\left[\sum_{i=1}^n%20\omega_i%20\boldsymbol{X}_i\boldsymbol{X}_i\right]\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}

Bref, tout dépend fondamentalement de la spécification de la variance conditionnelle. Sous R, c’est la régression binomiale négative de type 1 qui est considérée, i.e.

https://latex.codecogs.com/gif.latex?\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i)=\phi\cdot%20\mathbb{E}(Y_i|\boldsymbol{X}_i)=\phi%20\cdot%20\widehat{Y}_i

On toujours une relation de la forme

https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\beta}}_{QP}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{QP}))

avec (en simplifiant un peu)

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}

aussi, on a

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot%20V_\infty(\widehat{\boldsymbol{\beta}}_{P})

Mais comme annoncé en cours, des points étaient données pour ceux qui se contentaient d’affirmer que la variance des estimateurs était plus grande s’il y avait sur-dispersion.

Happy St Patrick’s Day

I love Saint Patrick’s Day for, at least, two reasons. The first one is that, on March 17th, you can play out loud The Pogues, the second one is that it’s the only day in the year when I really enjoy getting a Guiness in a pub. And Guiness is important in statistical science (I did mention a couple of hours ago – on this blog –  that beers were important for social reasons in the academic world, but that was for other reasons…)

> theta=seq(0,pi/2,length=101)
> leaf=sin(2*theta)+.25*sin(6*theta)
> for(k in 0:3)
+ polygon(leaf*cos(theta+k*pi/2),leaf*sin(theta+k*pi/2),col="green")

As mentioned in all my statistics and econometrics courses, the history of statistics (I mean here mathematical statistics) is closely related to Guinness.

A long time ago, there was a Guinness Brewing Company of Dublin, which – as its name suggests – was an Irish brewing company. And the boss, who was to inherit the family business, decided to attract young students, trained in chemistry at Cambridge or Oxford.

In 1899, William Sealy Gosset, who had obtained a double degree in math and chemistry, left Oxford to Dublin. And to be quite honest, being graduate in maths meant when he had studied differential equations and astronomy. Basically, mathematics were useless for Guinness, and he got there with his expertise in chemistry. In fact, William turned out to be also a very good administrator, but this has nothing to do with our story.

William had good memories of his studies in math, and he wondered if he could find a problem to look at. He started studies on workmanship, noting that conditions vary so much (temperature, from hops, malt, manufacturing conditions …) that there were only few consistent data. The “law of errors”  (the central limit theorem) can not apply under these conditions.

In short, Bill (now we know each other a little, we’ll call him Bill) took many measurements, and noticed that the Poisson distribution could be an interesting model to work with. To make the story short, Bill managed to use statistical techniques to control the variance of the production, meaning that he was able to lower losses in the production of beer.

A nice application like this one deserved publication in a scientific journal … Well, of course the Poisson distribution has long been known (it was 1904 and a few months before, Von Bortkiewicz found elegant applications of this law, as discussed in a post  a few weeks ago). But there was a disclosure issue there: Bill’s contract prohibited him from disclosing secrets to the competitors.

Meanwhile, Bill had met Karl Pearson, who was then editor of Biometrika, and encouraged him to publish his results. In 1906, Bill who had helped Guiness to gain a lot of money – doing applied mathematics can be usefull – managed to take a sabbatical to work with Pearson to Galton Laboratory biometrics. Bill and Karl decided to publish the work under a pseudonym “Student.” The legend claims that they had hesitated to use “pupil.”

And for almost 30 years, “Mr Gosset” honorable employee Guinness led a dissolute life by publishing in statistical journals (after work in the brewery) always under the pseudonym “Student”. Of course, it might not be that simple. I mean, Bill had a family life, too. And his wife was the captain of the national Hockey team. So I hardly imagine Bill playing the smart ass and doing mathematical computations, when it was time to wash the dishes or iron his shirt…

In 1908, he wrote a remarkable “the probable error of the mean” remarked, at least, by Ronald Fisher. In fact, Bill found that there was a interesting law, but – as the normal – it was difficult to manipulate to obtain confidence intervals. Without a computer, he had the idea of ​​using monte carlo methods to tabulate quantiles and construct its tables. And he was probably the first one to look carefully at the problem of small samples, unlike Karl Pearson, who always put focus on the asymptotic case.

In fact, looking at his small sample, he saw the denominator magnitudes very close to those specifically manipulated Karl, in particular a square root of chi-square law. Well, of course, remained the normality assumption, but at least we had some results for finite samples !

For the story, William Gosset suggested to use letter z for its statistics, the ratio between the mean and (empirical) standard deviation. But a few years later, statisticians became accustomed to use this letter for Gaussian distribution (i.e. when the variance is known), and it became the standard to use the letter t. Hence finally the present name of “Student-t distribution” and in regression outputs, we have the “t-test”.

A legend (told by Harold Hotelling in his memoirs) claims that the Guinness family discovered this double life on the day of the death of William Gosset in 1937 when mathematicians requested financial assistance to print a volume of the works of their employee. But another legend claims that Mr Guinness himself would have suggested his nickname when he had expressed his intention to publish his research… So I guess we’ll never know. But at least, I’ll think about Bill when I’ll get my first Guiness tonight (but I will probably not be able to tell this story anymore when I’ll reach the fourth…)

Bloguer, et autres activités d’un enseignant-chercheur

J’avais promis un billet bloguer c’est nul dans mon dernier billet sur la recherche et l’activité de blogger. Il va devoir attendre un peu… Dans ce billet, j’avais fait un diagramme pour expliquer la vie et le partage du temps pour un enseignant-chercheur et blogger. Bien entendu, c’est plus compliqué ! Déjà pour un enseignant-chercheur normal, la vie se partage en

  • 50% recherche
  • 50% enseignement
  • 30% tâches administratives diverses

(et oui, c’est dur à représenter dans un diagramme, Tom Roud ou Mix la Malice avaient tenté un décompte du temps passé dans chacune des activités, je renvoie à leurs billets respectifs). Dans ces tâches divers, on retrouve dedans se trouvent les taches officielles, qui figurent sur le CV, comme faire partie du comité bibliothèque ou du comité de programme,voire être responsable d’une programme de maîtrise, mais surtout les tâches non-officielles (en tout cas que parfois on évitera de mettre sur son CV), comme mettre du scotch-tape sur ses factures et discuter avec les services financiers pour comprendre pourquoi le remboursement d’un repas lors d’une conférence est refusé alors que le sandwich à l’aéroport est remboursé ou remplir un sondage en ligne pendant 40 minutes, commandité par un organisme subventionnaire qui veut savoir si le financement de la recherche fonctionne bien (que vous faites car vous voulez expliquer qu’il faut continuer à financer la recherche, mais vous ne trouvez pas la section où vous voulez insister sur le fait que perdre 40 minutes à remplir un sondage ou 2 heures pour vous faire rembourser 6.37 dollars pour un sandwich, c’est un peu du temps pris sur le temps qu’on veut passer à faire de la recherche). Je ne vais pas passer des heures à constituer un bestiaires de ces activités diverses, je renvoie à http://laviedemix.over-blog.com/, le merveilleux blog de mixlamalice, qui fait cela très bien !

Une autre activité (que je compte dans ces tâches diverses) que j’ai pratiqué plusieurs années lorsque j’étais (co)responsable de Master 1 en France, c’est faire jouer les réseaux. Pour faire des recommandations pour des étudiant(e)s qui souhaitent partir dans des programmes plus sélectifs (j’utilise le terme recommandation car il n’y a pas que la lettre que l’on envoie, il y a aussi souvent le coup de fil – ou le courriel – que l’on passe si  on pense vraiment que l’étudiant(e) le mérite), ou pour aider un étudiant à avoir un stage. Et dans ce cas, avoir un réseau, et l’entretenir, ça sert. Pour les stages en entreprises, il y a les copains de promos (j’ai fait une école il y a quelques années), les anciens collègues (j’ai travaillé dans le privé), mais il y a aussi les formations professionnelles que l’on peut faire, les invitations à déjeuner (pour parler boulot, ou à prendre une bière ou deux si on veut en plus le faire sur un mode plus détendu). Dans le monde de la recherche, ce réseau se crée lors de conférences, lors de visites de collègues. Et c’est un peu là que l’on commence à dépasser les 100% de la vie professionnelle dans mon décompte introduisant ce billet: les bons réseaux sont constitués de collègues mais surtout d’ami(e)s. De collègues avec qui ont pourra passer des vacances, qui restent squatter à la maison quand on les invite parler au séminaire (voire faire une semaine de recherche). Ceux qui connaissent l’adresse du blog de mon fils verront dans les photos de vacances des responsables de master, des co-éditeurs de revues, des organisateurs de conférences (voire des moitiés de théorèmes), en maillot de bain ou en train de faire la cuisine… Mais je m’égare. Mon point était que les activité en ligne sont un bon outil pour avoir, ou entretenir, un réseau. Et n’ayant pas de page facebook, le blog est un outil merveilleux pour ça !

Petit exemple pas plus vieux que la semaine dernière. Ma femme (qui est une mathématicienne) recoit un courriel d’un ancien collègue: “j’ai un bon étudiant qui cherche un stage, de recherche, si possible au Canada“. Alors je vais dire un truc avec des pincettes pour éviter de froisser des collègues, mais au Canada, quand on demande un stage à un prof, c’est un peu comme faire la manche de le métro aux heures de pointe à Paris: “salut, j’ai été contactée par un étudiant qui cherche un stage et…“, “oh là là, ça va être dur, j’ai pas l’sous“. En France, ça serait “j’ai pas le temps” mais au  Canada, c’est les sous que l’on mentionne avant tout (cela dit, en France non plus les sous ne sont pas là , mais les règles sont claires et un(e) étudiant(e) qui fait un stage à l’université se doute bien qu’il n’aura rien). Bref, les prises de contact de ma femme n’ont pas été très concluantes… Comme j’évite de parler boulot avec ma femme, il a fallu un peu de temps pour qu’elle évoque cette histoire, en m’expliquant un peu les thèmes que voudraient aborder l’étudiant. Je n’insisterais pas sur le fait que ma première question a été “il a une page internet ?” (ma femme a du me rappeler que les gens normaux n’ont pas de page internet quand ils sont en maîtrise). Mais j’ai glissé “tu sais, il  y a un super prof à Vancouver qui bosse sur ces thèmes là… on se retweete l’un l’autre sur twitter, et il a un super blog… je peux lui envoyer un mot sur twitter, voir si ça l’intéresse…“. Bref, j’envoie un message (privé). Et vendredi, ma femme apprenait que l’étudiant avais pris contact, que le contact semblait bien passé, et donc il devrait être cet été face au Pacifique…

La morale de mon histoire ? L’étudiant avait un CV qui parlait pour lui, et je n’ai pas fait grand chose. Et sans moi, l’étudiant aurait probablement trouvé, malgré tout, un stage passionnant. Mais le blog permet d’entrer dans un réseau de blogueurs (ou de blogueurs-lecteurs de blogs), avec qui on peut discuter, qui permet d’entrer en contact avec des personnes peut-être plus ouvertes que d’autres enseignants-chercheurs. Peut-être est-ce une illusion, peut-être qu’à force de les lire, d’apprendre sur leur vie privée, on a l’illusion de les connaitre… Mais mes expériences ont toujours été bonnes… Les blogs permettent d’entrer en contact avec une communauté de gens passionnants. Les contacts restent souvent virtuels comme lorsque je discute avec David MonniauxBaptiste CoulmontGizmoAlexandre Delaigue, Joel Gombin ou Olivier Bouba-Olga (damned, je me rends compte que je vais en oublier un paquet) mais j’ai toujours du plaisir à aller manger avec Tom Roud, je dois toujours une bière à Reka, et un resto à Camille (entre autres afin de découvrir qui se cache derrière cette identité secrète, c’est mesquin, je sais, mais c’est pas moi qui ait commencé). Derrière les blogueurs se cachent souvent des gens curieux et passionnés (dans le meilleur des cas passionnants), et c’est pour  se rapprocher de cette communauté qu’il est intéressant de bloguer. Bon, j’en conviens, il y a aussi des blogueurs imbus d’eux-même et pénibles… Mais je ne vais pas commencer mon billet bloguer c’est nul aujourd’hui.

Somewhere else, part 38

Some interesting posts, here and there, discovered this week,

Du coté des blogs francophones (et des journaux en ligne)

et un moment de détente avec

Did I miss something this week ?

Des tablettes… aux blogs ?

 Je viens de finir ‘du scribe au savant‘ d’Yves Gingras, Peter Keating et Camille Limoges. Il m’aura fallu du temps car j’ai l’impression d’avoir de moins en moins de temps pour lire ces temps ci (ou alors des articles ou des rapports qu’il faut évaluer). Pourtant le livre est tout petit quand on y pense. Je veux dire par là qu’on part de la Mésopotamie à l’histoire récente (disons, jusqu’à 1850 environ), en parcourant l’histoire des maths, de la médecine, des sciences naturelles, etc. Forcément, c’est sobre et court, alors quand une section nous intéresse moins, on prend sur soi une page (ou deux maximum), et quand au contraire on se passionne, on reste un peu sur sa faim. J’y ai trouvé une porte ouverte vers plein d’ouvrages même s’il a fallu que je trouve des références par moi-mêmes: à la fin de chaque paragraphe, on attend une référence pour aller plus loin ! C’est un peu ce qui est fait dans le chapitre sur les romains (ce n’est pas le chapitre que j’ai préféré), mais les références sont plus sommaires dans les autres chapitres. J’ai adoré le premier chapitre, qui met en parallèle l’Égypte et de la Mésopotamie, et j’ai vraiment appris plein de choses. J’ai beaucoup aimé l’idée de mettre le scribe (spécialiste de l’écriture et attaché aux inventaires administratifs du royaume) à l’origine de la naissance de la science. Ce qui met un peu l’actuaire à l’origine de la connaissance et du savoir, « actuarius: scribe chargé de la rédaction des procès-verbaux », et j’aime quand on essaye de me raconter une histoire qui me caresse dans le sens du poil. J’ai adoré lire une page (voire deux) sur l’importance de la connaissance des régimes des vents (en fonction des saisons, et des hémisphères) pour les grandes explorations (point dont je n’avais jamais pris conscience, avant de le lire). Ou sur la naissance des universités en Europe (beaucoup moins religieuses que je ne le croyais, mais là encore, à peine deux références…). En plus les schémas sont sobres, et clairs. Ce livre parle du savoir, et de sa transmission (essentiellement en occident, là aussi je suis resté un peu sur ma fin, mai si le livre avait été plus gros, je pense que j’aurais eu peur). J’attends la suite – même si ça ne serait plus réellement un liste d’histoire – qui devrait arriver sur les modes contemporains de transmission du savoir. Avoir un point de vue d’historiens des sciences enrichirait la discussion. Bref, pendant 3 semaines, j’ai appris des choses, tous les jours, et j’avoue qu’on y prend gout !

Examen intra, régression logistique et de Poisson

L’examen intra du cours ACT2040 aura lieu mercredi matin, de 9:00 à 12:00. Aucun document autorisé, sauf les calculatrices (modèle standard, cf plan de cours), et les téléphones seront formellement interdits. Il y aura 34 questions portant sur la première partie du cours (jusqu’à la fin des modèles de comptage, sections 1 à 5 des transparents). 15 questions porteront sur la base décrite dans un précédant billet, sur le nombre de relations extra-conjugales. Il s’agira de décrire les sorties en ligne ici. Je laisse 36 heures pour prendre connaissance de ces sorties. Une version sera donnée lors de l’examen (imprimée 2 pages par feuille, comme dans la version en ligne: si quelqu’un a besoin d’un exemplaire imprimé plus gros, merci de me le faire savoir avant l’examen).