Arthur Charpentier

Aller au contenu | Aller au menu | Aller à la recherche

vendredi 12 mars 2010

On the return period of the 2003 heat wave

The paper on the 2003 heat wave has been submitted a few months ago, and now, the preprint can be downloaded (here for the preprint version, or here on Hal). The paper has a long story, that I will probably summarize here, some day...
"Extremal events are difficult to model since it is difficult to characterize formally those events. The 2003 heat wave in Europe was not characterized by very high temperatures, but mainly the fact that night temperature were no cool enough for a long period of time. Hence, simulation of several models (either with heavy tailed noise or long range dependence) yielddifferent estimations for the return period of that extremal event.  "
The paper has been presented in several conferences, in started a couple of days after the birth of Romane, in March 2005, in talks on dependencies in environmental series....

Reinsurance, ruin and solvency issues: some pitfalls

The paper on solvency, ruine and reinsurance has now been officially submitted (here for the preprint version, or here on Hal). The goal was to prove that the following statements, "if we choose ruin probability as a risk measure, the goal of reinsurance is to reduce this probability to a certain chosen level" (that can be found in the Actuarial models: the mathematics of insurance ), or "reinsurance is able to offer additional underwriting capacity for cedants, but also to reduce the probability of a direct insurer’s ruin" (from the Encyclopaedia of Financial Engineering and Risk Management) might be false. Reinsurance can actually increase ruin probability.
"In this paper, we consider optimal reinsurance from an insurer’s point of view. Given a (low) ruin probability target, insurers want to find he optimal risk transfer mechanism, i.e. either a proportional or a nonproportional reinsurance treaty. Since it is usually admitted that reinsurance should lower ruin probabilities, it should be easy to derive an efficient Monte Carlo algorithm to link ruin probability and reinsurance parameter. Unfortunately, if it is possible for proportional reinsurance, this is no longer the case in nonproportional reinsurance. Some examples where reinsurance might increase ruinprobabilities are given at the end, when claim arrival and claim size are not independent. "
Note that this paper is a revised version of the work presented at the RESIM (rare event simulation) conference 18 months ago (here).

Forme des prénoms des français(e)s

Allez, un dernier billet sur les prénoms des français(e)s....

  • Etude sur la taille des prénoms
On peut - pour s’amuser - regarder la taille des prénoms, l’intuition étant qu’il y a de plus en plus de prénoms courts. La longueur moyenne des prénoms est représentée ci-dessous (l’axe des ordonnées correspondant au nombre de lettres moyen dans le prénom)

avec également ci-dessous la proportion des prénoms de moins de 4 lettres (en rouge), et la proportion des prénoms de plus de 9 lettres (en bleu),

On retrouve effectivement la baisse des prénoms longs depuis les années 60, et surtout une explosion des prénoms courts depuis 30 ans...

  • Etude des lettres de début et de fin dans les prénoms

Pour conclure, on peut aussi se demander si la forme des prénoms a changé, en particulier sur les premières et dernières lettres. Ci dessous, la proportion des prénoms qui commencent par un E (en rouge), ou qui finissent par un E (en bleu),

on note que de moins en moins de prénoms finissent par un E aujourd’hui. En revanche, de plus en plus finissent par un A,

On peut assui observer qu’il y a une mode des prénoms commençant par un S dans les années 70,

mais bon, ce genre d’étude peut durer des heures....
  • Etude de la succession des lettres dans les prénoms
Au début du 20ème siècle Andrej Markov a tenté de formaliser des modèles de « probabilités en chaînes ». A deux occasions, il a proposé des applications simples, sur l’étude de la succession de voyelles et de consonnes dans des textes russes. En 1913, il publiait un exemple de recherche statistique sur le texte Eugénie Onéguine illustrant la liaison des épreuves en  chaîne, où il étudiait les 20 000 premières lettres de l’ouvrage de Pouchkine.10 ans plus tard, il travaillait sur les 100 000 premières lettres d’un roman d’Aksakov. Micheline Petruszewycza publié un livre sur le sujet en 1981.
Dans Eugénie Onéguine, Markov avait obtenu la matrice de transition suivante

voyelle consonne
voyelle 1104 7534
consonne 7534 3837
autrement dit, 7534 fois, une consonne suivait une voyelle. Ce genre d’étude a été reprise à plusieurs reprises, en particulier pour étudier les différences entre les langues (les langues qui disposent d’un alphabet). On peut aussi utiliser cette approche pour voir des changements dans le temps de la forme des mots. 
Sur le graphique ci-dessous, on regarde la probabilité qu’une voyelle soit suivie d’une consonne (en rouge) et  la probabilité qu’une voyelle soit suivie d’une voyelle (en bleu),

On notera que les prénoms ont eu tendance par le passer à jouer davantage sur la répétition de voyelles, mais que cette mode semble toutefois passé depuis un trentaine d’année,

On notera également que les prénoms présentant des successions de consonnes se font de plus en plus rares,

Bref, les prénoms français ont décidément beaucoup changés.... Et maintenant j’arrête de jouer avec mes données.

mercredi 10 mars 2010

Solvency II'newspeak and claims reserving

I should give a talk on thursday afternoon at the internal seminair of the AXA Group. The subject of the talk will be "Solvency II’ newspeak: one year uncertainty for IBNR". In an old paper, Claude Bébéar mentioned that actuaries loved newspeak (or novlangue in French) here or there. Classically when focusing on IBNR uncertainty, we calculate the conditional prediction mean squared error as a quantitative measure of uncertainty, as done in Mack (1993),

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.Diapositivecdr1_m.jpg
The idea in Solvency II is to look one year foreward, assuming that next year’s value was known,
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.DiapositiveCDR2_m.jpg
and thus, we consider the difference between the two predictors: the one obtained now, and the one obtained next year (or to me more specific the ones that we might obtain new year since the blue point is unknown, as at now)
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.Diapositivecdr3_m.jpg
The aim of the talk is to describe what the quantity means, in much more simplier cases than IBNR calculation, and to see how to use bootstap technique to estimate that quantity. The slides can be found here.

Inventaire et Zillmérisation

Au lieu de faire une réponse au commentaire sur les PM (ici), je vais faire un billet. La question était "cettte  provision correspond elle à la PM zillmérisée. Si non, comment la calcule-t-on ? ". Finalement, ça correspond aussi à la suite du billet d’hier sur les termes actuariels (ici). Avant de parler de la PM, il faut revenir un peu davantage sur les primes, et distinguer la prime pure (que l’on a évoquée jusqu’à présent) et la prime commerciale. Ensuite seulement, on pourra parler des calculs de PM.

  • Primes pure et commerciale
