actuariat - M2-09/10 - Arthur Charpentier

actuariat - M2-09/10

lundi 8 mars 2010

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

Par Arthur Charpentier le lundi 8 mars 2010, 10:37

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 là)

Par exemple pour l’indice de Gini à un noeud, on calcule
$http://perso.univ-rennes1.fr/arthur.charpentier/latex/tree-01.png$
où $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

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é....

un commentaire aucun rétrolien

samedi 16 janvier 2010

Vraisemblance, quasi-vraisemblance, ou pseudo-vraisemblance ?

Par Arthur Charpentier le samedi 16 janvier 2010, 01:56

Un petit point technique sur les GLM (que j’ai passé rapidement en cours). Pour revenir un peu en arrière, il existe plusieurs méthodes usuelles d’estimation de paramètre. Supposons que l’on dispose d’un échantillon tiré suivant une loi $http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais01.png$
La méthode du maximum de vraisemblance consiste à maximiser la vraisemblance c’est à dire la loi jointe du n-uplet. C’est une vraie loi. Dans le cas discret

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais02.png$

soit, en supposant les observations indépendantes et de même loi

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais03.png$

Dans le cas continu,

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais04.png$

On suppose que l’on connaît ici parfaitement la famille de loi, et que connaître la valeur $http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais0.png$ du paramètre suffit à caractériser complètement la loi. Le modèle est alors parfaitement spécifié. Parfois, on peut avoir du mal à spécifier complètement la loi (et en fait on n’en a pas vraiment besoin... On peut supposer que connaître le paramètre nous permet de caractériser les premiers moments (et pas complètement la loi). Ca sera l’idée de la pseudo-vraisemblance. Une pseudo vraisemblance adaptée à l’ordre 1 sera une fonction telle que l’espérance coïncide, i.e.

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais05.png$

Autrement dit la pseudo-vraisemblance n’est pas forcément la vraie loi des observations, mais elle nous garantie qu’ici, le premier moment est le bon. Sur la méthode du pseudo du pseudo-maximum de vraisemblance, je renvoie vers une très bon survey fait par Alain Trognon (ici) voilà une vingtaine d’année.
Une autre méthode est basée sur la quasi-vraisemblance. Ici, on ne fait pas toujours pas vraiment d’hypothèse de d’hypothèse de loi. Pour reprendre les notations des modèles GLM, on suppose juste que

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais06.png$

On suppose de plus qu’il existe un lien entre la variance et l’espérance, i.e.

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais07.png$

On introduit alors la loi de

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais08.png$

Alors

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais09.png$

optime une "quasi-vraisemblance", qui permet généralement d’obtenir un bon estimateur du paramètre (il est asymptotiquement sans biais et et optimal). Son avantage est de ne pas reposer sur la spécification d’une loi. Notons que l’on peut obtenir la variance de l’estimateur

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/vrais10.png$

Pour aller un peu plus loin sur ces histoires de quasi-vraisemblance, on notera qu’il existe, sous R, des lois quasi (quasi poisson ou quasi binomiale) parmi les familles de lois exponentielles. Les lois de Poisson et binomiales sont des lois à un paramètre, c’est à dire que l’on force \phi à être égal à 1. Ce paramètre est simplement un modèle de surdispersion (dont on parlait tout à l’heure, ici). Mais si \phi n’est pas égal à 1, la loi exponentielle n’est plus nécessairement une loi de probabilité.
> X=runif(100)
> Y=rpois(100,lambda=exp(1+X))
> summary(glm(Y~X,family=poisson))
Call:
glm(formula = Y ~ X, family = poisson)
Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.9263     0.1125   8.233 < 2e-16 ***
X             1.0987     0.1698   6.469 9.86e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 138.25 on 99 degrees of freedom
Residual deviance: 94.62 on 98 degrees of freedom
AIC: 421.21
Number of Fisher Scoring iterations: 5

> summary(glm(Y~X,family=quasipoisson))
Call:
glm(formula = Y ~ X, family = quasipoisson)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.9263     0.1078   8.593 1.36e-13 ***
X             1.0987     0.1627   6.752 1.04e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 0.917879)

    Null deviance: 138.25 on 99 degrees of freedom
Residual deviance: 94.62 on 98 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5
On note que dans ce dernier cas, le paramètre indique une sous-dispersion, avec une valeur de 0,9178.
Enfin, dans la terminologie, il existe aussi une vraisemblance empirique... histoire d’embrouiller encore davantage. Mais je n’en parlerais pas....
En revanche je reviens très bientôt pour finir l’économétrie de l’assurance auto...

aucun commentaire aucun rétrolien

Modéliser les nombres de sinistres

Par Arthur Charpentier le samedi 16 janvier 2010, 01:46

La surdispersion est souvent évoquée pour les lois binomiales ou de Poisson. J’avais évoqué le fait qu’empiriquement, la variance d’une variable de nombre est souvent plus élevée que l’espérance (alors que pour une loi de Poisson, les deux sont supposés égaux, j’en parlais en particulier ici). La sur-dispersion signifie que la variance de Y dépasse la variance attendue. Considérons un modèle binomial par exemple, i.e. Y prend deux valeurs, 0 ou 1.

On sait que

et que

Mais il est possible qu’en pratique, la variance soit différente (souvent supérieure) de cette valeur. Assez souvent, il existe des classes (clusters) non identifiées par les variables à notre disposition,

. Comme nous l’avions mentionné ici (par exemple), l’hétérogénéité non-observée entraine de la surdispersion.
Il existe des tests de surdispersion, détaillé dans le cas Poissonnien dans la section 3.4 du livre de Colin Cameron et Pravin Trivedi, Regression Analysis of Count Data. L’idée la plus simple est de noter que la loi de Poisson est un cas particulier de la loi binomiale négative, dans le cas où $http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset0.png$ . Il suffit de faire un simple test (ratio de vraisemblance par exemple) de nullité d’un coefficient. Le hic est qu’on teste si un coefficient est sur le bord des valeurs possibles, mais Moran (1971) par exemple propose des solutions. Sur l’estimation du paramètre de surdispersion, j’en ai parlé ici (en présentant la loi quasiPoisson).
Sinon parmi les notions classiques sur la modélisation des nombres de sinistres, il y a la variable offset. Pour expliquer cette notion, on suppose que l’on rajoute un terme au prédicteur linéaire, i.e.

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset01.png$

de telle sorte que

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset02.png$

L’exemple classique est obtenu dans le cas où le nombre suit une loi de Poisson et que l’on prend la fonction de lien canonique, i.e. la fonction log,

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset03.png$

Aussi, si l’on considère $http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset04.png$ on obtient un effet multiplicatif. Et $http://perso.univ-rennes1.fr/arthur.charpentier/latex/offset05.png$ est alors l’exposition, i.e. la durée de la période de couverture (une personne restant 6 mois au lieu d’un an devrait avoir deux fois moins de sinistres, en moyenne, qu’une personne ayant des caractéristiques identiques. Au niveau des unités, notre référence sera une exposition de 1. Il convient alors de normaliser l’exposition afin qu’une année corresponde à 1 (pour des contrats d’une période d’un an). On parlera parfois du "nombre d’années police" Sous R, on peut utiliser le code suivant
> glm(Y~X1+X2,family=poisson(link=log),offset=log(E))
Dans le cas d’une variable binomiale (y-a-t-il eu des sinistres, ou pas), par exemple avec un modèle logistique, on peut introduire la variable offset via la commande suivante sous R,
> glm(Y~X1+X2+offset(E),family=binomial)

aucun commentaire aucun rétrolien

vendredi 8 janvier 2010

Market valuation des garanties des produits d'assurance

Par Arthur Charpentier le vendredi 8 janvier 2010, 11:22

Pour le dernier cours d’actuariat 2, on parlera de market valuation, i.e. de valorisation financière. Je vais aussi en profiter pour éclairer Bruno (pour qui j’avais déjà commencé à poster des billets sur le contrôle optimal) sur "comment introduit-t-on la mesure risque neutre dans un passif ?" , sous-entendu en assurance-vie (en assurance non-vie c’est à mon avis plus compliqué, mais c’est la question que l’on se pose quand on valoriser des cat-bonds par exemple). Notons que ce problème s’est posé assez tôt dans la littérature économique, je renverrais par exemple au papier de Brennan et Schwartz en 1976.

Valorisation de produits financiers (en marché complet)

Je passerais un peu rapidement (en renvoyant ici ou là pour les grands principes), mais les techniques principales pour valoriser les actifs financiers sont les prix d’Arrow-Debreu, les probabilités risque-neutre, et les déflateurs. Mais auparavant, rappelons que le marché est composé de n actifs financiers (bien distincts) dont les prix et les payoffs sont parfaitement connus, et de n états de la natures, de telle sorte que la matrice des payoffs soit de plein rang. Introduisons le concept de portefeuillle de réplication. Compte tenu du fait que la matrice des payoffs est de plein rang, on a l’existence et l’unicité d’un tel portefeuille. Considérons un ensemble de payoff contingents aux états de la nature

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-17.png$

