Category Archives: Computer

Tails of copulas, une lecture graphique

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

  • Quantifier la dépendance de queue

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

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/toc3latex2png.2.php.png

soit

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/toc3latex2png.3.php.png

La seconde est liée à une idée que l’on retrouve dans les travaux de Janet Heffernan, Stuart Coles ou Jonathan Tawn. L’intuition est la suivante (on peut la retrouver en ligne ici). Si https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-2.2.php.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-3.2.php.png ont la même loi et que l’on suppose les variables indépendantes, alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-1.2.php.png

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

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-4.2.php.png

Aussi, on peut supposer qu’il existe un indice https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png tel que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-5.2.php.png

Le soucis est que le cas d’indépendance correspond à https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png=2, alors que le cas de dépendance forte correspond au cas https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png=1. Il est alors usuel de faire une transformation affine pour se ramener sur [0,1], et que la force de la dépendance soit croissante avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png, e.g.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-8.2.php.png

Posons alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toc2latex2png.2.php.png

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

  • Les fonctions de concentration dans les queues

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

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/Llatex2png.2.php.png

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

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/Clatex2png.2.php.png

pour la queue supérieure,où https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-12.2.php.png est la copule de survie associée à https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-13.2.php.png, au sens où
https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-14.2.php.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-15.2.php.png

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

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toc2latex2png.3.php.png

ou

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toc2latex2png.4.php.png
  • Application statistique

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

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

http://perso.univ-re<br /><br /> nnes1.fr/arthur.charpentier/latex/toclatex2png-18.2.php.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-19.2.php.png

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

Le code R ressemble à ça

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

avec le code suivant pour la version empirique,

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

et pour la version théorique,

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

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

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

et celui là pour la version empirique

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

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

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

On obtient ici

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

et enfin pour une copule de Gumbel,

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

ou pour une copule de Clayton,

alors qu’une copule de Gumbel donnerait

De la vitesse de convergence, partie 1

Samedi après midi, à l’heure de la sieste pour certains, nous discutions avec Magali, Jean Baptiste et Florent des machines qui pouvaient trier les cartes, de là, la discussion a dévié sur les algorithmes de tri, et sur la difficulté d’expliquer l’importance de la vitesse d’un algorithme, ou plus généralement, la vitesse de convergence. Bref, avant de parler de vitesse de convergence en statistique – ou peut être avant sur les méthodes de calcul d’intégrale (ou d’espérance) – je voulais faire un petit billet rapide sur les algorithmes de tri. Je commence par des choses simples. De mémoire c’est des choses que je me souviens avoir vu au lycée dans les années 90, dans les cours d’informatique… Mais l’exemple peut être intéressant pour comprendre la différence, par exemple, entre https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-02.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-03.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-04.png.
Le premier algorithme est appelé tri à bulles. C’est probablement l’algorithme le plus simple, et sa vitesse est en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-01.png. On va balayer autant de fois la base que l’on a d’observations, et à chaque passage, on récupère le maximum (ou le minimum, peu importe). L’algorithme ressemble à ça,

> n=25; s=1
> X=sample(1:n,size=n,replace=FALSE)
>          for(i in 1:(n-1)){
+          for(j in 1:(n-1)){
+          s=s+1
+          if(X[j] > X[j + 1]){X[c(j,j+1)]=X[c(j+1,j)]}
+          }}

Si on l’implémente, voilà ce que ça donne (tous les graphiques partent de la même permutation initiale, et tous contiennent autant d’images, de l’ordre de https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-01.png)

Bref, un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png.Bon, j’ai été un peu brutal dans mon algorithme. En fait, on voit que les derniers points sont rapidement bien classés, donc inutile que la seconde boucle aille jusqu’au bout… on peut donc faire du https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-06.png, en faisant

 