Jusqu’à présent, je n’ai pas parlé des frais de gestion ou des frais d’acquisition. Dans le second cas, il s’agit des dépenses engagées pour établir la police, et rémunérer l’agent qui rapporte le contrat. Dans le premier cas, il faut distinguer les frais durant la période où il faut encaisser les primes (période 1), et les frais durant la période où il faut garantir un paiement à l’assuré (période 2). Au cours de la première période, on trouvera les frais d’établissement des quittances, par exemple. Les deux peuvent bien entendu coïncider, auquel cas les frais s’additionne. On supposera que les frais sont proportionnels aux capitaux garantis. On peut aussi supposer qu’il y a des frais de prime commerciale (proportionnels à la prime pure)
Sous ces hypothèses, pour calculer la prime chargée (ou prime commerciale) de l’assuré, on suppose que la valeur actuelle probable des engagements de l’assuré est égale à la valeur actuelle probable des engagements de l’assureur, à laquelle on ajoute la valeur actuelle probable des frais de gestion.
Si on considère un prime unique, on a quelque chose comme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-01.png
où http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-02.png désigne la prime commerciale unique, http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-03.png désignent les frais d’acquisition, http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-04.png correspond aux frais de gestion des engagements de l’assureur (il n’y a pas de http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-05.png car la prime est versée au début, en une seule fois),http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-06.png l’annuité associée à la garantie proposée par l’assureur (i.e. avec un capital garanti de 1) et r les chargements sur la prime. On suppose ici que tout est payé en début d’année, histoire de simplifier...
Pour une prime annuelle, en notant en minuscule les primes annuelles, et http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-07.png l’annuité associée aux engagements de l’assuré (i.e. payer sa prime annuellement tant qu’il est en vie pendant une période prédéterminée),
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-08.png
soit encore
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-50.png
Pour compliquer un peu les choses, notons qu’il existe une contre-assurance de cette prime commerciale.En effet, certains contrats d’assurance en cas de vie (disons de la retraite pour fixer les choses) où une garantie prévoie un remboursement des primes commerciales versées en cas de décès de l’assuré avant le début de la retraite. Toujours pour fixer les choses, considérons une assurance de capital différé, de durée http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png, souscrite à l’âge http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-11.png. Comme auparavant, on note http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-02.png la prime unique commerciale contenant cette contre assurance, de telle sorte que
La valeur actuelle probable de la contre-assurance correspond à l’annuité d’une assurance temporaire décès de http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png années, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-12.png
qui peut s’écrire
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-13.png
soit
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-14.png
Si l’on considère une prime annuelle, payée tout au long des n années, on a quelque chose de la forme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-15.png
où pour rappel (je renvoie au cours à ce sujet) http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-16.png désigne une assurance temporaire décès dont le capital croît d’une unité chaque année, ce qui donne
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-17.png
Ce n’est pas forcément joli, mais comme toujours en actuariat de l’assurance-vie, il suffit d’écrire les choses calmement pour retomber sur des notions qui ont été prédéfinies.
  • PM pure, d’inventaire ou zillmérisée
De la même façon qu’il existe plusieurs formes de chargements et de frais pour les prime, il existera plusieurs sortes de provisions mathématiques. Telle que je l’ai présentée, i.e. différence entre la valeur actuelle probable des engagements de l’assureur, et la valeur actuelle probable des engagements de l’assuré (des primes), il s’agit de la provision mathématique pure. Pour un contrat souscrit à l’age x, et vu après k années, je l’avais notée http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-18.png (c’est la notation officielle). Mais ce n’est pas ce que la réglementation demande de calculer.
"Les provisions mathématiques constituées par les entreprises d’assurance-vie et de capitalisation sont calculées en tenant compte, dans la détermination de l’engagement de l’assuré ou du souscripteur, de la partie des primes devant être versée par l’intéressé et représentative des frais d’acquisition du contrat, lorsque ces frais ont été portés en charge déductible par l’entreprise avant la fin de l’exercice à la clôture duquel la provision est constituée." (Article L331-1 du Code des Assurance, ici)
L’idée est qu’une fois la période de paiement de la prime (voire dès le début en cas de prime unique), la PM pure ne sera pas suffisante, l’assureur ayant davantage de frais.
Il va falloir prendre en compte différentes sortes de frais ou de chargements (évoqués auparavant). Le cas le plus simple est celui des chargements dits de règlement: si pour payer 1€ de prestation, cela coûte http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-20.png€ à l’assureur, la prime annuelle payée par l’assureur doit être multipliée parhttp://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-20.png (de telle sorte que les VAP initiale soient égales). Dans ce cas, un rapide calcul montre que la PM doit être multipliée par http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-20.png également (les engagements de part et d’autre sont multipliés par ce même facteur).
Bon, pour les frais de gestion, c’est un peu plus compliqué. Si on reprend ce qui a été fait ici, on peut montrer que la différence entre  les VAP donne une provision de la forme suivante
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-21.png
si la durée de paiement des primes était initialement http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-24.png, et que http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png est la durée de la garantie. Si http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-24.png=http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-10.png, il n’y a alors pas de provision de gestion, le second terme étant nul. Cette provision existe dès lors que http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-23.png (ce qui arrive au final assez souvent me semble-t-il). On arrive ainsi à la provision d’inventaire, qui est souvent notée http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-22.png. Mais ce n’est pas fini, il y a aussi les frais d’acquisition. Généralement, ces frais sont occasionnés à la signature du contrat, puis amori tout au long de la durée du contrat, via un chargement (constant) sur la prime. Aussi,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-26.png
On pose alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/zill-27.png
On baisse ainsi la PM en prenant en compte les chargements pour frais d’acquisition. Cette méthode où l’on prend en compte la valeur actuelle des frais d’acquisition dans que la PM est positive s’appelle la zillmérisation.
Voilà en gros ce que j’ai compris (je suis loin d’être un expert sur le sujet).... Car pour ceux qui en douteraient encore, l’assurance-vie n’est pas du tout mon domaine de prédilection. Promis, c’est mon avant dernier billet sur les PM (le dernier étant une ébauche de correction pour le devoir maison... mais on attendra que tout le monde l’ait rendu). Maintenant s’il y en a qui ont des compléments ou des corrections à apporter, je suis preneur ! Je peux toutefois renvoyer à un billet de Cimon sur le sujet il y a quelque temps maintenant, ici. Pour les historiens, je renvoie en revanche à une traduction (en anglais) du papier d’August Zillmer présentant cette méthode, datant de 1863 () !