Il est alors possible de créer un portefeuille (dit de réplication) qui permette de toucher précisément ces payoffs à la date 1. Le prix (à la date 0) d’un tel actif est alors le prix du portefeuille de réplication.
Considérons les actifs rapportant 1 dans un état de la nature, et 0 sinon, à la date 1, i.e.

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-16.png$

appelés actifs d’Arrow-Debreu. On note $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-15.png$ le prix aujourd’hui d’un tel actif. Notons que le prix d’un actif sans risque, rapportant 1 demain dans tous les états du monde est $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-1.png$ . Aussi, si on note $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-14.png$ le taux sans risque, alors

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/correc.png$

Le prix de l’actif rapportant $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-12.png$ à la date 1 est alors $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-10.png$ .
Posons $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-09.png$ tel que $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-08.png$ soit une distribution de probabilité. Alors

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-09.png$

où l’espérance est prix sous la mesure de probabilité $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-pst.png$ appelée probabilité risque-neutre.
Je renvoie ici (par exemple) pour l’introdution de cette probabilité via les arbres binomiaux. Bref, c’est ainsi que l’on arrive au théorème fondamental de la valorisation d’actif en marché complet: le prix d’un actif contingent l’espérance sous probabilité risque neutre, des payoff actualisés.
Aussi, prenons l’exemple d’un put européen, de maturité T et de strike K. Le prix (à la date 0) du put est donné par

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-01.png$

Rappelons d’ailleurs que la formule de Black et Scholes donne un prix pour un tel produit financier (un put européen),

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/BSC01.png$

où

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/BSC02.png$

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/BSC03.png$

Voilà rapidement comment on valorise un produit financier en marché complet, en calculant la valeur d’un portefeuille de réplication.

Valorisation de produits d’assurance

Pour valoriser les produits d’assurance, classiquement, on calcule des primes pures, ou des valeurs actuelles probables. Par exemple en assurance-vie (car c’est essentiellement de ça dont on parle non ?), considérons simplement une garantie plancher, où l’on s’engage à verser 1€ au décès d’un souscripteur. La valeur actuelle probable (que l’on calcule sous la probabilité historique) est alors

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-02.png$

que l’on pourrra réécrire

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-03.png$

(où j’adopte une convention de temps continu pour les taux pour garder une analogie avec le cas financier, mais ça n’est pas vraiment important ici). Bref, par rapport à la valorisation financière on valorise toujours des flux futurs (payoffs) actualisés, mais en prenant l’espérance sous la probabilité historique. Et c’est généralement comme ça que l’on valorise un produit d’assurance: les risques ne sont pas échangés sur des marchés liquides, on n’a donc pas vraimentt de marché complet, et les assureurs espèrent simplement que la mutualisation (via la loi des grands nombres) leur permettra d’avoir égalité entre les flux entrants et les flux sortants.
Remarque: j’avais fait l’an dernier un rapide billet sur la transformée d’Esscher (ici) qui présentait une méthode conciliant les deux approches. Mais je suis toujours à la recherche du papier original (datant de 1932) pour proposer une nouvelle discussion sur ce sujet...

Valorisation des garanties planchers dans les contrats d’assurance vie, une approche heuristique

Là où les choses se compliquent c’est quand on mélange des garanties d’assurance (un décès) avec des garanties financières (sur des taux ou bien un actif financier pour des contats en unité de compte). Considérons pour cela un contrat décès, où l’assureur s’engage à payer, au moment du décès à la date $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-tau.png$ le maximum entre la somme investie en actif risqué à la signature du contrat, $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-07.png$ et le niveau atteint par l’actif lors du décès, $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-06.png$ . Autrement dit, l’assureur s’engage à verser $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-05.png$ . De manière plus général, disons qu’il existe une valeur plancher en dessous de laquelle on ne descendra pas, que l’on notera K.
Notons que l’on peut réécrire le payoff

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-04.png$

Le terme de gauche est la valeur de l’actif, et le second terme peut être interprété comme la valeur de la garantie plancher. Par la suite, on ne s’intéressera qu’à la valeur de cette garantie. Le flux versé (pour cette garantie) à la date t, pour un individu d’âge x à la souscription est alors tout simplement

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-04.png$

(les dates étant supposées discrètes ici). Si l’on cherche à calculer la valeur actuelle probable de l’ensemble des flux, on voit que deux espaces de probabilités intervienne: l’espace des probabilités actuarielles (sous lesquelles on couvre le risque de décès en supposant les risques mutualisables, i.e. l’univers historique $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-05.png$ ) et l’espace des probabilités financières (sous lesquelles on couvre le risque financiers, i.e. l’univers risque neutre $http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-06.png$ ). Autrement dit, la valeur actuelle probable s’écrit

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-07.png$

que l’on pourrrait réécrire

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-C-08.png$

où le terme à droite est le prix d’un call de strike K de maturité t.
Notons que la commission de contrôle s’était penché sur le sujet il y a quelques années, via le mémoire d’actuariat de Sylvain Merlus et Olivier Pequeux (ici), et c’est ce genre de formule qui avait été proposé, à savoir de faire somme somme de prix de put. Et c’est cette approche que la réglementation a imposé de valoriser cette garantie.
Pour aller plus loin, il avait aussi été précaunisé d’utiliser la formule de Black et Scholes pour valoriser ces puts, ce qui a donné la formule "standard"

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-03.png$

où

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-02.png$

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/market-valuation-01.png$

Voilà comment on valorise une telle garantie. Je n’ai parlé de la couverture, car on voit qu’elle ne serait pas parfaite, loin de là.
Comme toutes les formules standards, elle a l’avantage de la simplicité.
Bon formellement, si on veut faire ça proprement, ça risque d’êetre plus compliqué, et il faudrait parler des mesures risques neutre forward, mais on fera ça une autre fois. Sinon tout ça est repris dans une multitude de documents que l’on peut trouver ici, là, voire ici ou encore là. Et pour aller plus loin, il y a eu récemment pas mal de livres qui abordent ces sujets....

aucun commentaire aucun rétrolien

mardi 15 décembre 2009

Chain Ladder dans les triangles incomplets

Par Arthur Charpentier le mardi 15 décembre 2009, 15:37

Suite du billet (ici) sur le provisionnement en assurance dommage. On sait depuis le papier d’England et Verrall (1994) que la régression log-Poisson coïncide avec la méthode Chain-Ladder, en tous les cas en moyenne (on peut aller jeter un oeil ici pour les arguments théoriques). Pour les sceptiques, je renvoie au code suivant.
> library(ChainLadder)
> an <- 10; ligne = rep(1:an, each=an); colonne = rep(1:an, an)
> passe = (ligne + colonne - 1)<=an; n = sum(passe)
> PAID=GenIns; INC=PAID
> INC[,2:an]=PAID[,2:an]-PAID[,1:(an-1)]
> Y = as.vector(INC)
> lig = as.factor(ligne)
> col = as.factor(colonne)
> base = data.frame(Y,col,lig)
> reg=glm(Y~col+lig,data=base,family="poisson")
> sum(exp(predict(reg,newdata=base))[passe!=TRUE])
[1] 18680856
On peut vérifier avec la fonction MackChainLadder de library(ChainLadder), i.e.
> MackChainLadder(GenIns)
MackChainLadder(Triangle = GenIns)
      Latest Dev.To.Date Ultimate      IBNR Mack.S.E CV(IBNR)
1 3,901,463      1.0000 3,901,463         0         0      NaN
2 5,339,085      0.9826 5,433,719    94,634    71,835    0.759
3 4,909,315      0.9127 5,378,826   469,511   119,474    0.254
4 4,588,268      0.8661 5,297,906   709,638   131,573    0.185
5 3,873,311      0.7973 4,858,200   984,889   260,530    0.265
6 3,691,712      0.7223 5,111,171 1,419,459   410,407    0.289
7 3,483,130      0.6153 5,660,771 2,177,641   557,796    0.256
8 2,864,498      0.4222 6,784,799 3,920,301   874,882    0.223
9 1,363,294      0.2416 5,642,266 4,278,972   970,960    0.227
10   344,014      0.0692 4,969,825 4,625,811 1,362,981    0.295
                  Totals