+          for(j in 1:(n-i)){

Bref, la constante dans les https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png peut être relativement importante comme on le voit ici.

Mais on peut faire  mieux, voire beaucoup mieux comme on le verra ensuite.
Le second algorithme est le tri par insertion. C’est aussi un algorithme très simple:  on balaye l’ensemble en plaçant les valeurs rencontrées à leur place, un peu comme quand on tri ses cartes… L’algorithme est là aussi en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png. Ici l’algorithme converge un peu plus rapidement que le précédant, mais cela vient du fait que les données étaient partiellement ordonnées sur la fin. Autrement dit, un autre échantillon aurait pu donner autre chose… Mais généralement, c’est le cas.

>          for(i in 1:n){
+           x = X[i]
+          for(j in (i+1):2){
+          s=s+1
+         if(X[j - 1] < x){X[j] = X[j - 1];  j = j - 1;  X[j] = x}
+         }
+         }

(le tri est fait à l’envers ici…. mais c’est un détail).

Enfin, le tri fusion est un autre algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-07.png. L’idée est de couper, de trier, puis de fusionner. L’idée étant que trier un tableau dont la première moitié des éléments sont ordonnés, et où la dernière moitié est ordonné doit être plus rapide que quand tout est dans le désordre. A chaque étape, l’algorithme est simple, on compare le minimum des deux listes, puis on le vire à l’étape suivante. Sous R, ça donne, quand on commence

 

> for(i in 0:(trunc(n/2)-1)){
+ Y=X[c(1,2)+2*i]
+ if(Y[1]>Y[2]){X[c(1,2)+2*i]=X[c(2,1)+2*i]}
+ }

Bref, en commençant par classer les paires (1,2) puis (3,4), etc, puis les listes (1,2) versus (3,4). Le code est alors de la forme

> k=4
> for(i in 0:(trunc(n/4)-1)){
+ Y=X[1:k+k*i]
+ Y1=c(Y[1:(k/2)],rep(999999,k/2));Y2=c(Y[(k/2+1):k],rep(999999,k/2));
+ for(j in 1:k+k*i){
+     X[j]=min(c(Y1[1],Y2[1]),na.rm=TRUE)
+    if(X[j]==Y1[1]){Y1=Y1[
-1]}
+    if(X[j]==Y2[1]){Y2=Y2[-1]}
+ }}

(avec une petite correction pour le dernier bloc, qui ne compare pas forcément un bloc de taille 2 contre un autre de taille 2 (mais l’adaptation du code est infime). Puis on continue avec (1,2,3,4) versus (5,6,7,8), c’est relativement rapide de tout classer, comme le montre l’animation ci-dessous.

Il existe de très nombreux autres algorithmes de tri, le but étant de monter qu’un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-07.png va effectivement beaucoup plus vite qu’un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png. J’ai tout fait ici en R, mais pour les adeptes d’autres langages, il existe plein de code en ligne, ici par exemple. Pour conclure, le graphique ci-dessous montre les différences entre les vitesses de convergence en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-08.png (en noir), en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-09.png (en rouge) et en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-10..png (en bleu),

La suite sur la convergence des méthodes de monte carlo, puis des théorèmes limites en stats….

Acheter des pochettes Pokemon ? et puis quoi encore !

Il y a quelque temps, Olivier Gergaud était venu à Rennes présenter son papier Stardom and Perceived Talent.Nous en avions profité pour discuter un peu d’un ancien papier qu’il avait fait sur les Pokemon (ici). Il faut avouer que je suis un peu jaloux car il essaye de répondre à une question que je me posais: le prix des cartes Pokemon a-t-il un lien avec les “pouvoirs” de la carte ou sa rareté.

Mais aujourd’hui je voulais répondre à une question simple que je me pose depuis que mon fils me dit (sans grand succès) “Papa, achète moi des cartes Pokemon” (ou “Naruto“, ou je ne sais quoi encore). Mon argument, en tant que père rationnel, est que ces cartes couchent très chers. Aujourd’hui, je vais le prouver. Supposons qu’il existe 200 cartes distinctes, et que ces cartes sont vendus 5 euros le paquet de 10. Pour simplifier, supposons que ces cartes soient vendues individuellement dans des cachés hermétiques à 50 centimes pièce. On va supposer que les cartes soient équi-réparties. Supposons que je puisse acheter une vignette spécifique (par exemple en ligne), dont je connais le personnage, pour 1 euro (soit le double du prix si le personnage est inconnu à l’achat). Quelle stratégie doit-je adopter ? et combien ça va me coûter ?
Je note https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-17.png le nombre de carte total (ici 200 car je suppose qu’il existe 200 personnages). Je note https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-01.png le coût d’une carte cachée (50 centimes ici) et https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-02.png le coût d’une carte non-cachée (1 euro). Supposons que j’ai à en ma possession les cartes https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-03.png (et que je n’ai pas les suivantes). Si je veux récupérer une nouvelle carte (que je noterais https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-05.png), en achetant au hasard i.e. dans des sachets hermétiques, le nombre de cartes à acheter est https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-06.png, qui suit une loi géométrique, de paramètre https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-07.png, où

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-08.png

Par exemple, si j’ai https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-09.pngcartes, et qu’il m’en manque simplement une, alors la probabilité de la tirer est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-10.png

Pour revenir au tirage d’une nouvelle carte que je n’avais pas (alors que j’en https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-04.png), le coût moyen de ma stratégie est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-11.png

J’ai donc intérêt à acheter cette carte directement dès lors que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-12.png

ce qui devient intéressant dès lors que

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-13.png

Dans mon exemple numérique, https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-02.png étant le double de https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-01.png, je commence à acheter directement dès lors que j’ai en ma possession 100 cartes distinctes. Le coût moyen est alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-14.png

Aussi, dans le cas que j’évoquais, on obtient 171€,
> .5*sum(200/(200-(1:101)-1))+1*(200-101+1)
[1] 171.0933

Je peux comparer à deux stratégies que l’on connaît. La première est de tout acheter à un collectionneur (i.e. à découvert, ce qui limite l’intérêt du “jeu“, j’en conviens). Dans ce cas, le coût total sera

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-21.png

soit 200€. L’autre stratégie est d’acheter continuellement des paquets clos chez mon buraliste. Dans ce cas, le coût moyen serait

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-16.png

Le calcul exact donne ici 587,80€, comme  le montre le calcul ci-dessous,
> 0.5*200*sum(1/1:200)
[1] 587.8031

mais on peut aussi utiliser une approximation de la forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-22.png

(ce qui donnerait 529€ environ). Pour ceux qui ont besoin de voir pour croire, j’ai représenté ci-dessous le prix en fonction du moment où je décide de passer des cartes fermées aux cartes visibles.

Ça y est, j’ai enfin un argument rationnel à faire valoir à mon fils…. même si je ne suis pas certain de réussir à le convaincre aussi facilement…

Exchange rates and correlation matrices

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.euro_dollar_yen_s.jpg I wanted to upload here a small problem I started to work on…. unfortunatley, I could not find (yet) a proper answer. Any comments and suggestions are welcomed. The problem is simple: consider 3 random variables that can be interpreted as exchange rates. Thus X is the exchange rate USD versus EUR, Y is the exchange rate EUR versus YEN, while Z is the exchange rate USD versus YEN. We then have the following constraint XY=Z. Given that constraint, what are the implied constraints on the correlation matrix of random vector (X,Y,Z) ?

  • Bounds on correlation

There are classically standard constraints on correlation. For instance, from Cauchy-Schwarz’s inequality

https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-01.png

but those bounds are not necessarily sharp, given distributions of https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-02.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-03.png. From Fréchet-Hoeffding bounds, we know that

https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-04.png

where https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-05.png. Further, correlation matrices are necessarily positive-semidefinite matrices.
To go further to find possible bounds, recall that

https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-06.png

and therefore

https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-07.png

i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-08.png

Similarity, https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-09.png will be a function of https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-10.png.
Thus, given marginal distributions of https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-02.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-03.png, covariances (and also correlations https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-11.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-12.png) are functions of L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-13.png” ne peut être affichée car elle contient des erreurs.,https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-14.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-15.png.
Somehow, we have obtained 2 degrees of freedom: if there are no additional constraints on those two coefficients, there should be no other constraints than having positive-semidefinite matrices.
I guess that https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-14.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-15.png should not be connected. Those quantities are related to the co-skewness* and they should be different when dependence is asymmetric.
Given a 3 parameter copula for pair https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-16.png, respectively describing global dependence (that can be related to L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-13.png” ne peut être affichée car elle contient des erreurs. and tail asymmetry (that should be related either to https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-14.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-15.pngdepending on the corner) it should be possible to derive any kind of triplets.

  • Running simulations

I have tried to generate a lot of samples. The idea is to draw pairs from mixtures of Marshall-Olkins’ copula (wich are asymmetric), where parameters where generated randomly. I draw 250 values according to that copula. Then, the third variable is obtained taking the product of normalized pairs (or Studentized). And then calculate the associated correlation matrix.

  1.  draw https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-17.png from https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-18.png (but actually any positive distribution should work) and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-19.png
  2. do 250 times the following (index https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-20.png)
  • draw https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-21.png with a Bernoulli https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-22.png distribution
  • draw&nbs
    p;https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-23.png from copula https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-24.png if https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-25.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-27.png if https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-26.png
  • set https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-28.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-29.png, and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-30.png
  1. estimate the correlation of sample https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-31.png

In dimension 3, I have the following graph,

In order to get a better understanding, I look at slices, for instance when https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-32.png. On the graph on the right, to dark area are values that can not be obtained since correlation have to be positive-semidefinite. Dots are points that have been obtained in a particular scenario. I have also included the convex hull if the points obtained with all the scenarios (here 50,000 scenarios), since I assume that the set of positive correlation is necessarily convex. Graphs below have been obtained when https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-33.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-34.png

So far I have not been able to obtain all possible correlations, I guess…  If anyone has suggestions…
* Recall that the co-skewness matrix, in dimension 2, for a pair https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-35.png is simply

https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-36.png

where https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-37.png.

Les tables de mortalité

Chose promise, chose dûe, un court billet expliquant les principales tables utilisées,

  • Les tables TD et TV 88-90

Ces tables datent un peu, et si je continue à en parler en cours, c’est parce qu’elles sont simples à utiliser (et pour continuer à me faire croire que je n’ai pas vieilli depuis mes études). Cette table est d’ailleurs tellement sérieuse qu’on la retrouve dans la loi (ici), dans un arrêté d’avril 1993. La table dite TD 88-90 (pour Décès) a été établie par l’INSEE suite aux observations réalisées entre 1988 et 1990 sur une population d’hommes. Elle était appliquée pour le calcul des primes des contrats d’assurance décès. La table dite TV 88-90 (pour Vie) a été établie par l’INSEE suite aux observations réalisées entre 1988 et 1990 sur une population de femmes. Elle était appliquée pour le calcul des primes des contrats d’assurance en cas de vie. Ces tables peuvent se récupérer à l’aide des codes suivants,

> TD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
+ sep=";",header=TRUE)
> TV=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
+ sep=";",header=TRUE)

Ces tables ont été remplacées par les tables dites TH et TF, respectivement.

  • Les tables TH et TF 00-02

Ces tables ont été établies à partir des données INSEE de la population française entre 2000 et 2002 et ont été lissées. Ce sont des tables générationnelles, qui nécessitent un correctif d’âge pour tenir compte des écarts de mortalité entre générations. Elles sont applicables à partir du 1er janvier 2006. L’institut des actuaires a proposé une “notice d’utilisation” en ligne ici, et Cimon en avait parlé sur son blog ().

  • Les tables TPRV 95

La table TPRV 93 (pour Table Prospective de Rente Viagère) est un extrait de la table dite plancher pour la tarification des contrats de rente viagère. Elle a été publiée par l’arrêté du 28 juillet 1993 (ici sans les annexes), et correspond à une table prospective qui retrace la mortalité des générations 1887 à 1993 (les tables prospectives sont au programme de Master 2).

La TPRV 93 représente la table complète de la génération 1950. La table est en ligne ici (en csv) lisible sous R avec le code suivant,

> TD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TPRV.csv",
+ sep=";",header=TRUE)