mardi 9 mars 2010

Comprendre la langue des actuaires

Il y a fort longtemps, Claude Bébéar expliquait que les actuaires avaient le don d’inventer des termes, d’inventer une langue qui leur est propre (cf ici pour la référence précise). Effectivement, comme me le faisaient remarquer les étudiants, c’est parfois délicat de comprendre ce que veulent les praticiens lors des entretiens pour les stages1. Fort heureusement, les mousquetaires de SéPIA ont mis en ligne sur le site de l’institut pas mal de présentations utilisées lors des formations interne de l’institut des actuaires. On retrouvera des slides sur les normes IFRS (ici ou ), ou sur la "phase 1" (ici), sur Solvabilité II (ici), sur la notion de best estimate (ici), sur les contrats en UC (ici) ou enfin sur les modèles internes (). Beaucoup d’autres devraient suivre sur le site, ici.

1 à mon époque, tous les praticiens parlaient de value-at-risk, terme très à la mode à la fin des années 90 mais dont on ne savait pas très bien ce que ça voulait dire... jusqu’à ce qu’on apprenne que la VaR n’était qu’un simple quantile....

Les actuaires ne sont pas bien exotiques...

Chose promis, chose due, une rapide étude sur les prénoms des actuaires.... Il y a environ 3660 actuaires recensés en France, i.e. membres (voire ex-membres) de l’Institut des Actuaires. Pour des raisons techniques, je n’ai pas pu récupérer leur date de naissance, mais leur année de promotion. On supposera (j’en conviens, c’est une hypothèse un peu forte) que les actuaires obtiennent leur titre vers 25 ou 30 ans. Au niveau méthodologique, j’ai retenu les actuaires sortis des formations depuis 1985, et j’ai supposé qu’ils étaient nés entre 1960 et 1985. J’ai donc deux bases de prénoms, celle des actuaires sur environ 25 ans, et celle de la population française. Si l’on regarde parmi le top 25 des prénoms donnés globalement dans la population, on notera que certains sont sortis davantage que d’autre, mais globalement on retrouve aussi des prénoms classiques chez les actuaires,

Certains prénoms sont sur-représentés (mais tout est relatif), comme les Olivier ou les Nicolas, alors que d’autres sont peu donnés, comme les Sylvie ou les Sandrine, mais je pense qu’il s’agit simplement d’un effet sexe, indiquant que les filles sont sous représentée dans la population des actuaires (de l’ordre de 30% me semble-t-il). On notera aussi les Jean ou les Pascal parmi les prénoms masculins sous-représentés (voire les Arthur, mais je commence à aller voir dans les prénoms exotiques). Les prénoms sont d’ailleurs assez proches de ceux de l’ensemble de la population, avec des prénoms de l’ordre de 7 lettres (en moyenne), comme pour l’ensemble de la population.... Pour aller un peu plus loin, on peut aussi comparer les courbes de Lorenz (ici ou pour quelques rappels), et on notera que les prénoms sont beaucoup plus dispersés chez les actuaires (membres de l’institut des actuaires français, en rouge toujours) que parmi l’ensemble de la population française, née à la même époque (en bleu).

Je n’ai malheureusement pas des données assez précises (et surtout de population assez grande) pour faire une étude aussi poussée que celle de Baptiste sur les prénoms et les CSP (ici)...

lundi 8 mars 2010

Droits des images

Un court billet de réflexion sur les droits des marques, des icônes volées, etc...

Lillie versus Kolmogorov-Smirnov

Toujours pour répondre à une question posée par mail ("ça quoi ça sert le test lillie sous R ? "), un court billet. Avant de parler de Lillie1, il faut revenir sur Kolmogorov-Smirnov. L’idée est de tester

http://perso.univ-rennes1.fr/arthur.charpentier/tex/KS-01b.png
contre
http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-02.png
Je prends une forme paramétrique pour la fonction que l’on cherche à tester car c’est ce qui est souvent fait en pratique. Et surtout, ça permettra de mieux comprendre (je pense) l’idée du test lilllie. L’idée du test repose sur le théroème de Glivenko-Cantelli. En fait, Andreï Kolmogorov et Vladimir Smirnov ont montré que l’on pouvait encadrer la vitesse de convergence,
http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-03.png
ce qui se montre en notant que http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-04.png converge en loi (si l’hypothèse H0 était la bonnne) vers un pont brownien changé de temps par la fonction quantile http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-05.png.
La loi limite est appelée loi de Kolmogorv-Smirnov, et ne dépend pas de la loi sous-jacente. Par contre elle n’est qu’asymptotique. Pour récupérer la loi à distance finie, regardons les simulations suivantes. On notera en particulier que la loi de la différence ne dépend pas de la vraie loi des observations (même si c’est cette loi que l’on teste, comme toujours lorsque l’on cherche la loi sous l’hypothèse H0)
> n=20
> ns=100000
> SU=SN=SE=rep(NA,ns)
> for(i in 1:ns){
+ X=runif(n)
+ SU[i]=max(c(abs(rank(X)/n-X),abs((rank(X)-1)/n-X)))
+ X=rnorm(n,0,1)
+ SN[i]=max(c(abs(rank(X)/n-pnorm(X,0,1)),abs((rank(X)-1)/n-pnorm(X,0,1))))
+ X=rexp(n,1)
+ SE[i]=max(c(abs(rank(X)/n-pexp(X,1)),abs((rank(X)-1)/n-pexp(X,1))))
+ }
> plot(density(SN),col="red",lwd=2)
> lines(density(SU),col="blue",lwd=2)
> lines(density(SE),col="green",lwd=2)
avec la densité de la loi du maximum entre la fonction de répartition théorique et la fonction de répartition empirique en rouge pour l’échantillon Gaussien, en bleu pour l’échantillon uniforme et en vert pour l’échantillon exponentiel. Peu importe la loi sous-jacente, la statistique de test est la même.

On peut en particulier récupérer les quantiles de cette lois pour n=20 observations (la loi du test dépend simplement du nombre d’observations).

> quantile(SE,.95)
      95%