Latest:    34,358,090.00
Ultimate: 53,038,945.61
IBNR:      18,680,855.61
Mack S.E.: 2,441,364.13
CV(IBNR):           0.13
Bref, les deux coïncident. Tout en étant fondamentalement très différents, comme l’avait noté Mack et Venter (2000). Mais que se passe-t-il si on ne dispose pas d’un triangle de données, mais d’un losange (par exemple en enlevant les 6 données du coin en haut à gauche).
> M=GenIns
> M[1,1:3]=NA
> M[2,1:2]=NA
> M[3,1]=NA
> an <- 10; ligne = rep(1:an, each=an); colonne = rep(1:an, an)
> passe = (ligne + colonne - 1)<=an; n = sum(passe)
> PAID=M; INC=PAID
> INC[,2:an]=PAID[,2:an]-PAID[,1:(an-1)]
> Y = as.vector(INC)
> lig = as.factor(ligne)
> col = as.factor(colonne)
> base = data.frame(Y,col,lig)
> reg=glm(Y~col+lig,data=base,family="poisson")
> sum(exp(predict(reg,newdata=base))[passe!=TRUE])
[1] 17810100
Le montant de provision est ici inférieur au montant obtenu par Chain Ladder (et qui est correct pour avoir refait les calculs),
> MackChainLadder(M)
MackChainLadder(Triangle = M)
      Latest Dev.To.Date Ultimate      IBNR Mack.S.E CV(IBNR)
1 3,901,463      1.0000 3,901,463         0         0      NaN
2 5,339,085      0.9826 5,433,719    94,634    69,796    0.738
3 4,909,315      0.9127 5,378,826   469,511   118,278    0.252
4 4,588,268      0.8661 5,297,906   709,638   130,514    0.184
5 3,873,311      0.7973 4,858,200   984,889   260,065    0.264
6 3,691,712      0.7223 5,111,171 1,419,459   410,087    0.289
7 3,483,130      0.6153 5,660,771 2,177,641   557,519    0.256
8 2,864,498      0.4155 6,893,493 4,028,995   861,372    0.214
9 1,363,294      0.2341 5,824,117 4,460,823   998,750    0.224
10   344,014      0.0684 5,028,331 4,684,317 1,492,354    0.319
                  Totals
Latest:    34,358,090.00
Ultimate: 53,387,997.33
IBNR:      19,029,907.33
Mack S.E.: 2,533,277.54
CV(IBNR):           0.13
Bref, si l’on ne dispose que de données incomplètes, il faut faire attention quand on fait une régression log-Poisson...

aucun commentaire aucun rétrolien

vendredi 11 décembre 2009

Démographie avec R, suite

Par Arthur Charpentier le vendredi 11 décembre 2009, 01:31

Dans le dernier TD d’actuariat, nous avons utilisé la library(demography) développée par Rob Hyndman (ici). On peut reprendre les bases utilisées dans le premier TD (ici).
> library(forecast)
> library(demography)
> DEATH=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> EXPOSURE=read.table("http://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
> YEAR=unique(DEATH$Year);nC=length(YEAR)
> AGE =unique(DEATH$Age);nL=length(AGE)
> MUF =matrix(DEATH$Female/EXPOSURE$Female,nL,nC)
> MUH =matrix(DEATH$Male/EXPOSURE$Male,nL,nC)
> POPF=matrix(EXPOSURE$Female,nL,nC)
> POPH=matrix(EXPOSURE$Male,nL,nC)
On a alors les données prêtes à être transformées dans des données de démographie,
> BASEH=demogdata(data=MUH, pop=POPH,
+        ages=AGE, years=YEAR, type="mortality",
+       label="France", name="Hommes", lambda=1)
> BASEF=demogdata(data=MUF, pop=POPF,
+        ages=AGE, years=YEAR, type="mortality",
+       label="France", name="Femmes", lambda=1)
On peut alors utiliser les fonctions de démographie, dont la fonction permettant d’estimer les paramètres du modèle de Lee & Carter.
> LCH=lca(BASEH)
> plot(LCH$age,LCH$ax,col="blue")
> plot(LCH$age,LCH$bx,col="blue")

> LCHf=forecast(LCH,100)
> plot(LCH$year,LCH$kt,col="red",
+ xlim=c(1900,2100),ylim=c(-300,150))
> lines(LCH$kt[nC]+LCHf$kt.f$mean)
> lines(LCH$kt[nC]+LCHf$kt.f$lower,col="red")
> lines(LCH$kt[nC]+LCHf$kt.f$upper,col="red")

Sinon on restreint à la période après guerre,
> LCH0=lca(BASEH,years=1948:2005)
> plot(LCH$age,LCH$ax,col="blue")
> lines(LCH0$age,LCH0$ax,col="blue")
> plot(LCH$age,LCH$bx,col="blue")
> lines(LCH0$age,LCH0$bx,col="blue")

> LCH0f=forecast(LCH0,100)
> plot(LCH$year,LCH$kt,col="orange",
+ xlim=c(1900,2100),ylim=c(-300,150))
> points(LCH0$year,LCH0$kt,col="red")
> lines(LCH0$kt[58]+LCH0f$kt.f$mean)
> lines(LCH0$kt[58]+LCH0f$kt.f$lower,col="red")
> lines(LCH0$kt[58]+LCH0f$kt.f$upper,col="red")

> LCHT=lifetable(LCHf)
> LCHTu=lifetable(LCHf,"upper")
> LCHTl=lifetable(LCHf,"lower")
> plot(0:100,LCHT$ex[,5],type="l",col="red")
> lines(0:100,LCHTu$ex[,5],lty=2)
> lines(0:100,LCHTl$ex[,5],lty=2)
> plot(2007:2105,
+ LCHT$ex[1,2:100]-LCHT$ex[1,1:99])

Bon, le dernier graphe me laisse un peu sceptique: je pensais que l’on pourrait noter l’accroissement d’un trimestre d’espérance de vie tous les ans...
On peut aussi utiliser la library(gnm), et lancer une régression à l’aide d’un outils plus général (qu’avait présenté Heather à la conférence UseR! cet été, ici).
> library(gnm)
> Y=DEATH$Male
> E=EXPOSURE$Male
> Age= DEATH$Age
> Year=DEATH$Year
> I=(DEATH$Age<100)
> base=data.frame(Y=Y[I],E=E[I],
+ Age=Age[I],Year=Year[I])
> REG=gnm(Y~factor(Age)+
+ Mult(Exp(factor(Age)),factor(Year)),
+ data=base,offset=log(E),
+ family=quasipoisson)
Initialising
Running start-up iterations..
Running main iterations........................
Done
> plot(1:99,REG$coefficients[2:100],col="blue")
> plot(1899:2005,REG$coefficients[201:307],col="red")
Le code peut être un peu long à faire tourner, mais ce code permet d’implémenter n’importe quel modèle de démographie,

Je ferais plus tard un billet sur les modèles de durée (avec de la censure).

aucun commentaire aucun rétrolien

lundi 7 décembre 2009

Modélisation des tables de mortalité

Par Arthur Charpentier le lundi 7 décembre 2009, 22:27

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 là 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 là.
Sinon pour les données que nous utiliserons en TD, elles sont ici (pour l’exposition, i.e. le nombre de vivants) et là (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("http://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> EXPOSURE=read.table("http://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 là 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ées ici et là, et un peu de code pour des modèles classiques, comme du modèle de Cox (ici) ou du Kaplan Meier (là).
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]
+ }
Erreur : indice hors limites
> 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("http://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])
+ }}

> 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 où se trouvent les erreurs le plus grandes,

aucun commentaire aucun rétrolien

dimanche 6 décembre 2009

Crédibilité ou régression ?

Par Arthur Charpentier le dimanche 6 décembre 2009, 00:15

On présente généralement la crédibilité dans sa version dynamique, inspirée par une lecture bayésienne de mise à jour d’information (comme je le faisais ici par exemple). Mais cette méthode est vraiment beaucoup plus générale. Historiquement, elle a été introduite pour tarifer des grandes groupes (d’où l’utilisation classique sur les flottes), mais elle peut aussi de voir comme outils de lissage à la place des méthodes glm (si on les applique sur des variables par classes).

La crédibilité: pondérer une classe par rapport à un groupe

Pour rappels, dans une approche telle qu’elle a été présenté par Hans Bühlmann, on s’intéresse à

avec pour le premier terme la moyenne d’un groupe,

où

est le facteur de crédibilité, avec

, et où respectivement,