Mais ce n’est que la première étape. On utilise ensuite un décallage d’âge.

On a parlé de ce point en TD (mais sur les taux de mortalité), c’est aussi ce qui s’appelle l’hypothèse de Rueff, et traduit la diminution des taux de mortalité par un rajeunissement (i.e. un gain sur les âges), i.e.

Notons qu’en fait,  dépend du niveau des taux d’intérêt. L’ensemble des tables de décallage sont en ligne ici.

  • Les tables TGH et TGF 05

Il s’agit ici, là encore, d’une population de rentiers, comme pour la (ici pour les 38 pages du code). Aussi, ces tables ont été construites sur une population différente des tables TH et TF (construites sur l’ensemble de la population française). La méthodologie est décrite ici.

from two to three…

A short post to give more details about the final remark in the course of Financial Econometrics, and more precisely the formula that can be found in the book of Philip Jorion,

Note that this formula can be found (perhaps written with slight changes) in several papers, e.g. the following sentence (on the http://www.bis.org/website),

or the following formula, on documents from the Bank of England website,

I recently pulished (in French, here) a paper on the Value-at-Risk, including the following graph,

Usually, three times the average over 60 trading days is the larger component, but during the financial crisis, it turned out that the daily component was almost three times higher than the average value over the past the months (this fact was mention by Paul Embrechts in some conference in Paris on risk measures).
The interpreation of the multiplicative k coefficient (which is from 2 to 3 in some publications, or which exceeds 3 in others) has been proposed in a paper of Gerhard Stahl, entitled three cheers. The idea is to use the Bienaymé-Tchebychev inequality. For random variables with finite variance, then

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-01.png

Recall that this inequality is simply a corrolary of Markov’s inequality

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-02.png

or for any increasing function https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-99.png

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-03.png

(taking function https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-04.png, applied to https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-05.png). This upper bound can be far away from the true probability, see e.g. the gaussian case below, i.e. if  https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-06.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-07.png

 

> z = seq(0,3,by=.01)
> P = 2*dnorm(k)
> U = 1/z^2
> plot(z,P,type="l",lwd=2,col="red",xlab="",ylab="")

The ratio between the two is given below,

> plot(z,U/P,type="l",lwd=2,col="purple",xlab="",ylab="",ylim=c(0,10))

Note that it is possible to interprete the axis values as probabilities values, taking quantiles of the gaussian distribution

> plot(pnorm(z),U/P,type="l",lwd=2,col="purple",xlab="",
+ ylab="",ylim=c(0,10),xlim=c(.9,1))
> abline(h=3,lty=2)

The interpretation is that the upper bound is 3 times higher than the true probability in the Gaussian case when z is the quantile of the https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-06.png distribution associated with probability level 99%.
Note that

  • if z is the 95% quantile of the mathcal{N}(0,1) distribution, the ratio is 2 (1.92)
  • if z is the 99% quantile of the mathcal{N}(0,1) distribution, the ratio is 3 (3.04)
  • if z is the 99.55% quantile of the mathcal{N}(0,1) distribution, the ratio is almost 4 (3.88)
  • if z is the 99.75% quantile of the mathcal{N}(0,1) distribution, the ratio is 5 (5.04)

A more formal explaination is to assume that X is symmetric, and then

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-09.png

Thus, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-10.png, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-11.png, we have an upper bound for the  https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-12.pngValue-at-Risk,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-20.png

where the upper bound is the upper bound for the https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-12.png Value-at-Risk for any distribution with finite variance and centred.
If  https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-31.png, then https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-32.png, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-33.png.  But since, https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-33.png for a https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-36.png distribution, then

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-21.png

and further

https://perso.univ-rennes1.fr/arthur.charpentier/latex/BTC-22.png

Les modèles en réassurance

Publication d’un papier sur les modèles en réassurance dans la revue Risques (ici), suite à un questionnement sur la pertinence des modèles classiques utilisés par les réassureurs. La question initiale était partie de la constatation – que l’on retrouve ici ou – sur l’utilisation (ou la mauvaise utilisation) de modèles sophistiques en finance de marché, en essayant d’expliquer que – d’un point de vue épistémologique au moins – les modèles des réassureurs étaient plus robustes. En particulier, on notera que les modèles les plus anciens utilisés par les réassureurs (en particulier la loi de Pareto comme je l’avais évoqué ici) ont eu une légitimité pratique durant plusieurs décennies avant d’être justifiés par la théorie des valeurs extrêmes. Je ferais d’ailleurs bientôt un billet sur l’histoire des valeurs extrêmes en statistiques, en revenant en particulier sur les travaux de Gumbel ou de Fréchet.

Sinon le code utilisé dans le papier est en ligne ici, et la base

> sinpe = read.table("https://perso.univ-rennes1.fr/arthur.charpentier/sinpe.csv",header=TRUE,sep=";")
> head(sinpe)
      DSUR  MNTPE
1 19850206 240439
2 19851228 125674
3 19850504 488331
4 19851118 457347
5 19850220 990919
6 19851214 182939
> annee=as.numeric(substr(as.character(sinpe$DSUR),1,4))
> sinistres=sinpe$MNTPE[annee>1992]
> XS=sinistres/100000

On se limite ici aux sinistres survenus après 1992. Si l’on visualise ces montants de sinistres, on obtient

> datesur=as.Date(as.character(sinpe$DSUR),"%Y%m%d")
> jour=datesur[annee>1992]
> plot(jour,sinistres/100000,xlab="",ylab="Coût individuel",cex=.5,ylim=c(0,600))
> ded=50
> abline(h=ded)

On peut aussi faire le graphique Pareto-log-log,

> library(evir)
> n=length(X)
> plot(log(sort(X)),log((n:1)/(n+1)),
+ xlab="Coûts des sinistres (logarithme)",ylab="Fonction de survie (logarithme)",cex=.8)
> out <- gpd(XS, 15)
> XI=as.numeric(out$par.ests[1]); BETA=as.numeric(out$par.ests[2]) 
> x0=seq(2,8,.01)
> lines(x,-1/XI*(x-log(15)),col="red")

(l’ajustement de la loi de Pareto permettant de tracer la droite rouge) ou encore visualiser l’estimateur de Hill de l’indice de queue,

> hill(X)

Pour finir, le code suivant permet de calculer la prime pure (ou plus généralement une prime de Wang) soit de manière non-paramétrique (burning cost) ou bien en utilisant l’ajustement d’une loi de Pareto. Le papier étant un papier de vulgarisation, j’ai pris le seuil de manière arbitraire, sans aucune recherche de valeur “optimale”

> DEDUC = seq(10,50,by=5)
> lambda=0;  seuil=5
> WG1=WG2=rep(NA,length(DEDUC))
> for(k in 1:length(DEDUC)){
+ deductible=DEDUC[k]
+ out <- gpd(XS, seuil)
+ XI=as.numeric(out$par.ests[1]); BETA=as.numeric(out$par.ests[2]) 
+ G0=function(x){1-pgpd(x+seuil, xi = XI, mu = seuil, beta = BETA)}
+ G=function(x){(G0(x+deductible-seuil))/(G0(deductible-seuil))}
+ F=function(x){pnorm(qnorm(G(x))+lambda)}
+ (wang1=integrate(F, 0, Inf))
+ X=XS[XS>deductible]
+ n=length(X)
+ FS= function(z){
+ m=rep(NA,length(z))
+ for(i in 1:length(m)){
+ m[i]=sum(X>z[i]+deductible)/n}
+ return(m)
+ }
+ G=function(x){pnorm(qnorm(FS(x))+lambda)}
+ (wang2=sum(G(seq(0,800,.01))*.01))
+ WG1[k]=as.numeric(wang1$value)
+ WG2[k]=wang2
+ }
> plot(DEDUC,WG2-DEDUC,type='b',xlab="Niveau de la priorité ('00 000 euros)",ylab="Prime pure par sinistres réassuré",ylim=c(0,50))
> lines(DEDUC,WG1-DEDUC,type='b',col="red",pch=4)
> legend(10,50,c("Aujustement d'une loi de Pareto","'Burning cost'"),
+ col=c("red","black"),lty=1,cex=.8,pch=c(4,1))

Nikkei’s past experience vs. SP500 (in euros)

Following Michael’s idea (here), I wanted to go further, based on his intuition (and dataset that he kindly sent me, there). If we consider the two series of Nikkei index and SP500 index in euros, we have to following graph,

the code is simply the following (the merging function is simply here to avoid problem with different trading days: since we look at the index and not the return, it is the simplest way to deal with it).

> library(RODBC)
> base = odbcConnectExcel(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/spx_nky_eurusd.xls", 
+ readOnly = TRUE)
> series1 = sqlQuery(base,query="select * from [Tabelle1$A2:B8837]") # SPX
> series2 = sqlQuery(base,query="select * from [Tabelle1$D2:E8631]") # NKY
> series3 = sqlQuery(base,query="select * from [Tabelle1$G2:H8945]") # EURUSD
> odbcCloseAll()
> series4=merge(series1,series3)
> series4$SPEUR=series4$SPX/series4$EURUSD
> series5=merge(series4,series2)
> x=(as.Date(series5[,1])-as.Date("01/01/0000","%d/%m/%Y"))/365.25
> yl=range(series5[,4])
> xl=c(1975,2010)
> plot(x,series5[,4],axes=FALSE,xlab="",ylab="",type="l",
+ lwd=3,col="red",xlim=xl,ylim=yl)
> axis(1)
> axis(2, col="red")
> par(new=TRUE)
> yl=range(series5[,5])
> plot(x,series5[,5],axes=FALSE,xlab="",ylab="",type="l",
+ lwd=3,col="blue",xlim=xl,ylim=yl)
> axis(4, col="blue")
> mtext("SP500 in Euros", 2, line=2, col="red", cex=1.2)
> mtext("NKY", 4, line=2, col="blue", cex=1.2)

Those two series series seem to have a similar pattern, so an idea can be translate the SP500 on the left,

Interesting isn’t it ? Suppose that we want to forecast (or forsee ?) the SP500 in euro for the next 10 years…

People who enjoy charts would have here a nice tool…

Those two series are extremely correlated, with a correlation of 0.9572,

> X1=series5[2501:n,4]
> X2=series5[1:(n-2500),5]
> cor(X1,X2)
[1] 0.9572484

But are the two series cointegrated (see here, here or therefor material on cointegration) ? Well, using standard procedure, we first have to prove that the two series are integrated. First, let us look at the autocorrelograms,

At first sight, we confirm the economic intuition that those indices should be integrated. Standard tests confirm that intuition,

> acf(X2,lag=1000,col="light green")
> acf(X1,lag=1000,col="light green")
> library(tseries)
> adf.test(X1)
        Augmented Dickey-Fuller Test
data:  X1 
Dickey-Fuller = -1.0768, Lag order = 17, p-value = 0.9264
alternative hypothesis: stationary 
> adf.test(X2)
        Augmented Dickey-Fuller Test
data:  X2 
Dickey-Fuller = -1.2905, Lag order = 17, p-value = 0.8788
alternative hypothesis: stationary

But if we want to go further, we have to find the cointegration relationship between the two series. From an heuristic point of view, a linear regression should be a good proxy,

> reg=lm(X1~X2)
> plot(residuals(reg))

> acf(residuals(reg),lag=1000,col="light green")

> adf.test(residuals(reg))
        Augmented Dickey-Fuller Test
data:  residuals(reg) 
Dickey-Fuller = -5.176, Lag order = 17, p-value = 0.01
alternative hypothesis: stationary 
Message d'avis :
In adf.test(residuals(reg)) : p-value smaller than printed p-value
> pp.test(residuals(reg))
        Phillips-Perron Unit Root Test
data:  residuals(reg) 
Dickey-Fuller Z(alpha) = -46.9775, Truncation lag parameter = 11,
p-value = 0.01
alternative hypothesis: stationary 
Message d'avis :
In pp.test(residuals(reg)) : p-value smaller than printed p-value

When we look at the autocorrelation function, it looks like we do have a stationary series.
This idea is – more or less – the idea of Engle-Granger two step procedure. But actually, we can not directly use Dickey-Fuller’s test to see if residuals are integrated. This was proved in Phillips and Ouliaris (1990), who also proposed a test (see e.g. here),

> library(tseries); po.test(cbind(X1,X2))
        Phillips-Ouliaris Cointegration Test
data:  cbind(X1, X2) 
Phillips-Ouliaris demeaned = -53.1766, Truncation lag parameter = 57,
p-value = 0.01
Message d'avis :
In po.test(cbind(X1, X2)) : p-value smaller than printed p-value
Another similar function can be found in R
> library(urca)
> summary(ca.po(cbind(X1,X2)))
######################################## 
# Phillips and Ouliaris Unit Root Test # 
######################################## 
Test of type Pu 
detrending of series none 
Call:
lm(formula = z[, 1] ~ z[, -1] - 1)
Value of test-statistic is: 45.2032 
Critical values of Pu are:
                  10pct    5pct    1pct
critical values 20.3933 25.9711 38.3413

Thus, we has to admit that those series are cointegrated.

Based on that idea, it is possible to model the stationary component, and forecast it for the next ten years, based on the assumption that we know the behavior of one time series. Hence, if we add the confidence interval due to the stationary component uncertainty, we have the following graph,

 Of course, again, only uncertainty related to the stationary process is considered here….

Modélisation des tables de mortalité

En cours d’actuariat, nous commencerons à parler des tables de mortalité prospectives.

Côté théorie, je peux renvoyer vers quelques slides de Pierre Devolder, ici. Je renvoie ici pour un court papier de vulgarisation sur le sujet. Sur l’ajustement de lois de mortalité prospectives, je renvoie au très beau mémoire en ligne ici (ou  pour des choses très proches en anglais). Pour l’approche statique, je renvoie également ici pour une méthode simple d’ajustement de loi de Makeham. Pour les aspects dynamique, il y a des documents techniques ici ou .

Sinon pour les données que nous utiliserons en TD, elles sont ici (pour l’exposition, i.e. le nombre de vivants) et  (pour les décès). Il faut retravailler un peu les données, à cause d’une valeur “110+”, mais rien de plus simple,

> DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> EXPOSURE=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Exposures_1x1.txt",header=TRUE)
> DEATH$Age=as.numeric(as.character(DEATH$Age))
Message d'avis :
NAs introduits lors de la conversion automatique 
> DEATH$Age[is.na(DEATH$Age)]=110
> EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age))
Message d'avis :
NAs introduits lors de la conversion automatique 
> EXPOSURE$Age[is.na(EXPOSURE$Age)]=110

