Lundi 15, Michel Grun Rehomme avait présenté les travaux qu’il avait obtenu avec Noureddine Ben Lagha sur des problèmes d’écrêtement de coûts de sinistres, que l’on peut rapprocher de la détection de points atypiques, ou outliers. Ce matin, Alexander McNeil nous a rappelé que les bagplot pouvaient être très intéressants en gestion des risques (mentionné ici). Je ne peux pas m’empêcher de revenir un petit peu sur cette notion, en m’écartant des applications financières proposées par Alexander pour proposer des applications en démographie et en climatologie (inspiré par des travaux de Rob Hyndman).
- Le boxplot et la régression quantile
Un outil classique (en dimension 1) est le boxplot, proposé par John Tukey en 1977. Par exemple, sur les coûts de sinistres en assurance auto, on obtient
On peut aussi faire un boxplot par classes (ou tranches) d’âge,
qui correspond, en passant à la limite (i.e. en considérant l’âge comme une variable continue), à un ensemble de régressions quantiles.
On peut aussi regarder l’impact de l’ancienneté du véhicule (à la souscription du contrat) sur le coût individuel des sinistres,
avec en rouge la régression (linéaire) du log du coût sur l’âge du véhicule. Là aussi, on peut faire des régressions quantiles,
Le code pour obtenir les régressions quantiles lissées est assez simple, une fois installé les library(quantreg) et library(splines),
> library(quantreg); library(splines) > fit = rq(log(cout)~bs(agecond,df=8),data=BASECOUT,tau=.5) > xu=seq(10,110) > yu=predict(fit,newdata=data.frame(agecond=xu)) > lines(xu,yu,col="red",lwd=2)
Cette notion a été généralisée en multivariée sous le nom de bagplot. Par exemple, pour déterminer la “médiane” en dimension 2, on va recherche un ensemble (appelé bag) contenant 50% des observations.
- La profondeur
Cela se fait à l’aide du concept de profondeur statistique, introduite par John Tukey en 1975. En dimension 1, si
désigne la fonction de répartition, la “profondeur de demi-espace” s’écrit
qui est maximal pour la médiane (on obtient alors 
). Plus généralement, pour un quantile d’ordre 
,
Cette notion se généralise en dimension supérieure en posant
On notera que Fraiman et Muñiz ont proposé une extension en 2001. S’il n’y a pas d’ex-aequos, en dimension 1,
qui appartient à 
, et dont la médiane est la plus profonde. Liu a proposé une profondeur dite simpliciale en 1984,
D’un point de vue statistique, la fonction de répartition est inconnue, mais on dispose de 
. On note alors 
l’équivalent empirique de la profondeur. En dimension 2, ça sera 
où 
. La méthode de construction est détaillée ci-dessous,
- Construction d’un bagplot
Pour cela, on note 
le plus petit nombre d’observations contenu dans un demi plan passant par 
. On note 
l’ensemble des points tels que 
contienne au moins 
points. On notera qu’on peut montrer que 
, et que ces régions sont convexes. Rousseeuw et Ruts ont proposé un algorithme permettant de construire un tel ensemble de région en 1998. Sur l’exemple suivant, on considère 20 points en dimension 2,
Notons qu’en changeant de méthode, on obtient des résultats qui peuvent être très changeant.
> library(mnormt) > Z=rmnorm(20,mean=c(0,0),varcov=matrix(c(1,.6,.6,1),2,2)) > source("http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/bagplot/bagplot.R") > bagplot(Z,factor=1,dkmethod=1) > bagplot(Z,factor=1,dkmethod=2)
La méthodologie a été proposée par Rousseeuw, Ruts et Tukey en 1999 (ici).
- Copules Gaussiennes, versus copules de Clayton
On peut essayer de voir si les bags changent de forme lorsque la copule sous-jacente change. L’exemple suivant montre le bagplot d’un vecteur gaussien,
alors que l’exemple suivant montre le bagplot pour une copule duale de Clayton, avec des marges Gaussiennes.
Dans le premier cas, on retrouve que les bagplots ressemblent des ellipses (ce qui correspond à ce que j’avais raconté ici sur les ellipses de confiance). On retrouve cette ellipse sur le graphique dessous, avec 5000 points (simulés) au lieu de 200,
pour la copule Gaussienne, et
pour la copule de Clayton. La petite animation ci-dessous montre l’évolution des bag en fonction de la corrélation du vecteur Gaussien,
- Illustration sur les données Loss-ALAE
Considérons la base (célèbre) des coûts de sinistres, et des frais associés.
> library(evd); data(lossalae)
La figure ci-contre montre l’échantillon de 1500 points, avec le logarithme des coûts, et le logarithme des frais. Classiquement, on s’intéresse à la loi conditionnelle des frais sachant le coût du sinistre. On peut alors tracer des boxplot conditionnels (comme pour le coût des sinistres conditionnels à l’âge du véhicule),
Mais l’idée des bagplot est sensemblement différentes. On s’intéresse véritablement au couple de variables.
- Du bagplot fonctionnel
Rob Hyndman a également utilisé ces méthode de bagplot quand la taille devient infiniment grande, i.e. dans un cadre fonctionnel. Par exemple, on peut utiliser cette technique pour détecter des années atypiques sur des évolutions de température, ou dans un ajustement de modèle de type Lee-Carter en mortalité. Ces deux cas sont abordées par Rob Hyndman dans la library(rainbow). Pour faire simple et rapide, on peut définir une profondeur fonctionnelle,
Pour commencer, on peut regarder les température moyenne mensuelles de la température à la surface de la mer, entre 1951 et 2007, obtenue ici.
On note que 4 années se détachent des autres, 1982/1983 et 1997/1998, qui sont particulièrement élevées, mais pas seulement. Certaines semblent décallées de quelques mois. Le code est ici simplement
> library(rainbow) > fboxplot(data = ElNino, plot.type = "functional", type = "bag") > fboxplot(data = ElNino, plot.type = "bivariate", type = "bag")
Dans le cas des données démographiques, il faut bricoler un peu. Repartons des données français que nous avons étudiées plus d’une fois,
> DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE) > 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)) > DEATH$Age[is.na(DEATH$Age)]=110 > EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age)) > EXPOSURE$Age[is.na(EXPOSURE$Age)]=110 > YEAR=unique(DEATH$Year);nC=length(YEAR) > AGE =unique(DEATH$Age);nL=length(AGE) > rownames(MUH)=AGE > colnames(MUH)=YEAR > rownames(MUF)=AGE > colnames(MUF)=YEAR > MUH=MUH[1:90,] > MUF=MUF[1:90,] > MUHF=fts(x = AGE[1:90], y = log(MUH), xname = "Age", + yname = "Log Mortality Rate") > MUFF=fts(x = AGE[1:90], y = log(MUF), xname = "Age", + yname = "Log Mortality Rate") > fboxplot(data = MUHF, plot.type = "functional", type = "bag") > fboxplot(data = MUHF, plot.type = "bivariate", type = "bag")
avec un code similaire pour les femmes. Le graphique suivant montre l’évolution des taux de mortalité des hommes en fonctions de leur âge, pour différentes années,
On observe que certaines années se détachent assez nettement, en particulier sur les âges allant de 15 à 40 ans. Comme le montre le graphique à gauche, ce sont les années de guerre qui ressortent fortement (surtout la première guerre mondiale). Pour les femmes, l’évolution des taux de mortalités en fonction des âges semble en revanche beaucoup plus régulière, puisqu’aucune année atypique ne semble ressortir,
Ceci se confirme sur le graphique de droite, où, contrairement aux hommes, aucune année ne semble réellement se détacher. Notons qu’au lieu de chercher des années atypiques, on peut chercher faire une lecture duale, en cherchant les âges qui se ont eut un changement en terme d’évolution du taux de mortalité sur un siècle différent. Par exemple pour les femmes, les taux de mortalité ont l’évolution suivante au cours du siècle,
ce qui permet de repérer quelques âges un peu atypiques, i.e. les bébés, qui ont eu un gain en terme de taux de mortalité considérable, relativement aux autres âges. En revanche, pour les hommes, il n’y a pas particulièrement d’âge qui soit statistiquement différent des autres,
- Applications en terme de stress scenarios
Nous avons vu qu’il y avait presque toujours (pour les hommes comme pour les femmes) des outliers, i.e. des courbes atypiques. Il serait donc normal de s’attendre à en trouver dans les scénarios de prédiction. Les bagplot sont alors un outils intéressant. Si on regarde les prédictions obtenues par la méthode de Lee & Carter par exemple, ou plutôt des scénarios générés par la méthode de Lee & Carter,
> LCH <- fdm(BASEH,order=10) > LCHf=forecast(LCH,100,method="ets") > LCHg=simulate(LCHf,1) > LCHg1=as.matrix(LCHg[,,1]) > MUHf=LCHg1[1:90,] > MUHF=fts(x = AGE[1:90], y = log(MUHf), xname = "Age", + yname = "Log Mortality Rate") > library(rainbow) > fboxplot(data = MUHF, plot.type = "functional", type = "bag") > fboxplot(data = MUHF, plot.type = "bivariate", type = "bag")
On peut alors identifier clairement des années atypiques – pour un scenario donné. Bref, on a un outils pour vérifier que l’on génère effectivement des points atypiques, un peu extrêmes, mais surtout, on peut les identifier. On peut alors creuser un peu pour mieux comprendre ce qui se passe….
Pour aller plus loin sur ces histoires de bagplot et aux régressions quantiles, je renvoie ici ou là (Thierry Roncalli avait remis la régression quantile au goût du jour il y a quelques années, à l’époque où je programmais encore en Gauss, je renvoie aussi ici pour des applications financières).
- Les quantiles bivariés
En fait, le point de départ du billet était très lié d’une notion qui méritera un autre billet: comment définir des quantiles (et plus généralement des mesures de risques) dans un cadre multivarié….
, i.e. minimiser la distance entre ma boule et le cochonnet. On appelle ça 
, i.e. maximiser la distance entre la boule de mon adversaire et le cochonnet. On appelle ça 
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: 
and 
. From Fréchet-Hoeffding bounds, we know that
. Further, correlation matrices are necessarily positive-semidefinite matrices.
will be a function of 
.
and 
) are functions of 
,
and 
.
, respectively describing global dependence (that can be related to 
from 
(but actually any positive distribution should work) and 
)
with a Bernoulli 
distribution
from copula 
if 
and 
if 
and 
, and 
In order to get a better understanding, I look at slices, for instance when 
. 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 
and 
.
The paper on the 2003 heat wave has been submitted a few months ago, and now, the preprint can be downloaded (
Toujours pour répondre à une question posée par mail (“ça quoi ça sert le test lillie sous R ? “), un court billet. Avant de parler de Lillie1, il faut revenir sur Kolmogorov-Smirnov. L’idée est de tester
converge en loi (si l’hypothèse H0 était la bonne) vers un pont brownien changé de temps par la fonction quantile 
.
ici). En pratique, on ne connaît pas cette loi, mais on peut l’estimer, en particulier pour les lois paramétriques. On peut alors chercher à tester autre chose, comme
donne
donne
dépend du niveau des taux d’intérêt. L’ensemble des tables de décallage sont en ligne 
si les deux époux sont en vie, soit en payant une prime 
si l’un des époux décède. En contrepartie, l’assureur verse 10 000 euros si les deux époux sont en vie, et 60% de ce montant au conjoint survivant.
, et 
(i.e. les deux membres du couple sont en vie à la signature du contrat).
(les autres découlant de cette dernière). Pour importer les tables de mortalité, on pourra utiliser