où $http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/theta-latex.png$ est la variable caractérisant les groupes. Bref, pour faire des calculs, il suffit de savoir calculer des variances intra et des variances inter. De manière formelle, sous R, on pourra utiliser un code de la forme suivante
> XBAR = tapply(X, as.factor(THETA), mean)
> n = tapply(X, as.factor(THETA), length)
> MU = weighted.mean(XBAR,n)
> VXintra = tapply(X, as.factor(THETA), sd)^2
> VXinter = weighted.mean(XBAR^2,n)-MU^2
> Z = n/(n+(VXintra/VXinter)
> pred = Z*XBAR + (1-Z)*MU
Pour l’instant, seule l’utilisation des glm avait été évoquée en tarification, mais cette méthode devrait permettre de lisser un peu les calculs de primes, en particulier sur les petites classes.

Discriminer par une variable de zonage

Prenons comme variable de classe (c’est à dire $http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/theta-latex.png$ ) la variable de zonage.
> BASE=read.table("D:\\r-data\\base_finale.txt",header=TRUE,sep=";")
> I=is.na(BASE$NBSIN_08)
> I=is.na(BASE$NBSIN_DV_08)
> THETA=BASE$ZONIER[I==FALSE]
> X=BASE$NBSIN_08[I==FALSE]
> XBAR = tapply(X, as.factor(THETA), mean)
> n = tapply(X, as.factor(THETA), length)
> MU =     weighted.mean(XBAR,n)
> VXintra = tapply(X, as.factor(THETA), sd)^2
> VXinter = weighted.mean(XBAR^2,n)-MU^2
> Z = n/(n+(VXintra/VXinter)
> pred = Z*XBAR + (1-Z)*MU
Essayons de comparer les fréquences obtenues par crédibilité avec des modèles glm,
> Y=X; X1=THETA
> Base=data.frame(Y,X1)
> reglm=lm(Y~as.factor(X1)+0,Base)
> regglm=glm(Y~as.factor(X1)+0,data=Base,family=poisson)
> cbind(n,XBAR,coefficients(reglm),exp(coefficients(regglm)),pred)
      n XBAR              pred
1     85 0.023 0.023 0.023 0.027
2     67 0.029 0.029 0.029 0.030
3    388 0.028 0.028 0.028 0.028
7    268 0.011 0.011 0.011 0.015
8    243 0.020 0.020 0.020 0.023
9    451 0.019 0.019 0.019 0.022
11 4104 0.026 0.026 0.026 0.026
12   540 0.018 0.018 0.018 0.020
13   244 0.024 0.024 0.024 0.026
14 4236 0.025 0.025 0.025 0.026
16    64 0.015 0.015 0.015 0.024
22 8048 0.033 0.033 0.033 0.033
26 10227 0.032 0.032 0.032 0.032
27 26542 0.032 0.032 0.032 0.032
28 9982 0.025 0.025 0.025 0.025
On note que sur les gros effectifs, on obtient la même chose (ou presque) qu’un modèle glm, mais que sur les petites classes, on s’adapte en essayant de prendre une valeur plus proche que l’ensemble de la population. Bref, la crédibilité est une manière élégante pour lisser les valeurs prédites...

Discriminer par l’ancienneté du permis et le sexe

Mais plus généralement, on peut essayer de prendre davantage de classes. Regardons par exemple l’ancienneté du permis.
> Y=BASE$NBSIN_08[I==FALSE]
> X2=BASE$ANC_PERMIS_08[I==FALSE]
> X3 = paste(as.character(X2),as.character(BASE$SEXE_08[I==FALSE]),sep="-")
Avant d’aller plus loin, regardons (rapidement) comment la fréquence évolue avec l’ancienneté du permis.
> library(gam)
> reggam=gam(Y~s(X2),poisson)
> plot(reggam,se=TRUE,col="red")
Étrangement, sur cette base, les "jeunes" conducteurs (ou disons plutôt les novices) sont incroyablement bons, et plus le temps passe, plus ils ont de chances d’avoir un accident.... Mais cette irrégularité (explicable par des raisons que je ne peux dévoiler pour des raisons de confidentialité) n’enlève rien à la généralité du raisonnement,
> THETA=X3
> X=Y
> XBAR = tapply(X, as.factor(THETA), mean)
> n = tapply(X, as.factor(THETA), length)
> MU =     weighted.mean(XBAR,n)
> VXintra = tapply(X, as.factor(THETA), sd)^2
> VXinter = weighted.mean(XBAR^2,n)-MU^2
> Z = n/(n+(VXintra/VXinter)
> pred = Z*XBAR + (1-Z)*MU
Comparons comme toute à l’heure avec des modèles de régression,.
> reglm=lm(Y~as.factor(X3)+0,Base)
> regglm=glm(Y~as.factor(X3)+0,data=Base,family=poisson)
> cbind(n,coefficients(reglm),exp(coefficients(regglm)),pred)
        n              pred
10-F 500 0.018 0.018 0.018
10-M 565 0.008 0.008 0.009
11-F 492 0.018 0.018 0.018
11-M 576 0.012 0.012 0.012
12-F 483 0.026 0.026 0.027
12-M 580 0.022 0.022 0.022
13-F 477 0.025 0.025 0.025
13-M 552 0.019 0.019 0.020
[...]
33-F 603 0.028 0.028 0.028
33-M 981 0.032 0.032 0.032
34-F 662 0.049 0.049 0.048
34-M 1071 0.027 0.027 0.027
35-F 619 0.053 0.053 0.051
35-M 993 0.029 0.029 0.029
36-F 529 0.034 0.034 0.033
[...]
54-M 276 0.057 0.057 0.052
55-F   60 0.083 0.083 0.054
55-M 258 0.062 0.062 0.056
56-F   56 0.089 0.089 0.055
56-M 240 0.020 0.020 0.021
57-F   42 0.095 0.095 0.052
57-M 208 0.067 0.067 0.058
Bref, on retrouve que pour les grandes classes (à partir de 500 personnes) les ordres de grandeurs sont très proches. En dessous, l’estimateur par crédibilité corrige afin de prendre en compte l’effet taille (et permet de lisser les prédictions).

aucun commentaire aucun rétrolien

samedi 5 décembre 2009

De la surdispersion des nombres

Par Arthur Charpentier le samedi 5 décembre 2009, 23:08

J’avais déjà fait un billet (ici) pour rappeler que la loi de Poisson était très attirante, mais souvent peu adapté sur des données trop hétérogènes. En faisant des régressions, on essaye d’expliquer l’hétérogénéité avec des combinaisons linéaires de variables explicatives, autant que faire se peut. Et parfois, ce n’est pas suffisant, il reste de l’hétérogénéité résiduelle.

Effets fixes et effets aléatoires en économétrie

Quand on cherche é expliquer une variable (dite endogène, Y) par d’autres (dites exogènes), il y a celles dont on peut facilement disposer et qui sont observables (notées X), et celles qui sont cachées, et dont l’effet n’en est pas moins réel (notées Z).

Formellement, on peut définir une vraisemblance sur l’ensemble de variables exogènes

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-01.png

et celle définie seulement sur les variables observables

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-02.png

Notons que

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-03.png

En actuariat, on s’intéresse à des fréquences de sinistres, supposées suivre des lois de Poisson (du moins si l’on arrive é constituer des classes suffisamment homogènes). Suppsons que

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-04.png

alors

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-05.png

où

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-06.png

On a alors un modéle dit à effets fixes, au sens où

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-07.png

où

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-08.png

La vraisemblance se calcule alors par intégration sur les effets aléatoires. La variance de cet effet aléatoire va induire de l’hétérogénéité. Le cas le plus classique est obtenu lorsque l’on suppose que U suit une loi Gamma, auquel cas la loi de Y conditionnelle à X est une loi binomiale négative.

Si l’on suppose que

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-09.png

, alors

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-10.png

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-11.png

Alors, en posant

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-12.png

on en déduit que

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-13.png

et donc

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-14.png

En revanche, si on calcule la loi nonconditionnelle, alors

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-15.png

alors que

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-16.png

soit

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-17.png

Tester la surdispersion

Un test graphique peut être facilement construit, en représentant la variance empirique dans un sous groupe de mêmes caractéristiques par rapport à l’espérance, par sous-classe. Dans un modèle d’équidispersion, les données seraient alignés sur une droite de pente 1.
> BASE=read.table("D:\\r-data\\base.txt",header=TRUE,sep=";")
> I=is.na(BASE$NBSIN_DV_08)
> X1=BASE$ZONIER_08[I==FALSE]
> Y=BASE$NBSIN__08[I==FALSE]
> reglmqp=glm(Y~as.factor(X1)+0,family=quasipoisson)
> reglmp =glm(Y~as.factor(X1)+0,family=poisson)
> summary(reglmqp)
Call:
glm(formula = Y ~ as.factor(X1) + 0, family = quasipoisson)
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)
as.factor(X1)1   -3.74950    0.73230 -5.120 3.06e-07 ***
as.factor(X1)2   -3.51155    0.73230 -4.795 1.63e-06 ***
as.factor(X1)3   -3.56311    0.31225 -11.411 < 2e-16 ***
[...]
as.factor(X1)27 -3.43071    0.03534 -97.091 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 1.072523)

On note qu’en ajustant une loi quasiPoisson, R nous dit qu’il y a un paramètre de surdispersion de l’ordre de 1,07 (correspond à un effet de surdispersion). Il s’agit de la courbe bleue sur le graphique ci-dessus. La courbe verte correspond à la régression (par moindres carrés) des variances empiriques sur les moyennes empiriques (sans constante, avec une pente de l’ordre de 1,033).

> M = tapply(Y, as.factor(X1), mean)
> V = tapply(Y, as.factor(X1), sd)^2
> plot(M,V")
> abline(0,1,col="red")
> abline(0, 1.076690 ,col="blue")
> abline(lm(V~0+M),col="green")
Pour un test plus formel, il convient de spécifier davantage l’hypothèse alternative. On supposera par exemple que

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-18.png

qui correspond tout simplement à un modèle avec effet aléatoire, et

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-19.png

est simplement la variance de l’effet aléatoire.
On cherche alors à tester

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-20.png

contre

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-21.png

On peut alors considérer la statistique de test suivante

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-22.png

Cette statistique doit suivre, sous H₀ une loi normale centrée réduite. Parmil les autres statistiques classiques, on pourra aussi utiliser

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-23.png

Sous R, ces tests ont été implémentés dans la library(AER), via la fonction dispersiontest().
> library(AER)
> dispersiontest(reglmp)
        Overdispersion test
data: reglmp
z = 5.9049, p-value = 1.764e-09
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion
1.076690
> dispersiontest(reglmp,trafo = 2)
        Overdispersion test
data: reglmp
z = 6.0101, p-value = 9.27e-10
alternative hypothesis: true alpha is greater than 0
sample estimates:
   alpha
2.566681
Ces deux tests sont détaillés dans les deux livres Cameron et Trivedi (1998, 2005).

Un autre exemple pourrait être le "nombre d’heures, en moyenne, passées à regarder la télévision" (via une des enquêtes de l’insee, ici). On se doute que la "consommation" de télévision (pour poursuivre dans la syntaxe de Patrick Le Lay) est liée à l’âge, et a priori pas de manière linéaire. On peut donc tenter un modèle gam,

Quid de la sous dispersion ?

Il existe des cas où les nombres sont sous dispersés. Regardons par exemple le nombre de redoublements de classes par des enfants scolarisés.
On se doute que l’âge de l’enfant intervient. On peut alors régresser suivant l’âge, ou bien régresser sur une constante, et prendre l’âge (ou plutôt le logarithme de l’âge si l’on considère une fonction de lien logarithmique) en variable offset. A partir de la base fournie par l’insee (ici), on obtient

> library(foreign)
> parent1=read.dbf("D:\\r-data\\parent1.dbf")
> parent2=read.dbf("D:\\r-data\\parent2.dbf")
> parent3=read.dbf("D:\\r-data\\parent3.dbf")
> base=data.frame(parent1,parent2,parent3)
> base0=data.frame(Y=base$NBREDOUB,
+ E=base$AGEA)
> reg1=glm(Y~E,data=base0,family=quasipoisson)
> summary(reg1)
Call:
glm(formula = Y ~ E, family = quasipoisson, data = base0)
Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.474161   0.089512 -38.81   <2e-16 ***
E            0.158290   0.005087   31.12   <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.8086682)
ou dans le cas de la régression où l’âge est déclaré en variable offset
> reg2=glm(Y~1,offset=log(E),data=base0,family=quasipoisson)
> summary(reg2)
Call:
glm(formula = Y ~ 1, family = quasipoisson, data = base0, offset = log(E))
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.65357    0.02421 -150.9   <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.7639523)
Bref, les deux régression nous confortent dans l’idée d’une sous-dispersion (mais les tests ne sont pas implémentés pour ça).
Si on pense que l’effet de l’âge n’est pas spécialement linéaire, on peut aussi tenter un modèle gam.