Coté code, je renvoie à un précédant billet sur le sujet, ici, et  pour le package (non officiel) de démographie développé par Rob Hyndman.

Enfin, à la fin du cours, on parlera un peu de modèles de durée. J’essayerai de taper quelques choses dans les jours à venir, mais côté TD, on utilisera ce que fait Frédéric Planchet (qui connaît fort bien le sujet), et qui se trouve en ligne ici. Les données sont en ligne ici, avec un peu de code R pour retraiter les donnéesici et , et un peu de code pour des modèles classiques, comme du modèle de Cox (ici) ou du Kaplan Meier ().

Les codes R utilisés en TD sont les suivants (avec en plus les graphiques obtenus),

> D=DEATH[DEATH$Year==an,]
> E=EXPOSURE[EXPOSURE$Year==an,]
> MU = D[,3:5]/E[,3:5]
> plot(0:110,log(MU[,1]),type="l",col="red")
> lines(0:110,log(MU[,2]),col="blue")

Avec ci-dessus les log des taux de décès instantannés pour l’année 1986.

> PH=matrix(NA,111,111) # ligne X, colonne H
> PF=PH
> for(x in 0:110){
+ PH[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,2]))
+ PF[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,1]))
+ }
> x=0
> plot(1:111,PH[x+1,],ylim=c(0,1),type="l",col="blue")
> lines(1:111,PF[x+1,],col="red")

