# Chaire ActInfo dans l’Argus

Passionnante interview dans l’Argus de l’Assurance, à paraître demain, sur la chaire actinfo.

# MultiAttribute Copula Utility Functions

In June, with Olivier L’Haridon, we will organize a (small) conference, in Rennes, on risk models in a multi-attribute framework. In order to fully enjoy the workshop (more to come on the blog), we organized this year an internal workshop on that topic. A gave an oveview in September on multivariate distributions, with an emphasis on spherical / elliptical distributions, distributions on the simplex, and copulas. This time, following recent presentations made by Olivier, I will present Ali E. Abbas (recent) contributions on copula-type multriattribute utility functions. Slides are online, and the presentation will be this Thursday

As discussed in the introduction, one (nice) application can be the choice of a seat in a theatre, see

# Radial Graphs for Time Series

On How to: Weather Radials, there was a nice visualisation of temperatures. Since I am too old fashioned for ggplot2, I wanted to reproduce a similar graph with the old plot style.

Assume that daily temperature is in a vector X (e.g. temperature in Montréal, QC, in 2009). To get a radial plot, use

> n=length(X)
> theta=seq(0,1-1/n,length=n)*2*pi
> r=30+X
> plot(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",xlab="",ylab="",axes=FALSE)
> for(t in 1:n){
+   if(X[t]>0) CL=rgb(0,0,1,.4)
+   if(X[t]<0) CL=rgb(1,0,0,.4)
+   if(X[t]==0) CL="white"
+   segments((30+X[t])*cos(pi/2-theta[t]),(30+X[t])*sin(pi/2-theta[t]),30*cos(pi/2-theta[t]),30*sin(pi/2-theta[t]),col=CL)
+ }
> for(r in 10*seq(0,6)) lines(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",col="light blue")

# Prix Innovation Formation Digitale en Assurance

Ce jeudi, un prix “Innovation Formation Digitale en Assurance” sera décerné à la formation Actuaire Data Science de l’Institut du Risk Management, dont on est directeurs des études avec Romuald Elie et Jérémie Jakubowicz. Comme on a pas mal trimé pour monter le dossier, on va se lancer des fleurs… Le prix sera remis par Arnaud Chneiweiss, aka @achneiweiss, secrétaire général du GEMA et délégué général de l’Association Française de l’Assurance (AFA) dans le cadre d’un salon de l’enseignement supérieur et de la recherche. Allez, vive nous !

Dans 10 jours, la seconde promo attaque sa deuxième séance de cours (on avancera sur les notions plus avancées de R, avec les bases de données et les graphs, avant de faire les modèle de régression, le mois prochain) et dans 20 jours, ça sera les soutenances des mémoires de la première promo….

# Les Données et le Monde Libre, et Ouvert

La période des conférences va commencer… Même si j’ai du en refuser – puis en annuler – un bon nombre (pour des raisons personnelles), je serais cette semaine au Breizh Camp devant 500 développeurs, pour parler données, algorithmes et économie. Dans deux semaines, je ferais deux interventions à Barcelone, une sur les données et l’autre sur les algorithmes (dans un optique plus assurantielle). Et en mai, je suis invité à la Banque d’Angleterre pour parler au Chief Economists’ workshop – What can central bank policy- makers learn from other disciplines?. J’ai été invité avec Michael Lo et Andrew Gelman (excusez du peu) pour parler un peu de données, de recherche, d’algorithmes, etc. Et au milieu de ces conférences un peu moins pointues que les conférences (très) spécialisées que j’ai l’habitude de faire, pour présenter mes papiers de recherche devant mes pairs, j’ai également été invité à parler dans un contexte assez singulier, sur l’économie numérique,

Oui, dans le cadre de rendez-vous parlementaires, autrement dit, avec des élus. A priori très loin de mon auditoire traditionnel.

On m’a demandé de faire un intervention sur la donnée comme pierre angulaire de l’économie. J’avais accepté car j’avais entendu “économie” comme “sciences économiques“, et je trouvais la question pertinente, et digne d’intérêt. Malheureusement, en réalité, le débat portera davantage sur la valeur de la donnée, et probablement sur les données ouvertes. Comme c’est un autre sujet qui m’intéresse, je ne me suis pas dégonflé (même si je pense que mon voisin Simon aurait été plus légitime…)

Dans mon exposé pour Breizh Camp, j’évoque – à la toute fin – les données ouvertes, en expliquant – en reprenant l’exemple du paradoxe de Simpson – que les données n’ont de sens que si elles sont individuelles, et pas agrégées ! Bref, la majorité des données – dites ouvertes – que l’on trouve en France sont inutiles, car pas assez fines ! Inutile dans l’idée d’expliquer des comportements et des stratégies individuelles. Pour expliquer l’obésité, il faut des données individuelles, pas par département !

Mais ce n’est pas nouveau. Lorsque je travaillais sur le climat, il y a une quinzaine d’années, avoir des données avec une faible granularité était impossible : j’ai travaillé sur la modélisation des tempêtes avec des données irlandaises et néerlandaises, pour les inondations, on a travaillé sur des données canadiennes. Lorsque je travaillais sur la période de retour de la canicule telle qu’elle a été observée pendant l’été 2003, j’ai eu la chance de trouver des données (journalières) françaises sur un site néerlandais. Et lorsqu’on a travaillé avec Amadou sur le paradoxe épidémiologique sur le tabagisme et le poids des bébés à la naissance, on a travaillé sur les données américaines du CDC, où plusieurs millions de données individuelles sont disponibles ! Bref, des données utiles pour faire de la recherche !

On commence à trouver des données intéressantes, en libre accès en France. La France a toujours eu des années de retard, en matière de transparence et d’accessibilité de données, j’en ai parlé à de nombreuses reprises sur mon blog. Je repense à la dizaine des courriels que j’ai du envoyer pour récupérer, auprès du Rectorat de l’Académie de Rennes, les statistiques de réussite au brevet des collèges sur 10 ans, pour seulement 4 collèges (pour un TD de statistique avec mes étudiants) ! Le pire c’est que je me dits “mais tout le monde le sait que c’est impossible d’obtenir des données en France, j’enfonce des portes ouvertes là“. Mais je suis constamment surpris de voir que les portes ne sont pas si ouvertes que ça. Pire, que certains tentent de les refermer.

En fait, j’ai l’impression que le débat sur l’ouverture des données est proche de celui sur l’ouverture de la publication scientifique, et sur les logiciels libres. L’idée sous-jacente est que tout le monde, chercheur institutionnel ou simple citoyen, peut avoir envie de comprendre le monde qui l’entoure. C’est tout du moins sur ce postulat que j’ai lancé mon blog voilà presque dix ans. Pour ma part, j’ai besoin de trois types d’objets : des données, sur lesquelles je teste des hypothèses, des modèles économétriques qui permettent de décrire aussi simplement que possible un phénomène complexe, et des packages informatiques. Plus généralement, on peut parler d’un accès ouvert aux connaissances, pour reprendre le titre de la tribune du Monde,

Plusieurs tribunes ont été écrites sur le sujet . Lorsque j’ai publié mon billet sur l’importance des logiciels libres, j’expliquais que ça partait d’un constat pragmatique, et aucunement d’une forme de militantisme (loin de moi l’idée de militer pour quelque cause que ce soit…). Et je dois avouer que les tribunes qui vont en sens inverse me fatiguent, je pense en particulier à la dernière tribune de Philippe Minard (surtout quand il prétend être la voix “des chercheurs“). Personnellement, je ne peux faire de la recherche que si je peux avoir accès à des données aussi fines que possible (et oui, souvent des données individuelles, ce qui pose le problème de l’anonymisation, on pourra en reparler un autre jour), si j’ai des logiciels et des codes informatiques que je peux partager avec les co-auteurs (voire mettre en ligne sur mon blog pour en faire profiter tous ceux qui le souhaitent), et si je peux comprendre ce qui a été écrit sur le sujet, c’est à dire avoir accès aux articles qu’il faudrait que le lise, et que je pourrais alors citer. Les modèles économiques des logiciels, des données et de la publication scientifique doivent changer ! Les expériences que l’on a pu observer dans le domaine des logiciels montrent que ceux qui refusent de s’adapter disparaîtront… Personnellement, les articles que je lis, je les trouve via arxiv ou hal (et je suis toujours choqué quand un éditeur de demande de citer la version publiée, alors que souvent je ne l’ai pas lue, ce qui peut causer des soucis quand la version publiée est assez différente… je repense à un erratum que j’ai du publier il y a quelques années). Et quand je n’ai pas les données que je veux, avec un peu de chance, avec un peu de webscrapping, on peut les extraire malgré tout….

Bref, je pense que mes prochains exposés vont continuer à ouvrir ces portes ouvertes, car il se serait dommage que quelqu’un ai la mauvaise idée de les refermer… Et en ce sens, je rejoins pleinement ce qui a été écrit dans Pour une Science Ouverte à Tous, publié ce mois-ci, dans Le Monde.

# Classification on the German Credit Database

In our data science course, this morning, we’ve use random forrest to improve prediction on the German Credit Dataset. The dataset is

> url="http://freakonometrics.free.fr/german_credit.csv"
> credit=read.csv(url, header = TRUE, sep = ",")

Almost all variables are treated a numeric, but actually, most of them are factors,

> str(credit)
'data.frame':	1000 obs. of  21 variables:
$Creditability : int 1 1 1 1 1 1 1 1 1 1 ...$ Account.Balance : int  1 1 2 1 1 1 1 1 4 2 ...
$Duration : int 18 9 12 12 12 10 8 ...$ Purpose         : int  2 0 9 0 0 0 0 0 3 3 ...

(etc). Let us convert categorical variables as factors,

> F=c(1,2,4,5,7,8,9,10,11,12,13,15,16,17,18,19,20)
> for(i in F) credit[,i]=as.factor(credit[,i])

Let us now create our training/calibration and validation/testing datasets, with proportion 1/3-2/3

> i_test=sample(1:nrow(credit),size=333)
> i_calibration=(1:nrow(credit))[-i_test]

The first model we can fit is a logistic regression, on selected covariates

> LogisticModel <- glm(Creditability ~ Account.Balance + Payment.Status.of.Previous.Credit + Purpose +
Length.of.current.employment +
Sex...Marital.Status, family=binomial,
data = credit[i_calibration,])

Based on that model, it is possible to draw the ROC curve, and to compute the AUC (on ne validation dataset)

> fitLog <- predict(LogisticModel,type="response",
+                   newdata=credit[i_test,])
> library(ROCR)
> pred = prediction( fitLog, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCLog1=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCLog1,"\n") AUC: 0.7340997 An alternative is to consider a logistic regression on all explanatory variables > LogisticModel <- glm(Creditability ~ ., + family=binomial, + data = credit[i_calibration,]) We might overfit, here, and we should observe that on the ROC curve > fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) > pred = prediction( fitLog, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCLog2=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCLog2,"\n")
AUC:  0.7609792

There is a slight improvement here,  compared with the previous model, where only five explanatory variables were considered.

Consider now some regression tree (on all covariates)

> library(rpart)
> ArbreModel <- rpart(Creditability ~ .,
+  data = credit[i_calibration,])

We can visualize the tree using

> library(rpart.plot)
> prp(ArbreModel,type=2,extra=1)

The ROC curve for that model is

> fitArbre <- predict(ArbreModel,
+                     newdata=credit[i_test,],
+                     type="prob")[,2]
> pred = prediction( fitArbre, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCArbre=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCArbre,"\n") AUC: 0.7100323 As expected, a single has a lower performance, compared with a logistic regression. And a natural idea is to grow several trees using some boostrap procedure, and then to agregate those predictions. > library(randomForest) > RF <- randomForest(Creditability ~ ., + data = credit[i_calibration,]) > fitForet <- predict(RF, + newdata=credit[i_test,], + type="prob")[,2] > pred = prediction( fitForet, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCRF=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCRF,"\n")
AUC:  0.7682367

Here this model is (slightly) better than the logistic regression. Actually, if we create many training/validation samples, and compare the AUC, we can observe that – on average – random forests perform better than logistic regressions,

> AUC=function(i){
+   set.seed(i)
+   i_test=sample(1:nrow(credit),size=333)
+   i_calibration=(1:nrow(credit))[-i_test]
+   LogisticModel <- glm(Creditability ~ .,
+    family=binomial,
+    data = credit[i_calibration,])
+   summary(LogisticModel)
+   fitLog <- predict(LogisticModel,type="response",
+                     newdata=credit[i_test,])
+   library(ROCR)
+   pred = prediction( fitLog, credit$Creditability[i_test]) + AUCLog2=performance(pred, measure = "auc")@y.values[[1]] + RF <- randomForest(Creditability ~ ., + data = credit[i_calibration,]) + fitForet <- predict(RF, + newdata=credit[i_test,], + type="prob")[,2] + pred = prediction( fitForet, credit$Creditability[i_test])
+   AUCRF=performance(pred, measure = "auc")@y.values[[1]]
+   return(c(AUCLog2,AUCRF))
+ }
> A=Vectorize(AUC)(1:200)
> plot(t(A))

# Prévision de l’intérêt pour le chocolat, avec des modèles ARIMA

Hier en cours, on s’était lancé sur la prévision de la série obtenue sur Google Trends pour le mot clé “chocolat“. La série est la suivante

> url = "http://freakonometrics.free.fr/report-chocolat.csv"
> X=report$V2 > plot(X,type="l") Comme on l’avait évoqué, c’est plus simple de construire un modèle ARIMA sur des séries mensuelles qu’hebdomadaire. Donc on va mensualiser la série. Ah oui, on va aussi enlever les deux dernières années pour faire une comparaison de modèles (en confrontant les prévision avec les observations) > date1=substr(as.character(report$V1),1,10)
> date2=substr(as.character(report$V1),14,23) > d1=as.Date(date1,"%Y-%m-%d") > d2=as.Date(date2,"%Y-%m-%d") > Vect_J=Vect_M=Vect_X=Vect_A= NULL > for(i in 1:length(d1)){ + jours = seq(d1[i],d2[i],length=7) + mois = as.POSIXlt(jours)$mon+1
+   annees = as.POSIXlt(jours)$year+1900 + donnees = rep(X[i],7) + Vect_J = c(Vect_J,jours) + Vect_M = c(Vect_M,mois) + Vect_A = c(Vect_A,annees) + Vect_X = c(Vect_X,donnees) + } > base = data.frame(J = Vect_J,M = Vect_M, + A = Vect_A,X = Vect_X) > base$AM = base$A*100+base$M
> Y = tapply(base$X,as.factor(base$AM),mean)
> Z = ts(as.numeric(Y[1:(146-24)]),
+ start=c(2004,1),frequency=12)
> plot(Z)

Compte tenu de l’amplitude croissante des cycles, il peut être intéressant, ici, de passer au logarithme de la série. Et comme on observe une rupture dans la tendance (linéaire au début, et relativement constante ensuite), on peut ajuste une régression linéaire avec une rupture (mais continue)

> plot(as.numeric(Z),type="l")
> X=log(as.numeric(Z))
> T=1:length(X)
> db=data.frame(X,T)
> trend=lm(X~T+I((T-80)*(T>80)),data=db)
> plot(X,type="l")
> lines(predict(trend),col="red")

C’est le résidu de cette série qu’on va chercher à modéliser,

> Y=residuals(trend)
> plot(Y,type="l")

Pour modéliser cette série, on peut commencer par regarder sa série d’autocorrélation

>  plot(acf(Y,lag=36),lwd=5)

et d’autocorrélation partielles,

> plot(pacf(Y,lag=36),lwd=5)

La série est stationnaire, mais il y a une composante cyclique forte (peut être une racine unité saisonnière ?). On peut tenter un premier modèle AR, ou ARMA (avec un AR les résidus ne sont pas un bruit blanc).

> model1 = arima(Y,order=c(12,0,12),
+ seasonal = list(order = c(0, 0, 0), period = 12))

La série des résidus est ici un bruit blanc

> E=residuals(model1)
> plot(acf(E,lag=36),lwd=5)

On peut alors utiliser ce modèle pour faire des prévisions, sur la série initiale (avant de prendre le logarithme, et d’ajuster une tendance)

> futur=length(X)+1:24
> Yp=predict(model1,n.ahead=24)$pred+ + predict(trend,newdata=data.frame(T=futur)) > Xp=exp(as.numeric( + Yp+.5*predict(model1,n.ahead=24)$se^2))
> plot(as.numeric(Z),type="l",xlim=c(1,150))
> lines(futur,Xp,col="red")
> Y2=tapply(base$X,as.factor(base$AM),mean)
> obs_reel=as.numeric(Y2[(146-24)+1:24])
> lines(futur,obs_reel,col="blue")    

avec en rouge notre modèle, et en bleu les vraies observations. On peut alors calculer la somme des carrés des erreurs, sur ces 24 observations,

> sum( (obs_reel-Xp)^2 )
[1] 190.9722

Mais on peut tenter d’autres modèles, par exemple en changeant la tendance (on ne va pas le faire) ou en changeant le modèle ARIMA (en intégrant une racine unité saisonnière)

> model3=arima(Y,order=c(0,0,2),
+ seasonal = list(order = c(0, 1, 1), period = 12))
> E=residuals(model3)
> model3

Coefficients:
ma1     ma2     sma1
0.2246  0.3034  -0.9999
s.e.  0.0902  0.0925   0.3503

sigma^2 estimated as 0.002842:  log likelihood = 152.37,  aic = -296.75

On vérifie que les résidus sont effectivement un bruit blanc

> plot(acf(E,lag=36),lwd=5)
> Box.test(E,lag=12)

Box-Pierce test

data:  E
X-squared = 6.326, df = 12, p-value = 0.8988

On peut alors faire une prévision sur notre série brute,

> futur=length(X)+1:24
> Yp=predict(model3,n.ahead=24)$pred+ + predict(trend,newdata=data.frame(T=futur)) > Xp=exp(as.numeric( + Yp+.5*predict(model3,n.ahead=24)$se^2))
> plot(as.numeric(Z),type="l",xlim=c(1,150))
> lines(futur,as.numeric(Xp),col="red")
> Y2=tapply(base$X,as.factor(base$AM),mean)
> obs_reel=as.numeric(Y2[(146-24)+1:24])
> lines(futur,obs_reel,col="blue")

La somme des carrés des erreurs est ici un peu plus faible

> sum( (obs_reel-Xp)^2 )
[1] 173.8138

C’est à dire que ce second modèle est meilleur que le précédant, sur les deux dernières années. C’est probablement celui que l’on devrait utiliser pour faire des prévisions pour les années à venir (à moins qu’en continuant l’analyse on n’en trouve un meilleur).

# Forecasts with ARIMA Models

In our time series class this morning, I was discussing forecasts with ARIMA Models. Consider some simple stationnary AR(1) simulated time series

> n=95
> set.seed(1)
> E=rnorm(n)
> X=rep(0,n)
> phi=.85
> for(t in 2:n) X[t]=phi*X[t-1]+E[t]
> plot(X,type="l")

If we fit an AR(1) model,

> model=arima(X,order=c(1,0,0),
+             include.mean = FALSE)
> plot(P$pred) > lines(P$pred+2*P$se,col="red") > lines(P$pred-2*P$se,col="red") > abline(h=0,lty=2) > abline(h=2*P$se[20],lty=2,col="red")
> abline(h=-2*P$se[20],lty=2,col="red") we observe the exponential decay of the forecast towards 0, and the increasing confidence interval (where the variance increases, from the variance of the white noise to the variance of the stationnary time series). Plain lines are conditional forecast (given our latest observation, since the AR(1) is a first order Markov process), and dotted lines are unconditional. Let us store some values, to use them as benchmark > s=P$se[20]
> y=P$pred If we fit a MA(1) model > model=arima(X,order=c(0,0,1), + include.mean = FALSE) > P=predict(model,n.ahead=20) > plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

after two lags, the forecast is null, and the (conditional) variance remains constant. But if we consider a moving average process with a longer order,

> model=arima(X,order=c(0,0,14),
+             include.mean = FALSE)
> plot(P$pred) > lines(P$pred+2*P$se,col="red") > lines(P$pred-2*P$se,col="red") > abline(h=0,lty=2) > abline(h=2*s,lty=2,col="red") > abline(h=-2*s,lty=2,col="red") > lines(y,col="grey") we get an output that can be compared with the AR(1) processes. Which makes sense since our AR(1) process can also be seen as a MA(∞), with infinite order. But if we think that our time series is not stationary, an we fit an integrated model > model=arima(X,order=c(0,1,0), + include.mean = FALSE) > P=predict(model,n.ahead=20) > plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

we observe the (standard) martingale property: the forecast is flat, and the confidence interval keeps increasing, and actually, the variance increases towards infinity (at a linear rate). So one should be very careful when differentiation a time series… it will have a huge impact on the forecasts….

# Changement Climatique et Assurance

L’article co-rédigé avec AnneAlexis et Julien sur le changement climatique et l’assurance est paru dans le dernier numéro de Variances.

# Breizh Camp, Economics with Computers

I have been invited, as keynote speaker, at the 6th BreizhCamp, organized in Rennes, from March 23rd till March 25th, “la conférence des développeurs du grand ouest” as they call it. I am deeply honored, since it is a huge conference… The organizing committee asked me to give a (brief) talk on data, and big data. But data is just the visible tip of the iceberg, and I cannot talk about data without mentioning algorithms. So I will try to talk about algorithmics, econometrics, machine learning, and data (and big data, of course).

Slides are now online… More to come in the next days…

# Bonferroni, les anniversaires et les terroristes (ou pas)

Le paradoxe des anniversaires est un exercice de probabilité que l’on découvre assez vite. Mais assez souvent, on l’oublie, alors qu’il est en réalité assez profond. Le problème est le suivant: quel est le nombre de personnes qu’il faut réunir pour avoir (au moins) une chance sur 2 pour que deux personnes aient leur anniversaire le même jour. Si on calcule la probabilité que deux personnes dans un groupe de n fêtent leur anniversaire le même jour, on peut calculer la probabilité que personne ne soit né le même jour, parmi n personnes. Un calcul rapide montre qu’avec 23 personnes, la probabilité qu’au moins deux personnes fête son anniversaire le même jour est 50.73%. Autrement dit, il suffit d’une groupe de 23 personnes (c’était la taille d’une classe d’élèves, du moins à mon époque…).

L’intuition, pour comprendre pourquoi nous avons l’impression que ce nombre est petit est que, naturellement, si on a un événement qui survient avec probabilité p, et qu’on répète l’expérience n fois, la probabilité que rien ne se produise est (1-p)n. Aussi, en première approximation (et j’insiste sur ce point parce que c’est avec ce genre d’argument fallacieux qu’on conclue que la probabilité d’avoir une catastrophe nucléaire dépasse 100%, conduisant à du “certitude statistique“), on en déduit que la probabilité qu’un événement survienne est de l’ordre de np. Sauf qu’il, on compte des paires. Et avec n personnes, on constitue environ n2/2 paires. Et donc la probabilité devient proche de n2p/2… ce qui est bien plus grand que np. C’est grossièrement ce qui explique que notre intuition se sent trompée : avec seulement 23 personnes, on peut constituer plus de 250 paires. La probabilité que toutes les paires correspondent à des dates d’anniversaire différente n’est pas si petit qu’il y parait.

Jeff Ullman parle ainsi de principe de Bonferroni,

“if you look in more places for interesting patterns than your amount of data will support, you are bound to find crap”

Et cette intuition faussée quand on travaille sur des paires se retrouve dans les techniques utilisées par les services de renseigngement, tel que le montre Anand Rajaraman dans Mining Massive DatasetsSupposons qu’il existe un milliard d’individus qui pourraient potentiellement être des terroristes. Que les gens, quels qu’il soient vont à l’hôtel 1 jour sur 100. Un hôtel donné peut accueillir 100 personnes. On suppose qu’il existe 100,000 personnes, de sorte qu’en moyenne ils soient plein tous les jours. Et supposons qu’on regarde sur un historique de 3 ans, i.e. 1000 jours pour avoir des calculs ronds. Supposons que pour traquer les terroristes, on cherche les gens qui ont séjourné, ensemble, dans un hôtel, le même jour, au moins deux fois… On imagine que c’est effectivement rare, et donc suspicieux…

La probabilité que deux personnes données dorment à l’hôtel le même jour est faible, 1/100×1/100, soit 0,0001. Et comme il y a 100,000 hôtel, la probabilité que les deux personnes descendent, en plus, dans le même hôtel, il faut diviser 0,0001 par 100,000. Ce qui donne 10-9 . Mais ça, c’est la probabilité pour qu’ils se croisent une fois ! Si on regarde la probabilité qu’ils se croisent deux fois, en trois ans, on arrive à une probabilité de 10-18 Pour ceux qui ont du mal, ça doit faire environ

0,00000000000000001 %

Et encore. Ici il peut en plus se croiser la seconde fois dans un autre hôtel. La probabilité qu’il se croisent dans le même hôtel que la première fois devient encore plus petite !

Supposons maintenant qu’un organisme de surveillance avisé décider de lister les paires de personnes qui, sur des paires de jours, se sont retrouvé dans le même hôtel. Ça doit être suffisamment rare pour être suspicieux. Rappelons que quand on considère n individus, le nombre de paires possibles est de l’ordre de n2/2. Autrement dit, avec un million de personnes, le nombre de paires est de l’ordre de 5×1017. Et pareil, le nombre de paires de jours est de l’ordre de 5×105. Aussi, le nombre d’événements suspicieux que l’on observe, sur ces trois années est de l’ordre de

(5×1017 ) x (5×10) x10-18 = 250000

Autrement dit, 1/4 de millions de paires de personnes sont des terroristes potentiels ! Soit 1/2 millions d’individus.

Maintenant, supposons qu’il y ait 10 terroristes. S’ils sont assez stupides pour se retrouver dans les mêmes hôtels, les mêmes jours, ça veut dire que sur les 500000 personnes qui vont être surveillées, 10 sont effectivement des terroristes, et on a 499990 faux positifs… Enfin, “faux positifs”… ils ont forcément quelque chose d’autre à se reprocher, non ?

# Econometrics vs. Machine Learning with Temporal Patterns

A few months ago, I did publish a (long) post entitled ‘some thoughts on economics, mathematics, econometrics, machine learning, etc‘. In that post, I was discussing possible differences between foundations of econometrics, and machine learning. I wanted to get back today on an important point, related to training/sampling datasets, when we have temporal data.

I was discussing this morning, with a student of the Data Science for Actuaries program, an interesting point related to claim frequency models, for insurance ratemaking. Since the goal is to predict claims frequency (to assess the level of the insurance premium), he suggested to use old data to train the model, and more recent one to test it. The problem is that the model did not incorporate any temporal pattern, and we got surprising results.

Consider here a simple dataset,

> set.seed(1)
> n=50000
> X1=runif(n)
> T=sample(2000:2015,size=n,replace=TRUE)
> L=exp(-3+X1-(T-2000)/20)
> E=rbeta(n,5,1)
> Y=rpois(n,L*E)
> B=data.frame(Y,X1,L,T,E)

Claims frequency is driven by a Poisson process, with one covariate, X1, and we assume that the intensity decreases (with an exponential rate). Consider here a standard linear regression, without any time effect

> reg=glm(Y~X1+offset(log(E)),data=B,
+ family=poisson)

We can also compute the empirical annualized claims frequency

> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B[abs(B$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp=Vectorize(p)(seq(.05,.95,by=.1)) and plot the two curves on the same graph, > plot(seq(.05,.95,by=.1),vp,type="b") > lines(u,exp(v),lty=2,col="red") This is what we usually do in econometrics. In machine learning, and more specifically to assess the quality of the model, and for model selection, it is common to split the dataset in two parts. A training sample, and a validation sample. Consider some randomized training/validation samples, then fit a model on the training sample, and finally use it to get a prediction, > idx=sample(1:nrow(B),size=nrow(B)*7/8) > B_a=B[idx,] > B_t=B[-idx,] > reg=glm(Y~X1+offset(log(E)),data=B_a, + family=poisson) > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B_a[abs(B_a$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_t=Vectorize(p)(seq(.05,.95,by=.1)) > lines(seq(.05,.95,by=.1),vp_t,col="red") The blue curve is the prediction on the training sample (as we usually do in econometrics), but then the red curve is the prediction on the testing sample. Here, volatility probably comes from the small size of the testing sample (1 observation out of 8). Now, what if we use the year as a splitting criteria : we fit a model on old years to fit a model, and we test it on recent years, > B_a=subset(B,T<2014) > B_t=subset(B,T>=2014) > reg=glm(Y~X1+offset(log(E)),data=B_a,family=poisson) > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B_a[abs(B_a$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_t=Vectorize(p)(seq(.05,.95,by=.1)) > lines(seq(.05,.95,by=.1),vp_t,col="red") Clearly, we miss something here… We were looking at such a graph this morning, and it took me some time to understand how training and validation samples were designed, and that there was a possible temporal effect (actually, this morning, it was based on a 3 year training sample, and a 1 year validation sample). Since there is a temporal pattern, let us capture it. As an econometrician, let me use a regression model > reg=glm(Y~X1+T+offset(log(E)),data=B, + family=poisson) > C=coefficients(reg) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C[1]+C[3]*(2000:2015))) > lines(u,v,lty=2,col="red") (I focus only on the evolution of the temporal variate on that graph). Here, we use a linear model, but there are usually no reason to assume linearity. So we might consider splines > library(splines) > reg=glm(Y~X1+bs(T)+offset(log(E)), + data=B,family=poisson) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > v2=predict(reg,newdata=data.frame(X1=0, + T=2000:2015,E=1)) > plot(2000:2015,exp(v2),type="b") > lines(u,v,lty=2,col="red") But here again, why should we assume that there is an underlying smooth function? There might be some ruptures… So let us consider a regression on factors > reg=glm(Y~0+X1+as.factor(T)+offset(log(E)), + data=B,family=poisson) > C=coefficients(reg) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C[2:17]),type="b") > lines(u,v,lty=2,col="red") An alternative might be to consider some more general model, like a regression tree > library(rpart) > reg=rpart(Y~X1+T+offset(log(E)),data=B, + method="poisson",cp=1e-4) > p=function(t){ + B=B[B$T==t,]
+   B$E=1 + mean(predict(reg,newdata=B)) + } > y_m=Vectorize(function(t) p(t))(2000:2015) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3+.5) > plot(2000:2015,y_m,ylim=c(.02,.085),type="b") > lines(u,v,lty=2,col="red") Here, it seems that something went wrong. I guess it’s coming from the exposure. So consider a simplier model, on the annualized frequency, and with weights that are related to the exposure > reg=rpart(Y/E~X1+T,data=B,weights=B$E,cp=1e-4)
> p=function(t){
+   B=B[B$T==t,] + B$E=1
+   mean(predict(reg,newdata=B))
+ }
> y_m=Vectorize(function(t) p(t))(2000:2015)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3+.5)
> plot(2000:2015,y_m,ylim=c(.02,.085),type="b")
> lines(u,v,lty=2,col="red")

That was for the econometrician perspective. With a machine learning perspective, consider a training sample (here based on old data) and a validation sample (based on more recent ones)

> B_a=subset(B,T<2014)
> B_t=subset(B,T>=2014)

If we consider a model, it is easy to get a prediction on recent years, even if the model was designed to model older ones,

> reg_a=glm(Y~X1+T+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[1]+C[3]*c(2000:2013,
+ NA,NA)),type="b")
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(C[1]+C[3]*2014:2015),
+ pch=19,col="blue")