0.2937023
> quantile(SE,.90)
     90%
0.264661
On peut d’ailleurs confirmer ces valeurs en regardant dans les tables (ici ou ), par exemple 0.294 à 95%,

Mais ce test ne sert que pour juger de l’adéquation à une loi précise (que nous avions noté http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-encore.png ici). En pratique, on ne connaît pas cette loi, mais on peut l’estimer, en particulier pour les lois paramétriques. On peut alors chercher à tester autre chose, comme
http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-06.png
contre
http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-07.png
qui s’apparente plus à un test d’ajustement à une famille de lois. Ce test étant différent du précédant, il n’y a pas de raison d’avoir la même loi pour la statistique de test (qui néanmoins, à  peut rester celle que nous avions considérée, à savoir la norme infinie entre les fonctions de répartition). Graphiquement par exemple, on cherchait aurparavant à tester l’ajustement de la loi rouge, et maintenant on va tester l’ajustement de la loi bleue,

qui va s’adapter à l’échantillon, contrairement au test initial,

Le code ici, dans le cas Gaussien par exemple (la loi sous-jacente n’est plus neutre dans ce cas) donne
> n=20
> ns=100000
> SN=SNML=rep(NA,ns)
> for(i in 1:ns){
+ X=rnorm(n,0,1)
+ SN[i]  =max(c(abs(rank(X)/n-pnorm(X,0,1)),
+               abs((rank(X)-1)/n-pnorm(X,0,1))))
+ SNML[i]=max(c(abs(rank(X)/n-pnorm(X,mean(X),sd(X))),
+               abs((rank(X)-1)/n-pnorm(X,mean(X),sd(X)))))
+ }
> plot(density(SN),col="blue",lty=2)
> lines(density(SNML),col="red",lwd=2)
Les densités sont sensiblement différentes,

et les quantiles aussi,
> quantile(SN,.95)
      95%
0.2943206
> quantile(SNML,.95)
      95%
0.1923295
Sous R, deux fonctions distinctes sont à utiliser. Par exemple l’ajustement à une loi http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-08.png donne
>  X=rnorm(20)
> mean(X)
[1] -0.553641
> sd(X)
[1] 0.8293167
>  ks.test(X,"pnorm", mean=0,sd=1)

        One-sample Kolmogorov-Smirnov test

data:  X
D = 0.3096, p-value = 0.03362
alternative hypothesis: two-sided
(on rejette ici l’ajustement à une loi normale centrée réduite) et l’ajustement à une loi http://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-09.png donne
> lillie.test(X)

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  X
D = 0.1519, p-value = 0.2629
avec une p-value de l’ordre de 26%. Notons que cette p-value est interprétable, contrairement au test suivant, qui n’est pas valide,
>  ks.test(X,"pnorm", mean=mean(X),sd=sd(X))

        One-sample Kolmogorov-Smirnov test

data:  X
D = 0.1519, p-value = 0.69
alternative hypothesis: two-sided
Encore une fois, la statistique de test est la même, sauf que la loi de la statistique de test sous H0 n’est pas juste dans le dernier cas. En particulier, n sera beaucoup plus exigent dans le test car on a choisi la meilleure loi normale possible (à l’aide d’un critère de maximisation de vraisemblance). Sinon je renvoie ici pour un document sur le même sujet.
1 Lillie pour Hubert Lilliefors qui, le premier, a proposé ce test

Quid des saints du calendrier ?

Pour poursuivre un peu mon étude sur les prénoms en France, commencée ici (oui, c’est un sujet qui me préoccupe en ce moment), on peut se demander quelle proportions des enfants qui naissent ont leur fête dans le calendrier ? Car mine de rien, tous ceux qui ont survolé un calendrier se sont un jour demandé "c’est vraiment un prénom, ça ?". Bref, si je prends le calendrier en ligne ici, retranscrit en csv , on obtient les proportions suivantes, allant de 45% en 1900 à moins de 20% en 2000,

Allons faire un tour en forêt....

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/arbre-01.PNG
Parmi les techniques de scoring et de discrimination, la méthode des arbres peut donner des pistes pour trouver un découpage "optimal" d’une variable continue en classes. Pas mal de documents sur le sujet traînent sur le net, comme ici pour le survey de Quilan (qui date un peu, mais donnne une belle introduction au sujet).
Le principe fondamental est de détercter des critères permettant de répartir une population en 2 classes (notées classiquement 0 et 1, mais que l’on peut généraliser en davantage de classes) prédéfinies. A chaque étape, on recherche la variable qui sépare le mieux les individus en chacune des classes. Une fois la meilleure séparation trouvée (suivant un critère que l’on a choisi a priori), on l’applique, puis on répète la même opération sur chaque noeud de manière à augmenter la discrimination. A chaque noeud, on crée ce qui s’appelle des "noeuds-fils". Chaque noeud-fils donne à son tour naissance à un ou plusieurs noeuds, etc. On arrête dès que les noeuds ne contiennent qu’une unique observation (ce qui est une classification un peu extrême) ou en se fixant un critère d’arrêt.
  • Un peu de théorie
Pour obtennir la meilleur séparation possible, on peut utiliser
  • un critère de type chi-deux,  lorsque les variables explicatives sont qualitatives, ou discrètes. On parle d’arbre CHAID, pour CHi-square Automatic Interaction Detection (ici par exemple)
  • le critère de Gini, pour toutes sortes de variables explicatives. On parle d’arbre CART, pour Classification And Regression Tree (ici)
  • l’entropie, là aussi pour toutes sortes de variables. Ce sont les arbres C4.5 ou C5, de Quinlan (ici ou )