pour la fonction de survie à la naissance (ci-dessus) ou à 50 ans (ci-dessous),

> somme=function(x){sum(x,na.rm=TRUE)}
> EH=EF=rep(NA,111)
> EH=apply(PH,1,somme)
> EF=apply(PF,1,somme)
> plot(0:110,EH,type="l",col="blue")
> lines(0:110,EF,col="red")

Ce qui correspond aux espérances de vie résiduelles.

> for(k in 1:3){
+ Z=MU[,k]
+ Z[is.nan(Z)==TRUE]=NA
+ MU[,k]=Z
+ }
> AGE=unique(DEATH$Age)
> ANNEE=unique(DEATH$Year)
> MUF=matrix(MU[,1],length(AGE),length(ANNEE))
> MUH=matrix(MU[,2],length(AGE),length(ANNEE))
> persp(AGE[1:100],ANNEE,log(MUH[1:100,]),
+ theta=-30,col="light green",shade=TRUE)
> persp(AGE[1:100],ANNEE,log(MUH[1:100,107:1]),
+ theta=-30,col="light green",shade=TRUE)
> image(AGE,ANNEE,log(MUH))

On a retourné le graph en mettant les années à l’envers, mais on peut aussi regarder la surface pour les femmes, ainsi que l’allure des courbes de niveau

> persp(AGE[1:100],ANNEE,log(MUF[1:100,107:1]),
+ theta=-30,col="light green",shade=TRUE)
> image(AGE,ANNEE,log(MUF))

Mais nous n’avons ici qu’une lecture statique des tables, et nous ne suivons pas vraiment un individu: en même temps que la personne vieilli, le temps passe…. Considérons une personne née en 1900, et regardons (ex post) l’évolution de sa probabilité de décès (en regardant non pas la probabilité de décéder à 18 d’une personne de 18 ans en 1900, mais en tenant compte du fait que la personne a fêté ses 18 ans en 1918).