> library(gam)
> reg3=gam(Y~s(E),data=base0,family=quasipoisson)
> summary(reg3)
Call: gam(formula = Y ~ s(E), family = quasipoisson, data = base0)
(Dispersion Parameter for quasipoisson family taken to be 0.8095)
Ce qui confirme cet effet de sousdispersion. Et là encore, un rapide test graphique confirme (plus ou moins) cette intuition. La courbe rouge correspond là encore à l’équidispersion, la courbe bleue est de pente 0,81 (obtenue avec deux régressions) et la courbe mauve de pente 0,76.
Un autre exemple pourrait être le nombre de mariages d’une personne. Là aussi, ce nombre doit dépendre de l’âge de la personne interrogée.

> p1=read.dbf("D:\\r-data\\hdv1.dbf")
> p2=read.dbf("D:\\r-data\\hdv2.dbf")
> p3=read.dbf("D:\\r-data\\hdv3.dbf")
> p4=read.dbf("D:\\r-data\\hdv4.dbf")
> base=data.frame(p1,p2,p3,p4)
> reg1=glm(BNBMAR~AGEE,
+ family=quasipoisson)
> summary(reg1)
Call:
glm(formula = BNBMAR ~ AGEE, family = quasipoisson, data = base)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.1685735 0.0230460 -50.71   <2e-16 ***
AGEE         0.0188916 0.0004085   46.25   <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.3370578)
> reg2=glm(BNBMAR~1,offset=log(AGEE),data=base,family=quasipoisson)
> summary(reg2)
Call:
glm(formula = BNBMAR ~ 1, family = quasipoisson, data = base,
    offset = log(AGEE))
Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.081337   0.006818 -598.6   <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.3178957)

Trop de valeurs nulles ?

Il est possible, dans certains car qu’un problème vienne d’un trop grand nombre de valeurs nulles pour accepter une hypothèse de loi de Poisson. La library(pscl) propose une fonction zeroinfl() assez simple d’utilisation. Au niveau de la synthaxe, après le | figurent les variables expliquant le surnombre des valeurs nulles. On peut tester par exemple ce problème sur les redoublements (ou l’on peut suspecter que le "problème" vient du fait que trop peu de personnes ne redoublent pas, ou de manière duale, qu’il existe une spirale du redoublement)
> base=data.frame(parent1,parent2,parent3)
> base0=data.frame(Y=base$NBREDOUB,E=base$AGEA)
> Y=base0$Y;X1=base0$E
> library(pscl)
> reg4 <- zeroinfl(Y ~ E | E,
+ data = base0, dist = "poisson")
> summary(reg4)
Call:
zeroinfl(formula = Y ~ E | E, data = base0, dist = "poisson")
Count model coefficients (poisson with log link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.18824    0.19465 -11.242   <2e-16 ***
E            0.09077    0.01040   8.727   <2e-16 ***
Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.02989    0.60910   9.900   <2e-16 ***
E           -0.65280    0.07936 -8.225   <2e-16 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Number of iterations in BFGS optimization: 25
Log-likelihood: -2466 on 4 Df
Rapidement, notons que tous les coefficients sont significatifs, ce qui nous donne envie d’accepter l’hypothèse que la sous dispersion serait expliquée par un trop fort nombre de valeurs nulles.

aucun commentaire aucun rétrolien

mardi 1 décembre 2009

Erreur(s) de prédiction et provisionnement

Par Arthur Charpentier le mardi 1 décembre 2009, 13:26

Dans les modèles de provisionnement, on souhaite avoir un ordre de grandeur de la qualité de notre ajustement, de notre prédiction. Formellement, on s’intéresse à

qui peut se décomposer sous la forme suivante

qui a - à peu près - égal à

Dans le cas de la régression linéaire, où nous avions déjà vu ce genre de choses, nous avions interprété les deux termes de la manière suivante

le premier terme correspond à l’erreur intrinsèque d’un modèle d’ajustement, qui n’est jamais parfait (c’est la variance du bruit)
le second terme correspond à l’erreur d’estimation des paramètres

On retrouve ainsi les graphiques classiques, en traduisant les variances sous forme d’intervalle de confiance.

Dans le cas d’un modèle GLM, par exemple un régression Poissonnienne comme on en utilise dans les modèles de provisionnement cf ici par exemple), on s’intéresse aux incréments dans les triangles de paiements, avec

La delta-method (que j’avais évoquée ici) nous garantie qu’asymptotiquement

or dans le cas d’une régression Poissonnienne avec un lien logarithmique

de telle sorte que s’il l’on considère une loi de Poisson surdispersée,

pour la partie inférieure du triangle (inconnue). On peut également montrer que

Le montant de provision que l’on cherche à estimer est alors la somme des prédictions de paiements à venir (c’est tout du moins ce qu’impose la réglèmentation)

En faisant un peu de calculs, on peut facilement montrer que