But if we use years as factors, things are more complicated.

> reg_a=glm(Y~0+X1+as.factor(T)+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> RMSE=function(A){
+   L=exp(C[1]*B_t$X1+ A[1]*(B_t$T==2014) + A[2]*(B_t$T==2015)) + Y_t=L*B_t$E
+   sum( (Y_t - B_t$Y )^2)} > i=optim(c(.4,.4),RMSE)$par
> plot(2000:2015,c(exp(C[2:15]),NA,NA),)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(i),pch=19,col="blue")

becase we need to get a prediction on levels that were not in our training sample. Here, we minimize the RMSE to quantify factor levels for recent years. And the output is not that bad.

So yes, it is possible to get a training dataset on older data, and test it on recent years. But one should be careful, and take into account, properly, temporal patterns.

# Where People Live

There was an interesting map on reddit this morning, with a visualisation of latitude and longituge of where people live, on Earth. So I tried to reproduce it. To compute the density, I used a kernel based approch

> library(maps)
> data("world.cities")
> X=world.cities[,c("lat","pop")]
> liss=function(x,h){
+   w=dnorm(x-X[,"lat"],0,h)
+   sum(X[,"pop"]*w)
+ }
> vx=seq(-80,80)
> vy=Vectorize(function(x) liss(x,1))(vx)
> vy=vy/max(vy)
> plot(world.cities$lon,world.cities$lat,)
> for(i in 1:length(vx))
+ abline(h=vx[i],col=rgb(1,0,0,vy[i]),lwd=2.7)