Par exemple pour l’indice de Gini à un noeud, on calcule
http://perso.univ-rennes1.fr/arthur.charpentier/latex/tree-01.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/tree-02.png correspond à la fréquence dans chacun des noeuds dans les classes (ici i=0,1). Aussi, si la variable est uniformément répartie dans les deux classes, l’indice de Gini sera élevé. On cherche alors à minimiser l’indice de Gini. Avec deux classes, l’indice est compris entre 0 et 1/2. Notons que l’indice de Gini mesure la probabilité que deux individus, tirés au hasard dans un noeud, appartiennent à deux classes différentes.
La séparation doit causer la plus grande baisse possible de l’indice de Gini. Dans le cas de l’entropie, on calcule
http://perso.univ-rennes1.fr/arthur.charpentier/latex/tree-03.png
La technique CART a été inventée en 1984 par Breiman, Friedman, Olshen et Stone. Comme on le voit sur le dessin ci-dessous, un des défauts de cette technique est d’obliger une rectangularisation des classes. On dipose de deux classes (les points rouges et les points verts), que l’on souhaite discriminer à l’aide de deux variables explicatives continues.
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.arbre-rpart-simulation_m.jpg
Le premier noeud à être constitué consiste à séparer l’espace en deux, suivant un axe vectical, i.e. les points tels que X1<1.77, et les points tels que X1>1.77. Dans un second temps, on discrimine pour les points à droite de l’axe en regardant la valeur de X2 (par rapport à 0.24). Pour la partie de gauche, on sépare à nouveau en deux verticalement, etc. Dans les arbres CHAID, proposés par Kass en 1980 à partir d’une idée de Morgan et Sonquist (datant de 1963) dans on remplace le test du chi-deux par un test de Fisher (d’analyse de la variance). On choisit les noeuds fils de manière à minimiser la variance intraclasse, et à maximiser la variance intraclasse. Classiquement, on scinde un noeud dès que le ratio variance intraclasse / variance interclasse atteint le seuil de 20%.
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/arbre-10.PNG
  • Application aux accidents cardiaques
Pour illustrer un peu mieux l’impact du critère de choix, dans la méthode CART, regardons cet exemple dont les données sont tirés d’un exemple emprunté à Gilbert Saporta),
> MYOCARDE=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/saporta.csv",head=TRUE,sep=";")
> head(MYOCARDE)
  FRCAR INCAR INSYS PRDIA PAPUL PVENT REPUL  PRONO
1    90  1.71  19.0    16  19.5  16.0   912 SURVIE
2    90  1.68  18.7    24  31.0  14.0  1476  DECES
3   120  1.40  11.7    23  29.0   8.0  1657  DECES
4    82  1.79  21.8    14  17.5  10.0   782 SURVIE
5    80  1.58  19.7    21  28.0  18.5  1418  DECES
6    80  1.13  14.1    18  23.5   9.0  1664  DECES

Il s’agit de victimes d’infarctus du myocarde, qui ont été observés à leur admission aux urgences, avec la fréquence cardiaque (FRCAR), un indcex cardiaque(INCAR), index systolique (INSYS), pression diastolique (PRDIA), pression arterielle pulmonaire (PAPUL), pression venticulaire (PVENT) et resistance pulmonaire (REPUL). Bref, à partir de toutes ces variables on souhaite mieux diagnostiquer les décès. Si on fait un peu de statistique descriptive, on obtient
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.myocarde-boxplot-1_m.jpg

et
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.myocarde-boxplot-2_m.jpg
On peut bien sûr faire une régression logisitque, ou une analyse discriminante, mais le but est ici de créer des catégories, i.e. de discrétiser de manière optimale ces classes.
Dans la méthode CART, on regarde l’ensemble des séparations possibles. Par exemple sur le graphique ci-dessous, on observe l’entropie obtenue en découpant à différents niveaux,
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.myocarde-arbr-plit-entropie_m.jpg
Le découpage optimale est obtenu en cherchant le minimum de l’entropie, i.e. lorsque l’index systolique est inférieur (ou supérieur) à 19. Notons qu’avec un critère de type Gini, on obtient quelque chose d’assez proche,
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.myocarde-arbr-plit-gini_m.jpg
où le découpage optimal est à nouveau basé sur un découpage de l’index systolique aux alentours de 19. A l’étape suivante, on sépare effectivement suivant ce critère, et on refait tourner l’algorithme sur chacune des sous-classes.
> tree(PRONO~FRCAR+INCAR+INSYS+PRDIA+PAPUL+PVENT+REPUL,data=MYOCARDE,split="gini")
node), split, n, deviance, yval, (yprob)
      * denotes terminal node

 1) root 71 96.030 SURVIE ( 0.40845 0.59155 ) 
   2) INSYS < 18.85 27 18.840 DECES ( 0.88889 0.11111 ) 
     4) PVENT < 15.25 22  8.136 DECES ( 0.95455 0.04545 ) 
       8) INCAR < 1.58 17  0.000 DECES ( 1.00000 0.00000 ) *
       9) INCAR > 1.58 5  5.004 DECES ( 0.80000 0.20000 ) *
     5) PVENT > 15.25 5  6.730 DECES ( 0.60000 0.40000 ) *
   3) INSYS > 18.85 44 31.160 SURVIE ( 0.11364 0.88636 ) 
     6) REPUL < 1094.5 37  9.195 SURVIE ( 0.02703 0.97297 ) 
      12) PVENT < 13 32  0.000 SURVIE ( 0.00000 1.00000 ) *
      13) PVENT > 13 5  5.004 SURVIE ( 0.20000 0.80000 ) *
     7) REPUL > 1094.5 7  9.561 DECES ( 0.57143 0.42857 ) *
On retrouve que le premier critère est effectivement suivant l’index systolique. Si l’on étude la population des indices elevés, on trouve que le critère optimal de partition est basé sur la resistance pulmonaire
> tree(PRONO~FRCAR+INCAR+INSYS+PRDIA+PAPUL+PVENT+REPUL,data=MYOCARDE[MYOCARDE$INSYS>18.85,],split="gini")
node), split, n, deviance, yval, (yprob)
      * denotes terminal node

1) root 44 31.160 SURVIE ( 0.11364 0.88636 ) 
  2) REPUL < 1094.5 37  9.195 SURVIE ( 0.02703 0.97297 ) 
    4) PVENT < 13 32  0.000 SURVIE ( 0.00000 1.00000 ) *
    5) PVENT > 13 5  5.004 SURVIE ( 0.20000 0.80000 ) *
  3) REPUL > 1094.5 7  9.561 DECES ( 0.57143 0.42857 ) *

...etc. Tout cela se résume dans l’arbre de décision suivant

Cette technique peut être utilisée pour découper en classes de manière intelligente (comme j’avais essayé de le faire ici), à condition que l’algorithme tourne. En effet, il n’est pas rare qu’aucune classe ne soit trouvé....

dimanche 7 mars 2010

De l'originalité dans les prénoms des français(e)s...