> MUFP=MUHP=rep(NA,111)    # MU PROSPECTIF
> for(i in 1:111){
+     MUHP[i]=MUH[i,colonne+i]
+     MUFP[i]=MUF[i,colonne+i]
+ }
> plot(0:110,log(MUHP),col="purple",type="l")
> D=DEATH[DEATH$Year==annee.naissance,]
> E=EXPOSURE[EXPOSURE$Year==annee.naissance,]
> MU = D[,3:5]/E[,3:5]
> lines(0:110,log(MU[,2]),col="blue")

On peut aussi l’ajustement par un modèle à la Lee-Carter (via la fonction fournie par LifeMetrics)

> source("https://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")
> DEATH=DEATH[DEATH$Age<90,] # virer age>90
> EXPOSURE=EXPOSURE[EXPOSURE$Age<90,] # virer age>90
> XV=unique(DEATH$Age)
> YV=unique(DEATH$Year)
> ETF=t(matrix(EXPOSURE[,3],length(XV),length(YV)))
> DTF=t(matrix(DEATH[,3],length(XV),length(YV)))
> ETH=t(matrix(EXPOSURE[,4],length(XV),length(YV)))
> DTH=t(matrix(DEATH[,4],length(XV),length(YV)))
> WA=matrix(1,length(YV),length(XV))
> LCF = fit701(xv=XV,yv=YV,etx=ETF,dtx=DTF,wa=WA)
-1166100 -> -263940.2  ->-204175.2 
-188375.9 -> -185822.4  ->-184403.5 
-183223.3 -> -183069.8  ->-182350.3 
-181701.1 -> -181600.6  ->-181145.5 
-180729.6 -> -180671.5  ->-180371.8 
-180101.9 -> -180068.6  ->-179871.2 
> LCH = fit701(xv=XV,yv=YV,etx=ETH,dtx=DTH,wa=WA)
-2564403 -> -1028776  ->-953629.9 
-881334.2 -> -868279.3  ->-834413.6 
-804946.8 -> -802817.2  ->-786343 
-772775 -> -772007.5  ->-764501.7 
-758653.8 -> -758399  ->-755212 
-752931 -> -752851  ->-751620.2 
-750811.1 -> -750786.7  ->-750351.4 
> plot(LCF$x,LCF$beta1,type="l",col="red")
> lines(LCH$x,LCH$beta1,col="blue")

> plot(LCF$x,LCF$beta2,type="l",col="red")
> lines(LCH$x,LCH$beta2,col="blue")

> plot(LCF$y,LCF$kappa2,type="l",col="red")
> lines(LCH$y,LCH$kappa2,col="blue")

> MUF      = DTF/ETF
> MUH      = DTH/ETH
> MUFMODEL = MUF; MUHMODEL = MUH
> for(i in 1:nrow(MUF)){
+ for(j in 1:ncol(MUF)){
+ MUFMODEL[i,j] = exp(LCF$beta1[j]
+ +LCF$beta2[j]*LCF$kappa2[i])
+ MUHMODEL[i,j] = exp(LCH$beta1[j]
+ +LCH$beta2[j]*LCH$kappa2[i])
+ }}

soit, pour avoir un dessin,

> persp(LCH$y,LCH$x,log(MUH),
+ theta=45,col="light green",shade=TRUE)
> persp(LCH$y,LCH$x,log(MUHMODEL),
+ theta=45,col="light blue",shade=TRUE)

On peut ainsi comparer les taux bruts (en vert) avec les taux estimés (en bleu)

mais aussi regarder  se trouvent les erreurs le plus grandes,

Dis moi comment tu t’appelles…

A chaque fois que je prends le métro à Paris, ou que je traverse la gare Montparnasse en descendant du train, je suis surpris par les affiches de publicités qui nous harcèlent. Mais ce qui me surprendra le plus, ce sont les compagnes où on nous présente un Stéphane, un Mathieu, voire un Jean-Pierre, qui doivent nous rendre les affiches plus “humaines” (enfin, je pense). Les gens qui font du marketing ne sont pas forcément stupides (même si personnellement, me montrer la gentille bouille de Stéphane ne va pas me faire aimer quelqu’un qui verbalise des usagers à deux-roues), et je me suis demandé si le choix du prénom (oui, je suis persuadé que le brave Stéphane de l’affiche ne s’appelle pas Stéphane… depuis quand le monde de la pub’ coïnciderait avec le monde réel ? même mes enfants savent que ce n’est pas le cas) était vraiment pertinent…. Bref, j’ai demandé au réseau Quetelet un accès à la base des prénoms…

Pour mon nouvel ami Stéphane, on note effectivement que les Stéphane ont, en moyenne, un peu plus de 36 ans (comme sur la photo), et si on regarde les quantiles, la moitié ont entre 32 ans (quantile à 25%) et 41 ans (quantile à 75%). Bref, on cible bien les trentennaires.

 

Maintenant, on  peut se demander un peu plus précisément qui est ciblé par cette campagne, on peut regarder d’où viennent les Stéphane,

On retrouve bien un flash pendant une petite dizaine d’années, plutôt dans la partie nord de la France… Bref, la campagne semble très ciblée….

Pour les Mathieu, on veut nous vendre des téléphones portables, je crois…Et effectivement, on est sur un âge moyen beaucoup plus faible, de l’ordre de 25 ans. En revanche, un peu plus étalé que tout à l’heure, avec la moitié ayant entre 15 et 28 ans.


Ici la distribution est plus étalée sur la France, dans le temps,

Bon, le dernier, c’est Jean-Pierre…. alors Jean-Pierre, on ne peut pas trop le voir. On voit qu’il captive des enfants1 en leur racontant son expérience, ce qui correspondrait à une personne un peu âgée, avec de l’expérience. Alors effectivement, si on regarde dans l’historique, les Jean-Pierre ont 55 ans en moyenne, avec 50% entre 50 et 62 ans.


Et là aussi, on retrouve une population proche des Mathieu quant à sa répartition géographique (des pays de cheminots ?).

1 au delà de la fascination exercée par notre ami Jean-Pierre sur ces élèves, pour ma part, ce sont les élèves qui me fascinent… une classe presque exclusivement constituée de gauchers ! peut être pour accentuer le fait qu’il ne faut pas croire que la pub’, c’est comme dans la vraie vie, je ne sais pas… Remarquez, je n’ai rien contre les gauchers sur les photos ! mais c’est un peu comme les auvergnats: un ça va, mais c’est quand il y en a trop que ça devient gênant….

Vendredi 13 et probabilités….