For the other axis, we use a miror technique, to ensure that -180 is close the +180

> Y=world.cities[,c("long","pop")]
> Ya=Y; Ya[,1]=Y[,1]-360
> Yb=Y; Yb[,1]=Y[,1]+360
> Y=rbind(Y,Ya,Yb)
> liss=function(y,h){
+   w=dnorm(y-Y[,"long"],0,h)
+   sum(Y[,"pop"]*w)
+ }
> vx=seq(-180,180)
> vy=Vectorize(function(x) liss(x,1))(vx)
> vy=vy/max(vy)
> plot(world.cities$lon,world.cities$lat,pch=19)
> for(i in 1:length(vx))
+ abline(v=vx[i],col=rgb(1,0,0,vy[i]),lwd=2.7)

Now we can add the two, on the same graph

# Macro vs. Micro Methods in Non-Life Claims Reserving (an Econometric Perspective)

With Mathieu Pigeon, we finalized a short paper on  ‘macro’ vs. ‘micro’ methods in claims reserving. The paper is now online on arxiv

Traditionally, actuaries have used run-off triangles to estimate reserve (“macro” models, on agregated data). But it is possible to model payments related to individual claims. If those models provide similar estimations, we investigate uncertainty related to reserves, with “macro” and “micro” models. We study theoretical properties of econometric models (Gaussian, Poisson and quasi-Poisson) on individual data, and clustered data. Finally, application on claims reserving are considered.