Lors d’une discussion avec quelques anciens collègues samedi midi, nous évoquions des prénoms, et du fait que les prénoms rares (pour ne pas dire bizarres) étaient de plus en plus fréquents. Je me suis demandé si les gens donnaient vraiment plus des prénoms originaux qu’avant. Était-ce une idée en l’air (du buzz médiatique, genre ici ou ), ou au contraire est-elle (statistiquement) fondée ?


Pour mesurer la concentration des prénoms (ou de n’importe quoi d’ailleurs), le plus simple est peut être d’utiliser la courbe de Lorenz. En statistique descriptive, on construit cette courbe de la manière suivante: à partir d’un échantillon que l’on supposera ordonné, i.e. , la courbe de Lorenz est le nuage de points et


La relecture probabiliste est que l’on trace


Je renvoie ici ou pour plus d’information sur la courbe de Lorenz.
Si je reviens à mon histoire de prénoms, on pourra noter qu’en 2003, par exemple1, il y a eu 795493 naissances, donnant lieu à 59481 prénoms (j’entends par à des prénoms d’orthographe différente),
> head(base)
    sexe preusuel annais nombre rang
32     1    AADIL   2003      3 8775
71     1    AARON   2003    135  764
100    1    ABASS   2003      3 8301
133    1    ABBAS   2003      5 6897
171    1    ABBES   2003      3 8816
228    1      ABD   2003     11 3982

Si on trie, on obtient le classement suivant,
> head(baser)
    sexe preusuel annais nombre rang
300101    1      LEA   2003   8988    1
110747    1    LUCAS   2003   8292   2
168477    1     THEO   2003   7855    3
77273    1     HUGO   2003   7386    4
318115    1    MANON   2003   6917    5

Autrement dit, le top 5 des prénoms représente 5% des naissances, et les 100 premiers représentent 43% des naissances. En 1900, les 5 premiers représentaient 22% de l’ensemble, et 80% pour les 100 premiers. Afin de résumer l’information, on notera, par exemple, que 2,5% des prénoms (les prénoms les plus donnés) représentaient 50% des naissances en 2003. C’est ceci que l’on récupère via la courbe de Lorenz.
Le graphique ci-dessous montre combien représentaient 10% des prénoms (les prénoms le plus donnés, en rouge) au cours du temps, et 2,5% des prénoms (en bleu). Autrefois, si l’on prenait les 2,5% des prénoms les plus donnés, ils étaient portés par 9O% de la population en 1900 (100-10), contre 50% de la population en 2004 (100-50) [des compléments sur la lecture des graphiques ont été apportés ici]

On peut alors effectivement parler de "course à l’originalité" depuis 1970. En 1900, seuls 5% des naissances donnaient lieu à des prénoms "originaux", contre 50% aujourd’hui. Étrangement, on notera qu’il y a eu des "pics d’originalité" pendant les guerres. Si quelqu’un a une explication ?... On peut aussi visualiser ci-dessous les courbes de Lorenz,

Les courbes montrent que les prénoms sont de moins en moins "concentrés" (pour reprendre une terminologie classique sur les courbes de Lorenz). On le retrouve également sur les prénoms à la mode qui furent très stables par le passé, et qui depuis les années 70 n’arrêtent pas de changer... Bref, il y a de plus en plus de prénoms originaux, certes, mais les prénoms classiques ou à la mode restent donnés par beaucoup de monde, l’indice de Gini restant très élevé. Et je ne suis pas assez calé en statistique textuelle pour voir si ce qui change, c’est le nombre de prénoms (avec un total de 6556 prénoms en 1900 contre 10 fois plus en 2004) ou simplement l’orthographe,
       sexe       preusuel annais nombre
 313765    1       MAEL   2004     37
313861    1      MAELE   2004     27
314056    1      MAELL   2004      3
314127    1     MAELLE   2004   1577
315525    1      MAHEL   2004      3

etc....
La prochaine fois, j’essayerais de reprendre une étude intéressante sur les liens entre les prénoms et les professions (ici) pour étudier les prénoms des actuaires.... oui, j’ai toute la base ! à suivre donc...

1 Je passe sur les soucis techniques de programmation. En effet, il existe un libellé étrange pour les "prénoms rares",
       sexe       preusuel annais nombre 
408961    1 _PRENOMS_RARES   2003  26253    
185715    1 _PRENOMS_RARES   2003  23718    

mais comme on a le tail on peut construire la courbe de Lorenz proprement... enfin, au moins par la partie de droite car j’ai postulé car les prénoms rares étaient donnés une unique fois (mais a priori c’est une ou deux, et je n’ai pas moyen de le savoir.... je n’ai que le nombre de naissances avec un prénom rare). Mais ça n’influence que la partie gauche qui peut être un peu différente....

vendredi 5 mars 2010

Effets nonlinéaires en tarification

Avant la formation SéPIA sur la tarification (ici), un petit mot pour répondre aux questions des étudiants sur la prise en compte d’effets nonlinéaires en tarification. Pour information complémentaire, la formation commencera par une introduction de Jean Pierre Verlé sur la tarification puis par un exposé de Michel Grun-Rehomme sur les méthodes zero-inflated, et sur l’écrètement. J’interviendrai par la suite précisément sur les méthodes pour prendre en compte les effets nonlinéaires.
Avant de revenir sur la prise en compte de variables continues discriminantes, un extrait de Popular Sciences, datant de 1968 (ici),

Attaquons la formalisation. Je ne mettrais pas le lien vers la base, mais je peux mettre le code utilisé par la suite. Considérons comme variable d’intérêt le fait d’avoir eu, ou pas d’accident dans l’année. La variable explicative est ici l’âge du conducteur.

Sur le graphique ci-dessus, on visualise en points noirs la moyenne empirique pour chaque âge. On fait également une régression logistique, qui nous dit que l’âge n’est pas une variable significative
> RGLM=glm(Y~age,data=B,family="binomial")
> summary(RGLM)
Call:
glm(formula = Y ~ age, family = "binomial", data = B)
Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-0.4799  -0.4686  -0.4599  -0.4513   2.1917 
Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept) -2.058649   0.092296 -22.305   <2e-16 ***
age         -0.002640   0.001812  -1.457    0.145   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 6809.7  on 10398  degrees of freedom
Residual deviance: 6807.6  on 10397  degrees of freedom
AIC: 6811.6
Number of Fisher Scoring iterations: 4
En fait, l’âge n’est pas significatif pris tel quel, i.e. en tant que variable continu, et de manière linéaire. Visiblement, on peut faire mieux. Tentons alors un modèle GAM,