Je ne sais pas si c’est la proximité avec Halloween, mais je trouvais ce vendredi 13 qui s’annonce avait quelque chose de spécial.. Et tant pis si après mes études sur les jours de naissance les jours de pleine lune, je passe définitivement pour un mystique. Disons plutôt que j’aime bien interroger les croyances de grands mères…. de la même manière que j’aime questionner les pratiques des actuaires (souvent pleines de bon sens, mais parfois un peu bancales).
Le problème avec les vendredi 13 est que l’on ne sait jamais si c’est un jour qui porte chance, ou au contraire un jour qui doit porter malheur… Dans le doute, je n’ai pas pu m’empêcher de mettre quelques trèfles à quatre feuilles, des brin de muguet et un pompon de marin.. Comme disait un très grand scientifique, “je n’y crois pas, mais il parait que ça marche même pour ceux qui n’y croient pas“. Je ne vais donc pas me priver !

Bon notons tout d’abord qu’avoir un vendredi 13 n’a rien d’exceptionnel… il y en a déjà au moins un par an. Mais ce qui est troublant c’est que le 13 d’un mois a plus de chances de tomber un vendredi qu’un autre jour. Pour avoir un cycle de semaines entières, il faut compter 400 ans, soit 20871 semaines (ce que l’on appelle un cycle grégorien). Mais en 400 ans, on a 4800 mois, soit 4800 “13”… qui n’est pas divisible par 7, ce qui tend à laisser penser que certains jours de la semaine tomberont (un peu) plus souvent que les autres… Et Brown a montré que c’était le cas (dès 1933, dans le American Mathemical Monthly1). Pour les amateurs, Jean-Luc Nothias avait fait un article en mars dernier sur le sujet (ici).
De plus, certains auteurs ont tenté des études statistiques sur le sujet. Par exemple, en 1993, le British Medical Journal avait publié un article sur le sujet (“Is Friday the 13th Bad for Your Health?” , ici) “Friday 13th is unlucky for some. The risk of hospital admission as a result of a transport accident may be increased by as much as 52 percent. Staying at home is recommended.” Bon, on va essayer de jeter un oeil….

  • Visites chez le médecin (généraliste) les vendredi 13

Les données a notre disposition (gentillement données par la Caisse d’Assurance Maladie) correspondent au nombre de consultations d’un médecin généraliste, par jour… Bon, ces données présentent un soucis technique: les données arrivent avec un (léger) délai, et les données du vendredi arrivent regroupées avec celles du samedi et du dimanche le mardi suivant. On ne peut donc pas comparer un vendredi au jeudi précédant par exemple… Pour éviter ce genre de soucis, la stratégie la plus simple consiste à comparer le nombre de consultations chez le médecin….

avec en rouge les vendredi 13, et en jaune la moyenne des vendredi 6 et 20…. Vu de loin, on consulte un peu plus les vendredi 13 que les vendredi avant et après. Si on fait un test de Student (quelle que soit l’hypothèse alternative) on rejette l’hypothèse comme quoi on conculte plus son médecin un vendredi 13…

> t.test(V13,V)
        Welch Two Sample t-test
data:  V13 and V
t = 0.6126, df = 9.052, p-value = 0.5552
alternative hypothesis: true difference in means is not equal to 0

Mais peut être est-ce par qu’une consultation chez le médecin généraliste n’est pas associé à de la “malchance” ?

  • Quid des catastrophes naturelles ?

Cherchons autre chose que l’on pourrait associer à cette prétendue “malchance“. Si on regarde les communes classées catastrophes naturelles en France, entre août 1982 et mars 2008, il y a eu 131 285 classements de communes comme touchées par une catastrophe naturelle.(et paru au Journal Officiel, ces données sont publiques). Seuls 393 classements ont été  fait un vendredi 13, ce qui pourrait sembler très faible… En fait, cette série est très ératique.

Le graphique montre l’évolution du nombre de classements les mois où il y a eu un vendredi 13. Statistiquement, on pourrait noter qu’il y a beaucoup plus de communes touchées par une catastrophe naturelle un mercredi 11 par exemple… Bref, la base semble inutilisable pour cette étude… simplement car un cinquième des classements concernent les tempêtes de décembre 1999 (ou plutôt les inondations qui ont suivies, car les tempêtes sont exclues du régime cat nat, en métropole).

  • Quid des gros sinistres de perte d’exploitation et des incendies entreprises…

J’aurais préféré avoir d’autres données d’assureurs, mais comme toujours, lorsque j’ai demandé les données, on m’a expliqué que c’étaient des données sensibles, et que je ne pourrais rien avoir ! Fort heureusement, il me reste quelques jeux de données…. en particulier une base de gros sinistres de sinistres en perte d’exploitation en France.
Cette fois, ma base va de début 1985 à fin 2000, avec en tout 2387 sinistres (cette fois l’interpréation devrait être plus simple car les données sont identifiées comme des sinistres). Au final 17 sinistres majeurs2 se sont produits un vendredi 13. Cette fois, le chiffre est significativement grand, car par exemple, seulement 4 se sont produit un jeudi 12, et 9 un samedi 14…

Là encore, si on regarde ce qui s’est passé les mois où il y avait un vendredi 13, on note qu’il y a eu beaucoup de sinistres ces jours là… Envie de faire un test ? On a en tout 2387 sinistres à faire tomber sur 5824 jours, parmi lesquels, on compte 27 vendredi 13. Bref, en supposant les sinistres indépendants, chaque sinistre à une probabilité de 0,463% de tomber un vendredi 13. Le nombre de sinistres tombant un vendredi 13 est donc supposé suivre une loi binomiale

On pourrait espérer avoir 11 sinistres qui tombent un vendredi 13 (et on en observe 17). Si on calcule le quantile à 95%, on tombe sur 17. En fait, on avait environ 1 chance sur 20 d’avoir autant de sinistres ces jours là…. ce qui correspond exactement au seuil que l’on retient pour juger de la significativité d’un test (comme je l’évoquais un peu ici). Bref, on ne peut pas vraiment conclure….
En fait, comme il s’agissait de sinistres d’entreprises, je pensais qu’il y aurais un biais, car (par définition) les vendredi 13 tombent un vendredi, et que les jours de semaines devaient avoir plus de sinistres que les jours de week end. Mais ce n’est pas le cas

Bref, on est dans la zone d’incertitude des tests statistiques….
Mais on peut aussi utiliser une base d’incendies intreprises que j’ai sous le coude depuis quelques années, avec des sinistres survenus entre 1982 et 1996. Cette base contient 9613 gros sinistres: En utilisant la même méthode que pour les sinistres de perte d’exploitation, on note que beaucoup d’incendies sont survenus un vendredi 13 dans cette base: 52. C’est un peu plus que ce que nous aurions pu espérer, à savoir 45,62. Là encore, en notant que le nombre d’incendies devrait suivre une loi binomiale, notons qu’il y avait 15,4% de chances d’atteindre 52 sinistres, ce qui n’est – cette fois-ci – pas particulièrement rare.