Bon, reste à implémenter tout ça sous R. Pour le début, on récupère notre triangle fétiche,
> source("http://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
> library(statmod)
> an <- 6; ligne = rep(1:an, each=an); colonne = rep(1:an, an)
> passe = (ligne + colonne - 1)<=an; np = sum(passe)
> futur = (ligne + colonne - 1)> an; nf = sum(passe)
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> Y = as.vector(INC)
> lig = as.factor(ligne)
> col = as.factor(colonne)
On fait ensuite une régression Poissonnienne (ainsi qu’une version où on remplace les valeurs manquantes par presque rien essentiellement pour récupérer des matrices de taille convenables)
> CL <- glm(Y~lig+col, family=quasipoisson)
> Y2=Y; Y2[is.na(Y)]=.001
> CL2 <- glm(Y2~lig+col, family=quasipoisson)
Maintenant on est prêt, on peut faire tourner les calculs
> p = 2*6-1;
> phi.P = sum(residuals(CL,"pearson")^2)/(np-p)
> Sig = vcov(CL)
> X = model.matrix(CL2)
> Cov.eta = X%*%Sig%*%t(X)
> mu.hat = exp(predict(CL,newdata=data.frame(lig,col)))*futur
> pe2 = phi.P * sum(mu.hat) + t(mu.hat) %*% Cov.eta %*% mu.hat
> cat("Total reserve =", sum(mu.hat), "prediction error =", sqrt(pe2),"\n")
Total reserve = 2426.985 prediction error = 131.7726
Si on revient sur le code pour faire du bootstrap sur le triangle des incréments modélisée par une régression Poissonnienne. Mais auparavant, construisons une fonction de type Chain Ladder pour calculer le montant de provisions à partir d’un triangle,
> CL=function(triangle){
+ n=nrow(triangle)
+ LAMBDA=rep(NA,n-1)
+ for(i in 1:(n-1)){
+ LAMBDA[i]=sum(triangle[1:(n-i),i+1])/
+ sum(triangle[1:(n-i),i]) }
+ DIAG=diag(triangle[,n:1])
+ TOTO=c(1,rev(LAMBDA))
+ return(sum(cumprod(TOTO)*DIAG-DIAG)) }
On peut alors construire des faux triangles d’incréments par bootstrap des résidus (normalisés), et ensuite, on utilise une fonction de type Chain Ladder pour calculer les réserves,
> base=data.frame(Y,lig,col)
> REG=glm(Y~lig+col,family=quasipoisson)
> PROVISION=rep(NA,50000)
> for(k in 1:50000){
+ simE=sample(E,size=36,replace=TRUE)
+ bruit=simE*sqrt(exp(YP))
+ INCsim=exp(YP)+bruit
+ INCM=matrix(INCsim,6,6)
+ CUMM=INCM
+ for(j in 2:6){CUMM[,j]=CUMM[,j-1]+INCM[,j]}
+ PROVISION[k]=CL(CUMM)}
> plot(density(PROVISION))
> quantile(PROVISION,.95)
95%
2551.724
> sd(PROVISION)
[1] 67.57396

Bref, l’écart-type obtenu est ici inférieur à l’erreur de prédiction obtenu auparavant, mais c’est normal car il ne s’agit ici que de mesure l’incertitude sur la prédiction. Pareil pour la formule de Mack,
> library(ChainLadder)
> MackChainLadder(PAID)
MackChainLadder(Triangle = PAID)

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

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

aucun commentaire aucun rétrolien

Les modèles Tweedie en actuariat (1)

Par Arthur Charpentier le mardi 1 décembre 2009, 12:04

Présentation des lois Tweedie

On suppose ici que la fonvtion variance soit de la forme suivante,

Le cas 0 correspond au cas Gaussien, le cas 1 correspond à la loi de Poisson, et le cas 2 à la loi Gamma. L’idée ici est de considérer les cas intermédiaire, et plus précisément entre 1 et 2.
On notera que ce sont les cas classique d’hétéroscédasticité, où l’on suppose que soit que la variance est proportionnelle aux variables explicatives (p=1) soit que l’écart-type est proportionnelle aux variables explicative (que l’on peut rapprocher de p=2).
Exemple: la loi Poisson composées
Considérons le cas où l’on cherche à modéliser une charge totale pour une police d’assurance donnée, où le nombre de sinistre suit une loi de Poisson

et où les coûts de sinistres (individuels) sont une loi Gamma

. Alors la variable composée

vérifie

En posant

alors on peut réécrire les deux premiers moments,

Bref, la loi Poisson composée semble être un cas particulier de loi Tweedie. En posant

de telle sorte que

on retrouve que cette loi est bien dans la famille exponentielle.

Ajustement d’une loi Tweedie

La difficulté est ici de trouver le paramètre puissance. On peut utiliser des méthodes de type maximum de vraisemblance. La "densité" une telle loi serait de la forme suivante
> ftweedie = function(y,p,mu,psi){
+ if(p==2){f = dgamma(y, 1/psi, 1/(psi*mu))} else
+ if(p==1){f = dpois(y/psi, mu/psi)} else
+ {lambda = mu^(2-p)/psi /(2-p)
+ if(y==0){ f = exp(-lambda)} else
+ { alpha = (2-p)/(p-1)
+    beta = 1 / (psi * (p-1) * mu^(p-1))
+    k = max(10, ceiling(lambda + 7*sqrt(lambda)))
+    f = sum(dpois(1:k,lambda) * dgamma(y,alpha*(1:k),beta))
+ }}
+ return(f)
+ }
On peut alors lancer une procédure de type maximum de vraisemblance profilée pour trouver la puissance optimale. On peut utiliser cette fonction pas seulement en tarification, mais aussi en provisionnement. Pour aller un peu plus vite, on peut utiliser la library(statmod), mais pour des aspects numériques, il convient de ne pas avoir de valeurs manquantes dans le triangle,
> source("http://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
> library(statmod)
> an <- 6; ligne = rep(1:an, each=an); colonne = rep(1:an, an)
> passe = (ligne + colonne - 1)<=an; n = sum(passe)
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> Y = as.vector(INC)
> lig = as.factor(ligne)
> col = as.factor(colonne)
> y = Y[passe]
> Y[is.na(Y)]=.01
On utilise alors le code suivant (je reviendrais un jour plus longuement sur les procedures d’affichage de sorties),
> pltweedie <- function(pow){
+ regt = glm(Y~lig+col, tweedie(pow,0))
+ reserve = sum(fitted.values(regt)[!passe])
+ dev = deviance(regt)
+ phi.hat = dev/n
+ mu = fitted.values(regt)[passe]
+ hat.logL = 0
+ for (k in 1:length(y)){
+    hat.logL <- hat.logL + log(ftweedie(y[k], pow, mu[k], phi.hat)) }
+ cat("Puissance =", round(pow,3), "phi =", round(phi.hat,2),
+ "Reserve (tot) =", round(reserve), "logL =", round(hat.logL,3),)
+ hat.logL}
Sur le triangle dont on dispose, on peut regarde la logvraisemblance profilée,
> for(pow in c(1,1.25,1.5,1.75,2)){pltweedie(pow)}
Puissance = 1   phi = 166.95    Reserve (tot) = 1345    logL = -Inf
Puissance = 1.25        phi = 42.92     Reserve (tot) = 1216    logL = -151.72
Puissance = 1.5         phi = 15.8      Reserve (tot) = 996     logL = -145.232
Puissance = 1.75        phi = 9.02      Reserve (tot) = 609     logL = -153.997
Puissance = 2   phi = 6.78      Reserve (tot) = 125     logL = -170.614
Il y a eu 22 avis (utilisez warnings() pour les visionner)
mais aussi lancer une routine d’optimisation,
>  optimize(pltweedie, c(1.01,1.99), tol=1e-4,maximum = TRUE)
Puissance = 1.384       phi = 23.78     Reserve (tot) = 1114    logL = -144.902
Puissance = 1.616       phi = 11.56     Reserve (tot) = 841     logL = -148.172
Puissance = 1.241       phi = 44.76     Reserve (tot) = 1221    logL = -152.711
Puissance = 1.462       phi = 17.9      Reserve (tot) = 1038    logL = -144.765
Puissance = 1.432       phi = 19.89     Reserve (tot) = 1069    logL = -144.625
Puissance = 1.429       phi = 20.04     Reserve (tot) = 1071    logL = -144.624
Puissance = 1.428       phi = 20.17     Reserve (tot) = 1073    logL = -144.624
Puissance = 1.428       phi = 20.15     Reserve (tot) = 1073    logL = -144.624
Puissance = 1.428       phi = 20.15     Reserve (tot) = 1073    logL = -144.624
Puissance = 1.428       phi = 20.16     Reserve (tot) = 1073    logL = -144.624
Puissance = 1.428       phi = 20.15     Reserve (tot) = 1073    logL = -144.624
$maximum
[1] 1.427873

$objective
[1] -144.6237

avec en bleu le montant de provision (total, toutes années de survenance confondues) et en rouge la logvraisemblance profilée (en fonction du niveau du paramètre de la fonction puissance). La suite lorsque j’aurais réussi à récupérer les données sur mon pc qui a crashé la semaine dernière...

aucun commentaire aucun rétrolien

mercredi 18 novembre 2009

Provisions pour sinistres restant à payer (1)

Par Arthur Charpentier le mercredi 18 novembre 2009, 13:13

Suite à la formation de début de semaine, et en préparation du cours de provisionnement non-vie, quelques lignes de code......

Les méthodes classiques de provisionnement

Pour finir (?) sur le provisionnement en assurance non-vie, quelques lignes de code R pour obtenir les principaux estimateurs utilisés en pratique. Tout d’abord l’estimateur Chain Ladder.
On suppose que

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-01.png$

on pose alors, à partir d’estimation des coefficients de transition

$L'image “http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-02.png” ne peut être affichée car elle contient des erreurs.$

> source("http://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
> LAMBDA=rep(NA,6)
> for(i in 1:5){LAMBDA[i]=sum(PAID[1:(6-i),i+1])/sum(PAID[1:(6-i),i])}
> LAMBDA[6]=1DIAGP=diag(PAID[,6:1])
> C.CL=diag(PAID[,6:1])*cumprod(rev(LAMBDA))
L’estimateur de Bornhueter-Ferguson est obtenu en considérant le modèle suivant (en notant 1 la première colonne)

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-03.png$

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-04.png$

On peut alors écrire $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-05.png$ et on suppose $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-06.png$ , i.e. $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-07.png$ . Aussi, $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-08.png$ est une estimation de la charge ultime, et $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-09.png$ est une cadence de paiements (allant de 0 à 1). On appelera estimateur de Bornhuetter-Ferguson l’estimateur

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-11.png$

en notant qu’un estimateur naturel de $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-12.png$ est

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-13.png$

On pourra noter que si

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-14.png$

on retombe sur l’estimateur Chain-Ladder. Il s’agit plutôt de prendre en compte un estimateur exogène de la charge ultime. On parlera d’avis d’expert. L’estimateur le plus utilisé est simplement basé sur un loss-ratio (ratio sinistres à primes) cible. Par exemple, si on suppose un loss ratio cible de 105%, on obtient les estimateurs suivants
> MU=1.05*PREMIUM
> BETA=cumprod(1/rev(LAMBDA))
> C.BF=DIAGP+(1-BETA)*MU
L’estimateur de Benktander-Hovinen est obtenu en considérant une méthode de type crédibilité, i.e. en posant

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-15.png$

On dérive alors la formule suivante

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-16.png$

Une réécriture plus simple donne

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-17.png$

Enfin, l’estimateur de Cape-Cod , tel qu’il a été proposé par Hans Bühlmann (1983), est obtenu en supposant que

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-18.png$

où $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-09.png$ est la cadence de paiement., et où $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-19.png$ est généralement interprété comme la prime acquise pour l’année i. Pour estimer le coefficient $http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-20.png$ on utilise

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-21.png$

Un estimateur global robuste est obtenu en calculant

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-22.png$

Bref, à partir de là, on dérive un estimateur pour le terme diagonal

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/CL-23.png$

Le code R est alors de la forme suivante
> kappa=DIAGP/BETA/PREMIUM
> K=weighted.mean(kappa,w=BETA*PREMIUM)
> C.CC=DIAGP+(1-BETA)*K*PREMIUM

aucun commentaire aucun rétrolien

jeudi 12 novembre 2009

De la pratique du provisionnement (non-vie)

Par Arthur Charpentier le jeudi 12 novembre 2009, 10:09

Quelques illustrations pour le cours sur le provisionnement que nous avons commencé juste avant les vacances dans le cours de statistique de l’actuariat 2. Afin d’illustrer les cadences de paiments, quelques courbes de déroulé, avec sur chaque graphique, le montant des paiements cumulés (en bas), et le montant de la charge estimée par les gestionnaires de sinistres (dite dossier-dossier). Il y a sur tous les graphiques 3 années de survenance représentées (avec forcément moins d’information sur les années les plus récentes). Il ne s’agit pas des montants à proprement parler, mais des montants exprimés en base 100 (la référence étant la charge ultime estimée). Il ne s’agit donc pas d’une méthode pour estimer la charge finale (cela a été fait en amont pour tracer ces graphiques), mais juste illustrer l’idée des cadences de paiements, et surtout de présenter les différences entre les types de risques considérés (car il s’agit de vraies données d’assureur).
Pour commencer, un peu d’assurance automobile de particuliers,

On peut séparer en deux risques assez distincts en automobiles, assurance automobile des dommages matériels,

et l’assurance automobile des dommages coroporels,

Pour les particulier, une autre assurance assez classique est l’assurance multirisque-habitation,

Mais on peut aussi regarder les risques d’entreprise, avec par exemple l’assurance risque industriel,

ou un peu plus compliqué, l’assurance responsabilité civile entreprise,

Un cas particulier qui a beaucoup fait parlé est celui de l’assurance responsabilité civile médicale,

Et pour conclure, un dernier risque relativement compliqué à modéliser, l’assurance construction,

Sur le calcul des provisions dites "pour sinistres restant à payer", je renvoie au code des assurances, en particulier l’article R331-15 (ici par exemple), qui sert de base à toute la modélisation. "La provision pour sinistres à payer est calculée exercice par exercice. Sans préjudice de l’application des règles spécifiques à certaines branches prévues à la présente section, l’évaluation des sinistres connus est effectuée dossier par dossier, le coût d’un dossier comprenant toutes les charges externes individualisables ; elle est augmentée d’une estimation du coût des sinistres survenus mais non déclarés. La provision pour sinistres à payer doit toujours être calculée pour son montant brut, sans tenir compte des recours à exercer ; les recours à recevoir font l’objet d’une évaluation distincte. Par dérogation aux dispositions du deuxième alinéa du présent article, l’entreprise peut, avec l’accord de l’Autorité de contrôle des assurances et des mutuelles, utiliser des méthodes statistiques pour l’estimation des sinistres survenus au cours des deux derniers exercices."

aucun commentaire aucun rétrolien

mercredi 28 octobre 2009

Suite de la tarification a posteriori: les méchanismes bonus-malus

Par Arthur Charpentier le mercredi 28 octobre 2009, 16:29

Les chaînes de Markov

L’outil de base pour l’étude des systèmes bonus-malus est l’utilisation des chaînes de Markov, dès lors que l’on travaille sur des classes de bonus. Pour plus de compléments sur les chaînes de Markov, je renvoie à des vieilles notes de cours (inachevée, et sûrement remplis de coquillettes) ici.
Pour reprendre les grandes idées, on considère une suite de variables aléatoires à valeurs dans $http://perso.univ-rennes1.fr/arthur.charpentier/latex/calE.png$ , $http://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png$ . Cette suite est sera une chaîne de Markov si pour tout $http://perso.univ-rennes1.fr/arthur.charpentier/latex/ngeq1.png$ , et pour tout

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/suitei.png$

telle que

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/markov-proba1.png$

on ait

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/proab1.png$
= $http://perso.univ-rennes1.fr/arthur.charpentier/latex/proab2.png$

On parlera d’absence de mémoire (au sens où seule la dernière valeure intervient dans le calcul de la loi conditionnelle, et pas l’intégralité de la trajectoire). Si $http://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png$ est une chaîne de Markov de distribution initiale $http://perso.univ-rennes1.fr/arthur.charpentier/latex/lambdaV.png$ et de matrice de transition P, pour tout $http://perso.univ-rennes1.fr/arthur.charpentier/latex/nkentier.png$ ,

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/proba-markov-statioonnaire.png$

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/pr135.png$

Pour étudier l’évolution sur le long terme, on peut utiliser l’équation de Chapman-Kolmogorov: si $http://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png$ est une chaîne de Markov de distribution initiale $http://perso.univ-rennes1.fr/arthur.charpentier/latex/lambdaV.png$ et de matrice de transition P, pour tout
$http://perso.univ-rennes1.fr/arthur.charpentier/latex/nkentier.png$ , pour tout $http://perso.univ-rennes1.fr/arthur.charpentier/latex/zvze54v6.png$ ,

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/chapman-kol.png$

Un corrolaire important est que si

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/ezjpjio.png$

alors

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/alcnaoli.png$

Bref, on peut facilement étudier l’évolution à long terme de la chaîne de Markov en faisant un peu de calcul matriciel, i.e. des calculs de puissances de matrices.
En fait, formellement, on peut soit étudier des matrices, soit des graphes. Par exemple, la matrice stochastique suivante
XXXX
peut se traduire par le graph suivant,

Si on veut étudier l’évolution sur deux périodes, il suffit de calculer le carré de la matrice, i.e.

Les systèmes bonus-malus par classes

Les tableaux ci-dessous sont tirés des annexes de la bible de Jean Lemaire, Bonus-Malus systems in automobile insurance. Le plus simple est probablement celui de Hong Kong, qui est Markovien (à l’ordre 1),

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-hong-kong-B-9_m.jpg

L’exemple ci-dessous est celui qui a existé en Allemagne dans les années 80,

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-germany-B-7_m.jpg

puis celui qui existait en Finlande,

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-finnish-B-6_m.jpg

et enfin le système italien du début des années 90,

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-italy-B-11_m.jpg

Pour tous ces sytèmes, une fois spécifiée la loi du nombre de sinistres par an, il est possible de construire la matrice de transition des classes de bonus sur un an, ou plus.
En particulier, l’étude de la loi invariante de ce genre de sytème (que je mentionnais en introduction) permet d’étudier la viabilité du système à plus ou moins long terme, et la vitesse de la convergence vers l’état stationnnaire.

Le système bonus-malus français

En France, le mécanisme de bonus-malus peut paraître un peu plus compliqué, et il fonctionne de la manière suivante.

Le calcul du bonus-malus s’effectue en fonction des sinistres survenus (ou non) pendant la période de douze mois consécutifs précédant de deux mois l’échéance annuelle du contrat d’assurance.
Seuls les accidents responsables donnent lieu à malus. C’est à dire qu’en cas de foce majeur, ou si l’accident est entièrement imputable à un tiers, alors le conducteur n’aura pas de malus. Notons aussi que Le malus est rattaché au contrat d’assurance. Autrement dit, si quelqu’un cause un accident avec un véhicule assuré par un tiers, il se supportera pas le malus.
La prime de référence est la prime qui sera utilisée pour ensuite actualiser le montant de prime en tenant compte du coefficient de majoration ou minoration. Il s’agit de la prime technique calculée par l’assureur en fonction de charactéristiques a priori, comme la puissance du véhicule véhicule, la zone géographique de circulation, l’usage socioprofessionnel ou encore le kilométrage parcouru....
A partir de cette prime de référence, on calcule une prime qui tient compte de l’expérience passée, via un coefficient multiplicatif.
La valeur de référence est de 1 pour le taux de majoration-minoration. Ce taux subie une réduction - corrrespondant à un bonus - pour une année sans accident, d’un niveau de 5 %. Toutefois, la valeur du coefficient est arrêté à la deuxième décimale et arrondi par défaut. Au final, le taux prend donc un nombre fini de valeurs. Pour résumer schématiquement

1 année sans accident	100 x 0.95 = 0.95 de Bonus
2 années sans accident	0.95 x 0.95 = 0.90 de Bonus
3 années sans accident	0.90 x 0.95 = 0.85 de Bonus
4 années sans accident	0.85 x 0.95 = 0.80 de Bonus
5 années sans accident	0.80 x 0.95 = 0.76 de Bonus
6 années sans accident	0.76 x 0.95 = 0.72 de Bonus
7 années sans accident	0.72 x 0.95 = 0.68 de Bonus
8 années sans accident	0.68 x 0.95 = 0.64 de Bonus
9 années sans accident	0.64 x 0.95 = 0.60 de Bonus
10 années sans accident	0.60 x 0.95 = 0.57 de Bonus
11 années sans accident	0.57 x 0.95 = 0.54 de Bonus
12 années sans accident	0.54 x 0.95 = 0.51 de Bonus
13 années sans accident	0.51 x 0.95 = 0.48 plafonné à 0.50

La réduction maximale est donc de 0.50, soit 50% de réduction sur la prime de référence.
En revanche, chaque accident responsable donne lieu à une majoration - on palera de malus - de 25 % du coefficient par accident. Toutefois, cette majoration est réduite de moitié (12,5 %) si la responsabilité n’est que partiellement engagée.Là encore, le coefficient obtenu est arrêté à la deuxième décimale et arrondi par défaut. Aussi, avec deux accidents responsables dans l’année, le coefficient sera multiplié par un coefficient de 1,56. Afin d’éviter l’absence d’assurance, on notera qu’en aucun cas, le coefficient de réduction-majoration ne peut excéder 3,50.
Pour les jeunes conducteurs, l’article A. 335-9-1 du code des assurances prévoit que l’assureur a la possibilité de majorer la prime de référence (dans des limites fixées par la loi là encore).

Notons enfin qu’il existe des avantages pour les bons conducteurs: si un conducteur est au bonus maximal depuis au moins trois ans, il conserve son bonus après le premier accident responsable. Notons que c’est le genre de disposition qui figure dans la loi et que certains assureurs utilisent à des fins commerciales (comme le montre l’affiche ci-dessus ou ci-dessous, qui a été analysée sur le blog de Gizmo (ici)).

Il existe aussi une règle dite de descente rapide qui peut atténuer les effets des malus: après deux années consécutives sans accident, votre coefficient ne peut être supérieur à 1.
Ce qui rend le système intéressant, c’est qu’il est obligatoire, en France. Aussi, en cas de changement d’assureur, le coefficient de bonus-malus sera recalculé à partir d’un relevé d’information: l’assureur est tenu de fournir ce document à la résiliation du contrat (contenant des informations sur le conducteur et le véhicule, ainsi que le nombre des sinistres survenus au cours des cinq années précédant l’établissement du relevé d’information, ainsi que leur nature et la part de responsabilité retenue)
Davantage d’information se trouve ici, là (ce document prpopose d’ailleurs une petite simulation à la fin sur le coût d’un accident) ou là encore (pour une approche plus actuarielle), voire là (pour une approche plus économique).

aucun commentaire aucun rétrolien

jeudi 1 octobre 2009

De l'interprétation d'un effet nonlinéaire en régression

Par Arthur Charpentier le jeudi 1 octobre 2009, 11:10

Je vais poursuivre un peu le TD de statistique de l’actuariat 2, où nous avions parlé tarification et régression Poissonnienne.

Analyser un effet nonlinéaire

Dans un premier temps, nous avions regardé un modèle simple, où l’on voulait voir si l’âge du conducteur expliquait la fréquence de sinistres. Bref, on peut utiliser une régression log-Poisson, qui donne la prédiction suivante, en fonction de l’âge,

Mais l’âge est une variable un peu particulière car c’est une fausse variable continue: seule les valeurs entières apparaissent. Comme on chercher à calculer

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-1.png$

notons qu’il est possible de considérer l’estimateur empirique naturel de cette grandeur, i.e.

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-2.png$

où

$http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-3.png$

Graphiquement, on obtient la courbe suivante,

Le côté ératique de la série à droite et à à gauche vient de la faible représentation dans la base, comme en atteste l’histogramme de l’âge du conducteur,

Je peux en profiter pour revenir à un commentaire posté ici sur les modèles GAM. En effet, au lieu de regarder la fréquence en fonction de l’âge, on peut s’intéresser à la la différence par rapport à la fréquence globale.
> m=mean(base$nbsin)
> Y=predict(REG,newdata=data.frame(age=X,expo=E),type = "response")
> lines(X,Y/m-1,col="red")
On obtient alors le graphique suivant (en superposant à l’ajustement obtenu par GAM)

Notons que si l’on ajuste un modèle log-Poisson uniquement sur la tranche d’âge [50-70], on est très proche de ce que propose le modèle GAM, d’où l’interprétation de l’ajustement GAM comme un ajustement local.

Changer de modèle (et les comparer)

On a choisi un modèle log-Poisson, mais on se doute que d’autres modèles seraient possible. En particulier dans notre base, les assurés ont eu 0 ou 1 sinistre l’année passée. Autrement dit, des lois binomiales peuvent être ajustée,
> REG0 = glm(nbsin~age,base,
+ family=poisson(link="log"),
+ offset=expo)
> REG1a = glm(nbsin~age,base,
+ family=binomial(link = "logit"),
+ offset=expo)
> REG1b = glm(nbsin~age,base,
+ family=binomial(link = "probit"),
+ offset=expo)
> REG1c = glm(nbsin~age,base,
+ family=gaussian(link = "identity"),
+ offset=expo)
Le critère AIC est souvent invoqué pour juger de l’adéquation d’un modèle. Il vaut ici
> AIC(REG0)
[1] 36470.7
> AIC(REG1a)
[1] 36044.45
> AIC(REG1b)
[1] 37304.02
> AIC(REG1c)
[1] 32549.83
Mais on peut également regarder le critère BIC,
> AIC(REG0, k = log(nrow(base)))
[1] 36488.81
> AIC(REG1a, k = log(nrow(base)))
[1] 36062.57
> AIC(REG1b, k = log(nrow(base)))
[1] 37322.14
> AIC(REG1c, k = log(nrow(base)))
[1] 32577
Autrement dit, en utilisant le règle "the smaller, the better" on retiendrait le dernier modèle qui est un modèle Gaussien, manifestement assez mauvais: sur le graphique ci-dessous, la régression log-Poisson est en mauve, presque confondu avec la régression logistique¹, en bleu, la régression probit en vert et le modèle linéaire (Gaussien) en rouge,

Bref, ce dernier modèle n’est pas convainquant... Par contre, notons que si l’on compare le modèle log-Poisson et logit, ce dernier conduit à tarifer moins cher les âges extrêmes,

Autrement dit, au lieu de comparer des modèles actuariels sur la base d’outils statistiques, et uniquement de ces outils, il peut être instructif d’essayer de comprendre qui verra sa prime baisser avec tel ou tel modèle... Mais je reviendrais plus longuement, une autre fois, sur les histoires de choix de modèles.....

aucun commentaire aucun rétrolien

« billets précédents - page 1 de 2

Arthur Charpentier