ou en changeant les paramètres de lissage,

On repère ici un effet non linéaire ! Pour constituer des classes, on peut alors tenter d’utiliser les arbres de régression (plutôt que de tenter un découpage manuel)
> library(tree)
> TREEg=tree(Y~age,data=B)
>  plot(TREEg)
>  text(TREEg,col="red")
> TREEg=tree(Y~age,data=B,split="gini")
>  plot(TREEg)
>  text(TREEg,col="blue")

Je ferais bientôt un billet sur les arbres, mais en attendant, je peux renvoyer ici ou . On peut aussi utiliser une autre librairie de fonctions, qui donne sensiblement les mêmes choses
> library(rpart)
> fit = rpart(Y ~ age, method="anova", data=B)
> summary(fit) 
Call:
rpart(formula = Y ~ age, data = B, method = "anova")
  n= 10399
          CP nsplit rel error    xerror       xstd
1 0.01748133      0 1.0000000 1.0001955 0.02596249
2 0.01000000      2 0.9650373 0.9687768 0.02468425
Node number 1: 10399 observations,    complexity param=0.01748133
  mean=0.1010674, MSE=0.09085279
  left son=2 (8873 obs) right son=3 (1526 obs)
  Primary splits:
      age < 27.5 to the right, improve=0.01588077, (0 missing)
Node number 2: 8873 observations,    complexity param=0.01748133
  mean=0.085315, MSE=0.07803635
  left son=4 (7935 obs) right son=5 (938 obs)
  Primary splits:
      age < 74.5 to the left,  improve=0.02603656, (0 missing)
Node number 3: 1526 observations
  mean=0.1926606, MSE=0.1555425
Node number 4: 7935 observations
  mean=0.06981727, MSE=0.06494281
Node number 5: 938 observations
  mean=0.2164179, MSE=0.1695812
On peut alors utiliser les classes proposées ici pour créer les classes tarifaires, à savoir
[18-27] - [28-74] - [74-100]
La fonction sous R est ici
> brk=c(0,27.5,74.5,120)
> cut(X[1:20],brk)
 [1] (27.5,74.5] (27.5,74.5] (0,27.5]    (27.5,74.5] (27.5,74.5] (27.5,74.5]
 [7] (27.5,74.5] (74.5,120]  (27.5,74.5] (74.5,120]  (27.5,74.5] (27.5,74.5]
[13] (27.5,74.5] (27.5,74.5] (27.5,74.5] (27.5,74.5] (27.5,74.5] (27.5,74.5]
[19] (27.5,74.5] (27.5,74.5]
Levels: (0,27.5] (27.5,74.5] (74.5,120]
On peut utiliser cette fonction directement dans la régression logistique, et faire ensuite une prédiction
> RGLMb=glm(Y~cut(age,brk),data=B,family="binomial")
> summary(RGLMb)
Call:
glm(formula = Y ~ cut(age, brk), family = "binomial", data = B)
Coefficients:
                         Estimate Std. Error z value Pr(>|z|)   
(Intercept)              -1.43281    0.06491 -22.075   <2e-16 ***
cut(age, brk)(27.5,74.5] -1.15669    0.07844 -14.745   <2e-16 ***
cut(age, brk)(74.5,120]   0.14615    0.10247   1.426    0.154   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 6809.7  on 10398  degrees of freedom
Residual deviance: 6493.3  on 10396  degrees of freedom
AIC: 6499.3
Number of Fisher Scoring iterations: 5
On notera qu’ici, on nous suggère de regrouper la première classe (ici considérée comme classe de référence) et la dernière (qui n’est pas significativement différente de la classe de référence). La prédicition donne

ce qui n’est pas complètement absurde... On peut essayer d’aller un peu plus loin, en creusant davantage dans les classes construites dans la construction des artbres. Par exemple,
> brk=c(15,20.5,24.5,32.5,47.5,52.5,74.5,120)
> RGLMb=glm(Y~cut(age,brk),data=B,family="binomial")
> summary(RGLMb)
Call:
glm(formula = Y ~ cut(age, brk), family = "binomial", data = B)
Coefficients:
                         Estimate Std. Error z value Pr(>|z|)   
(Intercept)               -1.0515     0.1079  -9.750  < 2e-16 ***
cut(age, brk)(20.5,24.5]  -0.5071     0.1537  -3.300 0.000967 ***
cut(age, brk)(24.5,32.5]  -0.9681     0.1353  -7.157 8.25e-13 ***
cut(age, brk)(32.5,47.5]  -1.7104     0.1353 -12.637  < 2e-16 ***
cut(age, brk)(47.5,52.5]  -1.1391     0.1505  -7.570 3.73e-14 ***
cut(age, brk)(52.5,74.5]  -1.6843     0.1301 -12.945  < 2e-16 ***
cut(age, brk)(74.5,120]   -0.2351     0.1339  -1.756 0.079017 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 6809.7  on 10398  degrees of freedom
Residual deviance: 6453.5  on 10392  degrees of freedom
AIC: 6467.5
Number of Fisher Scoring iterations: 5

Après, il convient bien entendu de faire éventuellement des regroupements de classes, via plusieurs tests de Fisher...
En revanche, sur des bases plus volumineuse, il est délicat de trouver la partition à retenir, car il faut arbitrairement couper dans l’arbre (surtout que la base précédante était simulée). Sur des vraies données (plus de 100000 contrats), on obtient

La classification est plus complexe à lire


et ça sera l’objet d’un autre billet....

Copules et dépendance (2)

Je poursuis le billet que j’avais commencé il y a quelque temps (ici) sur les copules. En effet, un étudiant me dit avoir utilisé "l’algorithme de Devroye" qui est donné comme exercice dans le livre de Roger Nelsen,

  • Un algorithme de génération de la copule de Clayton
Je n’aime pas faire les exos à la place des élèves (surtout quand ce ne sont pas les miens), il faut faire un peu de calculs. En fait, c’est comme ça que Cook et Johson ont obtenu, en 1981, une distribution qui correspondait avec celle obtenue par Clayton.

On notera que ce n’est pas l’utilisation de la représentation frailty de la copule de Clayton...
  • Travail sur un algorithme faux, et tests d’ajustement