Par contre, sur cette base, l’effet week end semble plus prononcé, ce qui semble aller dans le sens que ces 52 sinistres, c’est loin d’être exceptionnellement elévé pour un jour de semaine.

On pourrait aussi parler de karma, ou de chose dans le genre….

1 B. H. Brown, “Solution to Problem E36”, American Mathematical Monthly, vol. 40, issue 10 (1933), p. 607 Pour faire simple, un cycle grégorien correspond à 400 années, soit 136097 jours, ou encore 20 871 semaines. Si on regarde la distribution des jours on voit que samedi et jeudi tombent 684 fois, lund et mardi 685, mercredi et dimanche 687, et vendredi 688. Bon, d’accord, je pinaille…

2 pour les plus curieux, ou les plus sceptiques, je mets ici les codes pour ceux qui veulent jouer avec les données… c’est un exercice amusant pour jouer avec la manipulation de dates sous R…

> B1=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/data-cat-nat.txt",header=TRUE)
> B2=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/data-pe.txt",header=TRUE)
> B4=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/data-incendies.txt",header=TRUE)
> B3=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/date-v13.txt",header=TRUE) 
> B2$DATEPE=as.character(B2$DATEPE)
> x2=as.Date(B2$DATEPE,"%d/%m/%Y")   
> B3$date=as.character(B2$date)
> x3=as.Date(B3$date,"%d/%m/%Y")
> p2=sum((x3>=min(x2))&(x3<=max(x2)))/as.numeric(max(x2)-min(x2))
> n2=length(x2);  n2*p2
> qbinom(.05,n2,p2); qbinom(.95,n2,p2)
> pbinom(sum(x2%in%x3)  ,n2,p2)
> NOM=weekdays(x3[1]+1:7-5)
> S=rep(NA,7)
> for(i in 1:7){S[i]=sum(weekdays(x2)==NOM[i])}
> barplot(S,names.arg=NOM,col="light blue")


voilà pour les sinistres en perte d’exploitation… le reste n’est qu’une variation sur le même thème…

Séminaire Probabilité et Statistique, UBO, Brest

Talk at the statistical seminar at the Université de Bretagne Occidentale, in Brest, Wednesday May 6th Tuesday May 5th, 14h (in  10 days), on “multivariate extremes. Slides can be found here.

The talk will give a detailed introduction on multivariate extremes and related concepts. Then the case of Archimedean copula will be fully described (following the paper with Johan Segers).

[04/05/2009]: some applications in risk management will be shown at the end of talk, as well as some news things on spatial correlation.

and in order to illustrate tail convergence of Archimedean copulas, I have uploaded two animations, with tail independence below,

with tail dependence (or asymptotic dependence),

Short course on catastrophes: the Mexican earthquake

One case study of catastrophes covers will be the Mexican cat bond. The dataset of past historical events in Mexico can be downloaded here.

Since everything we need is in this dataset, everything can be done using R.

mex=read.table("https://blogperso.univ-rennes1.fr/arthur.charpentier/public/conferences/EQ-mex-6-5.txt") 
mex$D=mex$V1*10000+mex$V2*100+mex$V3 
mex$Date=as.Date(as.character(mex$D),"%Y%m%d")

 

Convolution(s) et lois composées, aspects numériques

Nous avons vu en cours que l’outil de base pour travailler sur les sommes de variables aléatoires était la convolution,

(f * g )(t)  , stackrel{mathrm{def}}{=} int_{-infty}^{infty} f(tau)cdot g(t - tau), dtau
= int_{-infty}^{infty} f(t-tau)cdot g(tau), dtau.

permettant d’obtenir la densité ou la fonction de répartition de la somme de deux variables indépendantes.

Cet outils permet en théorie d’obtenir la forme d’une loi composée, i.e. d’une somme dont le nombre de termes est aléatoire.

Pour les calculs numériques de la convolée, Vincent Goulet a programmé pas  mal de choses dans la librairie actuar. Numériquement, il faut discrétiser la loi des coûts de sinistres. Si on suppose par exemple que les coûts individuels sont discrets, on peut obtenir la distribution simplement

> fx <- c(0, 0.15, 0.2, 0.25, 0.125, 0.075,
+         0.05, 0.05, 0.05, 0.025, 0.025)
> pn <- c(0.05, 0.1, 0.15, 0.2, 0.25, 0.15, 0.06, 0.03, 0.01)
> Fs <- aggregateDist("convolution", model.freq = pn,
+                     model.sev = fx, x.scale = 25)
> summary(Fs)
Aggregate Claim Amount Empirical CDF:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0   150.0   275.0   314.5   400.0  2000.0 
> c(Fs(0), diff(Fs(25 * 0:21))) # probability mass function
 [1] 0.05000000 0.01500000 0.02337500 0.03467500 0.03257656 0.03578639
 [7] 0.03980787 0.04356232 0.04751800 0.04903380 0.05189806 0.05137886
[13] 0.05118691 0.05030486 0.04818189 0.04575882 0.04280890 0.03937836
[19] 0.03574568 0.03196808 0.02832446 0.02478833
> plot(Fs)

On peut aussi utiliser des méthodes d’approximation récursives,

> Fs <- aggregateDist(“recursive”, model.freq = “poisson”,
+                     model.sev = fx, lambda = 3, x.scale = 25)
> plot(Fs)

Ce qui donne le graphique suivant

Des méthodes d’approximation de type normal power sont aussi programmées, ainsi que les techniques de simulations. Par exemple ici pour une survenance Poissonienne avec des coûts Gamma, on écrit

> Fs <- aggregateDist("recursive", model.freq = "poisson",
+ model.sev = fx, lambda = 3, x.scale = 25)
> plot(Fs)

Enfin, pour lier avec le DM que je donnerais mercredi, Vincent Goulet a aussi programmé une méthode d’estimation de quantile,

> model.freq <- expression(data = rpois(3))
> model.sev <- expression(data = rgamma(100, 2))
> Fs <- aggregateDist("simulation", nb.simul = 1000,
+                     model.freq, model.sev)

Estimation de quantile par noyau beta

Le papier sur l’estimation de quantile par noyau beta, coécrit avec Abder Oulidi, est accepté pour publication dans Statistics and Computing, http://link.springer.com/…

In this paper we propose several nonparametric estimators of quantiles based on Beta kernel and applied to transformed data by the generalized Champernowne distribution initially fitted to the data. A Monte-Carlo based study will show that those estimators improve the efficiency of a traditional ones, not only for light tailed distributions, but also heavy tails, when the probability level is close to 1. We also compare these estimators with the Extreme Value Theory Quantile applying to Danish data on large fire insurance losses.