Malheureusement, dans l’algorithme envoyé, l’étudiant avait permutté les facteurs de scale et de shape de la loi Gamma, ce qui donne des choses assez différentes. Mais sur lesquelles je vais rebondir. On va faire comme si on n’avait pas vu la faute dans l’algorithme, et que l’on se demande si on a généré la copule de Clayton que l’on pense avoir simulée....
Considérons l’algorithme suivant, qui génére effectivment une distribution à marges dans [0,1], on se demande si la copule associé est, ou pas, la copule de Clayton (donc on prend les rangs)
> n=10000
> x1=rexp(n,1)
> x2=rexp(n,1)
> theta=1/2
> x=rgamma(n, shape=1, scale=theta)
> u1=(1+x1/x)^(-theta)
> u2=(1+x2/x)^(-theta)
> x1=rank(u1)/(n+1)
> x2=rank(u2)/(n+1)
> plot(x1,x2)
Vu de loin, cela ressemble à une copule de Clayton.

On peut d’ailleurs essayer de récupérer la densité de la copule,
>  beta.kernel.copula.surface <- function (u,v,bx,by,p) {
+  s <- seq(1/p, len=(p-1), by=1/p)
+  mat <- matrix(0,nrow = p-1, ncol = p-1)
+ for (i in 1:(p-1)) {
+    a <- s[i]
+    for (j in 1:(p-1)) {
+        b <- s[j]
+    mat[i,j] <- sum(dbeta(a,u/bx,(1-u)/bx) * dbeta(b,v/by,(1-v)/by)) / length(u)
+    } }
+ return(data.matrix(mat)) }
> Z= beta.kernel.copula.surface(x1,x2,bx=.1,by=.1,p=26)
>  p=26
> u=seq(1/p, len=(p-1), by=1/p)
> persp(u,u,Z,theta=30,col="green",shade=TRUE)

ou en changeant la fenêtre de l’estimation à noyau
> Z= beta.kernel.copula.surface(x1,x2,bx=.01,by=.01,p=26)

Bref, ça ressemble à du Clayton.... Si on regarde la section diagonale de la copule empirique, on note que l’on n’obtiens ni la copule de Clayton de paramètre 2, ni celle de paramètre 1/2,

> U=seq(.01,1,by=.01)
> D=rep(NA,length(U))
> for(i in 1:length(U)){
+ D[i]=mean((x1<U[i])&(x2<U[i]))
+ }
> U=c(0,U)
> plot(U,c(0,D),cex=.5)
> lines(U,(2*U^(-1/theta)-1)^(-theta),col="red")
> lines(U,(2*U^(-theta)-1)^(-1/theta),col="blue")
Bref, ça ne semble pas correspondre à une des deux copules de Clayton recherchées,

On peut aussi des grandeurs comme la fonction
http://perso.univ-rennes1.fr/arthur.charpentier/latex/genest-1.png
qui est l’extension en dimension 2 de l’integral probability transform (cf ici par exemple). Ou plutôt la fonction http://perso.univ-rennes1.fr/arthur.charpentier/latex/genest-3.png définie comme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/genest-2.png
> V=rep(NA,n)
> for(i in 1:n){
+ V[i]=sum((x1<x1[i])&(x2<x2[i]))/(n-1)
+ }
> lambda=sort(V)-(1:n)/(n+1)
> plot(sort(V),lambda)
> u=seq(0,1,by=0.01)
> l=(u^(-theta)-1)/(-theta*u^(-theta-1))
> lines(u,l,col="blue")
> l=(u^(-1/theta)-1)/(-1/theta*u^(-1/theta-1))
> lines(u,l,col="red")

Pour conclure définitivement, je reprendrais des choses que j’avais faites dans ma thèse. En fait, la copule de Clayton est la seule copule continue qui soit invariante par troncature. On peut alors regarder les tau de Kendall ou Rho de Spearman conditionnels. Normalement, ils devraient être constants,
> U=seq(.01,1,by=.01)
> R=K=rep(NA,length(U))
> for(i in 1:length(U)){
+ I=(x1<U[i])&(x2<U[i])
+ R[i]=cor(x1[I],x2[I],method = "spearman")
+ K[i]=cor(x1[I],x2[I],method = "kendall")
+ }
> plot(U,R,cex=.75,type="b")
> abline(h=R[length(R)],col="green")

> plot(U,K,cex=.75,type="b")
> abline(h=K[length(R)],col="purple")

Moralité (même si je n’ai pas tracé les régions de confiance des tests) je pense pouvoir conclure que l’algorithme tel qu’il est là ne génêre pas une copule de Clayton de paramètre 2 ou 1/2.

Ecart absolue ou écart type ?

Un petit billet très court pour répondre à une question d’élèves de licence. L’erreur quadratique moyen est définie comme

http://perso.univ-rennes1.fr/arthur.charpentier/latex/eqt-01.png
alors que l’erreur absolue moyenne est définie comme
http://perso.univ-rennes1.fr/arthur.charpentier/latex/eqt-02.png
(que je définie ici comme l’écart à la moyenne, certain considérant la distance à la médiane, qui minimise cette quantitié, alors que la moyenne minimise l’écart quadratique, mais je reviendrais là dessus bientôt, sur l’histoire de la loi normale).
Souvent ces quantités sont présentées comme des mesures de la même quantité, que l’on appelerait "dispersion" autour de la moyenne.
Si l’on prend deux observations, ces quantités sont rigoureusement identiques, car
http://perso.univ-rennes1.fr/arthur.charpentier/latex/eqt-03.png
compte tenu du fait que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/eqt-04.png
En revanche, comme le montre le graphique ci-dessous,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/eqt-00-col.png

Aussi, avec davantage de termes, la différence peut devenir significativement différente. En particulier, comme le notent Goldstein et Taleb (2007), en ligne ici, dans le cas Gaussien centré réduit,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/eqt-05.png
Autrement dit, ces quantités peuvent parfaitement être différentes. Sur l’exemple ci-dessous, on considère un groupe de 4 individus, trois ayant le même revenu, et le dernier ayant un revenu plus important. On s’arrange pour que le revenu moyen soit constant (100 ici) et on regarde comment se comportent ces deux quantités lorsque le revenu du dernier individu s’accroit (de celui des autres).

(j’ai autorisé les revenus négatifs pour la simplicité des calculs et pour avoir un joli dessin). On en conclue qu’un écart significativement différent entre l’écart-type (L2) et l’écart absolu (L1) signifie simplement qu’il y a probablement des points aberrants dans la base (aberrants au sens "sensiblement différents", des "outliers").

- page 1 de 40