# All posts by Arthur Charpentier

Arthur Charpentier, professor in Montréal, in Actuarial Science. Former professor-assistant at ENSAE Paristech, associate professor at Ecole Polytechnique and assistant professor in Economics at Université de Rennes 1.  Graduated from ENSAE, Master in Mathematical Economics (Paris Dauphine), PhD in Mathematics (KU Leuven), and Fellow of the French Institute of Actuaries.

# (Advanced) R Crash Course, for Actuaries

The fourth year of the Data Science for Actuaries program started this morning. I will be there for the introduction to R. The slides are available online (created with slidify, the .Rmd file is also available)

A (standard) markdown is also available (as well as the .Rmd file). I have to thank Ewen for his help on slidify (especially for the online quizz, and the integration of leaflet maps or the rgl animated graph….)

# Visualizing effects of a categorical explanatory variable in a regression

Recently, I’ve been working on two problems that might be related to semiotic issues in predictive modeling (i.e. instead of a standard regression table, how can we plot coefficient values in a regression model). To be more specific, I have a variable of interest $Y$ that is observed for several individuals $i$, with explanatory variables $\mathbf{x}_i$, year $t$, in a specific region $z_i\in\{A,B,C,D,E\}$. Suppose that we have a simple (standard) linear model (forget about time here) $$y_i=\beta_0+\beta_1x_{1,i}+\cdots+\beta_kx_{k,i}+\sum_j \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i$$

Let us forget the temporal effect to focus on the spatial effect today. And consider some simulated dataset. There will be only one (continuous) explanatory variable. And I will generate correlated covariates, just to be more realistic.

n=1000 library(mnormt) r=.5 Sigma &lt;- matrix(c(1,r,r,1), 2, 2) set.seed(1) X=rmnorm(n,c(0,0),Sigma) X1=cut(X[,1],c(-100,quantile(X[,1],c(.1,.4,.7,.85)), 100),labels=LETTERS[1:5]) X2=X[,2] Y=5+X[,1]-X[,2]+rnorm(n)/2 db=data.frame(Y,X1,X2)

Here we have $$y_i=\beta_0+\beta_1x_{1,i}+\sum_{j\in\{A,B,C,D,E\}} \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i$$ The goal here is to get to graph to visualize the vector $\hat\alpha=(\hat\alpha_A,\cdots,\hat\alpha_E)$. Let us run the linear regression

reg1=lm(Y~X1+X2,data=db) idx=which(substr(names(reg1$coefficients), 1,2)=="X1") v1=reg1$coefficients[idx] names(v1)=LETTERS[2:5] barplot(v1,col=rgb(0,0,1,.4))

Note that it is possible to add some sort of “confidence interval” to discuss significance (or to avoid to spend hours discussing differences in bar heights that are not significantly different)

library(Hmisc) sv1=summary(reg1)$coefficients[idx,2] (bp1=barplot(v1,ylim=range(c(0,v1+2*sv1)))) errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE) My main concern here is the “reference” that is considered. Should $A$ be the reference? Why not $B$ db$X1=relevel(db$X1,"B") reg1=lm(Y~X1+X2,data=db) idx=which(substr(names(reg1$coefficients),1,2)=="X1") v1=reg1$coefficients[idx] names(v1)=LETTERS[c(1,3:5)] library(Hmisc) sv1=summary(reg1)$coefficients[idx,2] (bp1=barplot(v1) errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE)

Why not the smallest one? Why not the largest one?… What if there is no simple way to choose. Furthermore, let us get back to the original point, which is that there might be some temporal aspects. More precisely, we can have $\hat\alpha^{(t)}=(\hat\alpha_A^{(t)},\cdots,\hat\alpha_E^{(t)})$. If we have also $\hat\alpha^{(t+1)}$ and we get another plot, how do we interpret it. If for $E$ the bar is taller, it means that relative to $A$, the difference has increased. I have the feeling that the interpretation is more complicated because we do not see, on that graph, changes in $\hat\alpha^{(t)}_A$.

Let us try something else. First, let us get back to the original setting

db$X1=relevel(db$X1,"A")

Consider here the regression without the intercept, so that all values remain

reg1=lm(Y~0+X1+X2,data=db) idx=which(substr(names(reg1$coefficients),1,2)=="X1") v1=reg2$coefficients[idx] names(v1)=LETTERS[1:5] barplot(v1)

It can be hard to read, especially if $Y$ takes (very) large values, and you think that barplots should start at 0. But still, having those 5 values is nice. Why not rescale that graph?

A natural idea my be to consider the case where no spatial component is considered, and to look at the difference with that reference.

 reg1=lm(Y~1+X2,data=db) reg2=lm(Y~0+X1+X2,data=db) idx=which(substr(names(reg2$coefficients),1,2)=="X1") v1=reg2$coefficients[idx] v2=v1-reg1$coefficients["(Intercept)"] barplot(v2,col=rgb(0,0,1,.4)) sv2=summary(reg2)$coefficients[idx,2] (bp2=barplot(v2,ylim=range(c(v2-2*sv2,v2+2*sv2)))) errbar(bp2[,1],v2,v2-2*sv2,v2+2*sv2,add=TRUE)

I like that graph, I should admit it. Now, I still have some remaining questions. For instance, can we insure that when only the intercept is considered, the value of $\hat\beta_0$ is somewhere between $\hat\beta_A,\cdots,\hat\beta_E$? Is it possible that $\hat\beta_A-\hat\beta_0,\cdots,\hat\beta_E-\hat\beta_0$ are all positive? In that case, I would find that hard to interpret.

Actually, if I really want values that can be seen as compared to some average, why not consider a (weighted) average of $\hat\beta_A,\cdots,\hat\beta_E$? (weights being here proportion in each class, in each region)

 w=table(db$X1) v3=v1-sum(w*v1)/sum(w) (bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3)))) errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE) I like that one. But what if, instead of normalizing at the end, we normalize the original dependent variable. By “normalize”, I mean “rescale”, to have a centered variable. db$Y0=db$Y-mean(db$Y) reg3=lm(Y0~0+X1+X2,data=db) sv3=summary(reg3)$coefficients[idx,2] (bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3)))) errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE) This one is nice, because it is extremely simple to explain. But what if instead of a linear regression, we add a logistic one (with $Y\in\{0,1\}$)? or a Poisson regression… So maybe it cannot be the best solution here. Let us try something else… In insurance ratemaking, people like to use “zonier“. It is a two-stage regression. The idea is to run a regression without any spatial components, first. Then, consider the regression of residuals on spatial variables. Here, it would be something like reg1=lm(Y~1+X2,data=db) reg2=lm(Y~0+X1+X2,data=db) Since we focus on residuals, those are centered, and we have an easy interpretation of respective values  sv4=summary(reg4)$coefficients[idx,2] v4=reg4$coefficients (bp4=barplot(v4,names.arg=LETTERS[1:5]))) errbar(bp4[,1],v4,v4-2*sv4,v4+2*sv4,add=TRUE) I guess that it can also be use in generalized linear models, with Pearson (or deviance) residuals. Another possible idea can be the following. Again, the goal is not to have the true values, but to visualize on a graph how regions can be different. Here, all of them are significantly different. And in region $A$, $Y$ is smaller, ceteris paribus (other things equal in the sense that we have taken into account $x_1$). And in region $E$ it is larger. Here, the graph helps to “see” those differences. Why not consider a completely different graph. What if we plot vector $a$ instead of $\alpha$, where $a_A$ can be interpreted as the value of the coefficient if we consider region $A$ against “not region $A$“. What if we consider 5 regressions where dichotomous versions of $Z$ are considered : $Z_j=\mathbf{1}_{Z=j}$. v5=sv5=rep(NA,5) names(v5)=LETTERS[1:5] for(k in 1:5){ reg=lm(Y~I(X1==LETTERS[k])+X2,data=db) v5[k]=reg$coefficients[2] sv5[k]=summary(reg)$coefficients[2,2]} We can plot that sequence of values, including some confidence intervals (that would be related to significance with respect to all other regions) (bp5=barplot(v5,ylim=range(c(v5-2*sv5,v5+2*sv5)))) errbar(bp5[,1],v5,v5-2*sv5,v5+2*sv5,add=TRUE) Looking at values does not give intuitive results, but I have the feeling that it is easy to explain what we plot (we compare each region to “the rest of the world”), and the ordering of $a$ seems to be consistent with $\alpha$ (but I could not prove it). Here are some ideas I got. I should be able to provide other graphs, but I would love to discuss with anyone on that topics, to find a proper and nice way to visualize effects of a categorical explanatory variable in a regression model (that can be a logistic one). Comments are open… # Holt-Winters with a Quantile Loss Function Exponential Smoothing is an old technique, but it can perform extremely well on real time series, as discussed in Hyndman, Koehler, Ord & Snyder (2008)), when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class. Furthermore, I like it because I think it has nice pedagogical features. Consider simple exponential smoothing, $$L_{t}=\alpha Y_{t}+(1-\alpha)L_{t-1}$$ where $\alpha\in(0,1)$ is the smoothing weight. It is locally constant, in the sense that ${}_{t}\hat Y_{t+h} = L_{t}$  library(datasets) X=as.numeric(Nile) SimpleSmooth = function(a){ T=length(X) L=rep(NA,T) L[1]=X[1] for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]} return(L) } plot(X,type="b",cex=.6) lines(SimpleSmooth(.2),col="red") When using the standard R function, we get hw=HoltWinters(X,beta=FALSE,gamma=FALSE, l.start=X[1]) hw$alpha [1] 0.2465579

Of course, one can replicate that optimal value

V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(erreur^2)) } optim(.5,V)$par [1] 0.2464844 Here, the optimal value for $\alpha$ is the one that minimizes the one-step prediction, for the $\ell_2$ loss function, i.e. $$\sum_{t=2}^n(Y_t-{}_{t-1}\hat Y_t)^2$$ where here ${}_{t-1}\hat Y_t = L_{t-1}$. But one can consider another loss function, for instance the quantile loss function, $$\ell_{\tau}(\varepsilon)=\varepsilon(\tau-\mathbb{I}_{\varepsilon\leq 0})$$. The optimal coefficient is then obtained using HWtau=function(tau){ loss=function(e) e*(tau-(e&lt;=0)*1) V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(loss(erreur))) } optim(.5,V)$par }

Here is the evolution of $\alpha^\star_\tau$ as a function of $\tau$ (the level of the quantile considered).

T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(T,HW,type="l") abline(h= hw$alpha,lty=2,col="red") Note that the optimal $\alpha$ is decreasing with $\tau$. I wonder how general this result can be… Of course, one can consider more general exponential smoothing, for instance the double one, with $$L_t=\alpha Y_t+(1-\alpha)[L_{t-1}+B_{t-1}]$$and$$B_t=\beta[L_t-L_{t-1}]+(1-\beta)B_{t-1}$$so that the prediction is now ${}_{t}\hat Y_{t+h} = L_{t}+hB_t$ (it is now locally linear – and no longer constant). hw=HoltWinters(X,gamma=FALSE,l.start=X[1]) hw$alpha alpha 0.4200241 hw$beta beta 0.05973389 The code to compute the smoothed series is the following DoubleSmooth = function(a,b){ T=length(X) L=B=rep(NA,T) L[1]=X[1]; B[1]=0 for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] } return(L+B) } Here also it is possible to replicate R using the $\ell_2$ loss function V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(erreur^2)) } optim(c(.5,.05),V)$par [1] 0.41904510 0.05988304

(up to numerical optimization approximation, I guess). But here also, a quantile loss function can be considered

HWtau=function(tau){ loss=function(e) e*(tau-(e&lt;=0)*1) V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(loss(erreur))) } optim(c(.5,.05),V)$par } and we can plot those values on a graph T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(HW[1,],HW[2,],type="l") abline(v= hw$alpha,lwd=.4,lty=2,col="red") abline(h= hw$beta,lwd=.4,lty=2,col="red") points(hw$alpha,hw$beta,pch=19,col="red") (with $\alpha$ on the $x$-axis, and $\beta$ on the $y$-axis). So here, it is extremely simple to change the loss function, but so far, it should be done manually. Of course, one do it also for the seasonal exponential smoothing model. # Time for a break # Justice “actuarielle”, algorithmes… et données Il y a un peu plus d’un an, Virginie Gautron m’envoyait plein de documents sur la “justice actuarielle”, concept que je découvrais alors. Pour comprendre un peu de quoi il s’agit, je peux renvoyer vers pénalité et gestion des risques : vers une justice « actuarielle » en Europe ? qui dresse un état de l’art, en francais. J’avoue avoir un peu mis de côté ensuite (faute de temps), et puis l’autre jour, j’ai dévoré surveillé et punir de Nicolas Bourgoin, qui contient un chapitre sur le sujet. Sans vouloir tomber dans le cliché des deux cultures (à la CP Snow), il faut bien admettre que pour un économètre, lire des écrits de juristes parlait de modèles prédictifs (ou actuariels) est une expérience marquante. Et ce qui me trouble (j’ai d’ailleurs le même soucis avec la psychologie – c’est un sujet proche) c’est que l’accent est souvent mis sur le modèle, l’algorithme, ce truc qui semble fasciner les juristes, alors qu’il me semble que le coeur du problème, ce sont les données. On nous dit, par exemple, que les algorithmes sont racistes (c’est un peu la thèse de Cathy O’Neil dans Weapons of Math Destruction), mais le problème ce n’est pas le modèle. C’est le fait qu’il y a un énorme biais de sélection dans les données ! Et selon l’adage “garbage in, garbage out“… dit de manière plus policée, je pourrais tenter un “purgamentum init, exit purgamentum” pour parler comme un juriste. Ou dit autrement, si on entre des mauvaises données en entrée, faut pas espérer grand chose en sortie… Pour comprendre le biais de sélection, on peut relire Mostly Harmless Econometrics de Joshua Angrist et Jörn-Steffen Pischke. En particulier, dans l’introduction, ils posent la question simple suivant : vaut-il mieux, suite à une opération, renvoyer les gens chez eux, ou bien les laisser en chambre de convalescence, à l’hôpital ? Le soucis est que les données sont biaisées : oui, si on regarde les données brutes, ceux qui sont rentré chez eux se sont remis plus vite… mais peut être a-t-on autorisé les gens à rentrer justement parce qu’ils allaient bien ? Et on retrouve ce biais de sélection un peu partout. Probablement même encore davantage dans les problèmes relatifs à la justice. Prenons un exemple simple. On va s’intéresser à la probabilité d’être coupable d’un crime, et d’être condamné par la justice. C’est ma variable $Y$. Pour avoir un modèle économétrique, il me faut des variables explicatives. J’ai une variable $X_3$ et une autre $X_4$, la première influençant vraiment la probabilité de commettre un crime, et l’autre étant juste du bruit. Maintenant, entrons dans le coeur du modèle. Supposons que j’ai une variable corrélée à ma richesse $\tilde{X}_2$, comme le lieu d’habitation, un critère racial, ou pour faire simple le diplôme. Notons $X_1$ le niveau de diplôme de la personne (mais encore une fois, on peut imaginer n’importe quoi, l’idée est que cette variable est corrélée à la richesse). Et surtout, cette variable n’est pas du tout une variable causale dans mon modèle. La probabilité de commettre un crime n’est pas influencée par $X_1$. Une personne ayant un diplôme de second cycle universitaire a autant de chance de commettre un crime qu’une personne ayant raté le brevet. On va même supposer que la probabilité de commettre un crime soit indépendante de la richesse. En fait, la richesse me permet de payer un bon avocat. Et avoir un bon avocat permet de diminuer la probabilité d’être condamné. Pour faire simple, j’ai deux niveaux de richesse, qui me donne soit un avocat commis d’office, soit un avocat qui s’impliquera peut être davantage (je simplifie à des fins pédagogiques). Bref, mon modèle aura la forme suivante : la probabilité d’être condamné pour un crime s’écrit $$p=(0.025+0.05\cdot X_3)\cdot[k\cdot \mathbb{1}(X_2=R)+\mathbb{1}(X_2=P)]$$ Autrement dit, ma probabilité est une fonction croissante (linéaire) de $X_3$, et si je suis riche (ce qu’on note $X_2=R$), j’ai moins de chance d’être condamné, avec un facteur $k\in[0,1]$, toutes choses étant égales par ailleurs, ceteris paribus, histoire de parler latin comme les juristes. Pour simuler un jeu de données, on considère une population d’un million de personnes, et je suppose que ma richesse est très corrélée avec mon niveau de diplôme. Pour info, la classe $A$ ce sont les gens les plus diplômés, et je suppose qu’un bon avocat divise par trois la chance d’être condamné.  n=1000000 r=.95 k=1/3 S=matrix(c(1,r,r,1),2,2) library(mnormt) set.seed(1) vectX=rmnorm(n,varcov=S) X1=cut(vectX[,1],breaks=c(-10,1,-1,10), c("C","B","A")) X1=relevel(X1,"A") X2=cut(vectX[,2],breaks=c(-10,1,10),c("P","R")) X3=runif(n) X4=runif(n) Y1=rbinom(n,size = 1,prob=(.025+.05*X3)*(k*(X2=="R")+(X2=="P"))) B0=data.frame(X1,X2,X3,X4,Y1) On peut faire un peu de statistiques pour regarder les caractéristiques des gens condamnés (uniquement) B1=data.frame(X1=X1[Y1==1],X2=X2[Y1==1],X3=X3[Y1==1],X4=X4[Y1==1]) Par exemple si je regarde la proportion de gens les plus diplômés: ils sont 15% dans la population totale, et 5% parmi les personnes condamnées  mean(B0$X1=="A") [1] 0.158734 mean(B1$X1=="A") [1] 0.05016072 t.test(B0$X1=="A",B1$X1=="A") data: B0$X1 == "A" and B1$X1 == "A" t = 219.5, df = 1283500, p-value &lt; 2.2e-16 95 percent confidence interval: 0.1076038 0.1095428 sample estimates: mean of x mean of y 0.15873400 0.05016072 Autrement dit, on a une différence significative entre ces deux pourcentages. Pareil si je regarde la valeur moyenne de $X_3$, qui était de 50% dans la population et de 58% dans la population condamnée  t.test(B0$X3,B1$X3) data: B0$X3 and B1$X3 t = -163.46, df = 844940, p-value &lt; 2.2e-16 95 percent confidence interval: -0.08455430 -0.08255058 sample estimates: mean of x mean of y 0.5000094 0.5835618 alors que pour $X_4$, on a 50% dans la population totale, et 50% dans la population condamnée. Mais on s’y attendait….  t.test(B0$X4,B1$X4) data: B0$X4 and B1$X4 t = -0.81161, df = 48645, p-value = 0.417 95 percent confidence interval: -0.003869879 0.001603452 sample estimates: mean of x mean of y 0.5001326 0.5012659 Passons à la seconde phase. Il semble y avoir un vrai intérêt dans la littérature sur les modèles actuariels pour la construction de modèles liés à la probabilité de récidive. Là encore, je vais faire une hypothèse forte : faire de la prison n’influence en rien la probabilité de récidiver. Autrement dit, pour cette population déjà condamnée, je vais refaire tourner mes simulations, avec une probabilité de commettre un crime qui dépend seulement de $X_3$, mais avec toujours notre effet richesse qui permet de payer un avocat qui fera baisser la probabilité d’être condamné.  set.seed(2) n1=nrow(B1) Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P"))) B2=cbind(B1,Y2) reg=glm(Y2~X1+X3+X4,data=B2,family=binomial)

J’ai construit ici mon modèle permettant de calculer une probabilité de récidive. Le modèle est construit sur le niveau de diplôme $X_1$ (qui n’influence en rien, encore une fois, la probabilité de commettre un crime), la variable $X_3$ (qui effectivement joue un rôle) et la variable de bruit $X_4$.

 set.seed(2) n1=nrow(B1) Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P"))) B2=cbind(B1,Y2) reg=glm(Y2~X1+X3+X4,data=B2,family=binomial) Je fais ici un score logistique car c’est le plus logique.  summary(reg) Coefficients: Estimate Std. Error z value Pr(|z|) Inter.) -4.04753 0.11592 -34.916 2e-16 *** X1C 0.61311 0.10914 5.618 1.93e-08 *** X1B 0.68807 0.09977 6.896 5.34e-12 *** X3 1.00610 0.08052 12.495 2e-16 *** X4 -0.13562 0.07286 -1.861 0.0627 . Premier élément à observer, la variable $X_4$ est à la limite d’être significative… Certains auteurs sur ces thèmes diraient que c’est un signal faible. Pourtant, en l’occurence, c’est du bruit ! Juste du bruit !! Et parmi les signaux forts, très clairement significatifs, on retrouve $X_1$, notre niveau de diplôme. Ici, si on est faiblement qualifié (niveaux $B$ ou $C$) alors on a une chance de récidive beaucoup plus importante ! C’est très fortement significatif ! Alors que pour rappel, la probabilé de commettre un crime, et de récidiver n’est aucunement causé par cette variable ! C’est juste une variable corrélée avec la variable richesse, qui permet de faire baisser la probabilité d’être condamné. Plus intéressant, si je regarde la valeur du coefficient, en fonction de $k$, autrement dit la valeur d’un bon avocat (que l’on peut s’offrir si on est riche) on obtient la courbe suivante. Autrement, sur la partie droite, même si un bon avocat ne baisse que de 5% la chance d’être condamné, le niveau de diplôme semble être significativement liée à la probabilté de récidive. Avec un pic autour de 30%. Si les avocats ont un pouvoir incroyable, par exemple en divisant par 20 la probabilité d’être condamné, étrangement, la variable de diplôme redevient très peu significative, mais je pense que c’est parce qu’on n’a alors plus personne de très diplômé dans la base. On pourrait me reprocher d’avoir une corrélation trop forte entre richesse et niveau de diplôme, ici. En mettant une corrélation à 70% (dans mon modèle Gaussien latent), j’obtiens qui donne un effet moindre, mais là encore, si un bon avocat permet de réduire la probabilité d’être condamné par un facteur compris entre 5% et 90%, alors la variable de diplôme semblera significative. Mon modèle n’est pas faux ici. C’est essentiellement que j’ai un biais dans mes données. Et ne pas en tenir compte fait que le modèle n’a aucun sens, et ne peut pas être interprété sans commettre d’erreurs (et ici je cumule les âneries, en voyant un signal faible alors qu’il y a du bruit, et en voyant un signal fort pour une variable qui n’a aucun lien causal). J’insiste lourdement, mais le diplôme n’intervient aucunement ici. Et on peut remplacer le diplôme par toute variable corrélée à la richesse, comme locataire et propriétaire, etc. Bref, l’utilisation de modèles économétrique dans des problèmes judiciaires n’est pas simple ! Je parle ici de biais de sélection, mais la formulation qui me semble la plus juste serait une mauvaise spécification du modèle (on n’utilise pas les bonnes variables, mais ces dernières sont peut-être non-observables – ou difficiles à justifier éthiquement : peut-on demander le revenu de la personne pour faire un score de récidive, et avoir un score qui dépende effectivement du revenu ?). Il y a aussi un problème d’identification, car on a du mal à savoir ce qui influence le fait de commettre un crime, et la condamnation. Ce qu’on observe, ce sont les condamnations, pas les actes (je ne parle pas des condamnations à tort ici, que j’avais pu évoquer dans un autre billet, sur le principe de précaution). Mais j’essayerais de revenir sur ce point après les vacances…. Dernière anecdote pour finir. Après avoir découvert ces modèles via Virginie Gautron, j’avais discuté avec une amie, présidente de tribunal, qui me parlait de la nullité des experts qu’elle avait sous la main. Et elle voulait faire un peu de ménage. Bref, elle me demandait si je pouvais faire un modèle attribuant un score à un expert. Et là encore, ma réponse a été que oui, faire un modèle c’est facile, à condition d’avoir des données. Comment puis-je avoir un score $Y$ d’un expert? Qui sait ce qu’est un bon expert? C’est quelqu’un contre qui on ne demande pas de contre-expertise? Malheureusement, la justice ne fonctionne pas ainsi, et la demande n’a rien à voir avec la qualité de l’expert. Et la qualité, vu du point de vue de la juge, est une notion qui peut être très différente de celle des deux parties qui s’opposent… Bref, je pourrais facilement faire un modèle. Mais en avoir un qui ait du sens, c’est une autre paire de manches ! # Actuariat de l’Assurance Non-Vie #11 Cette semaine, nous allons finir la modélisation du passif d’une compagnie d’assurance non-vie. Les slides sont en ligne. # Actuariat de l’Assurance Non-Vie #10 La semaine dernière, nous avions commencé la modélisation du passif d’une compagnie d’assurance non-vie. Les slides introductifs sont en ligne. # Optimal Portfolios Exam Pour l’examen du cours de Gestion d’Actifs (optimal portfolios) un cas pratique avec des vraies données, disponibles sur loc="http://freakonometrics.free.fr/base_portefeuille.csv" base=read.csv2(loc,sep=";") base[,1]=as.Date(as.character(base[,1]),"%Y-%m-%d") for(i in 2:10) base[,i]=as.numeric(as.character(base[,i])) Pour le rendu, je demande un markdown dont la trame est en ligne ici, exam_portfolio.Rmd, et la sortie est exam_portfolio.html. Il faut télécharger le fichier Rmd, l’ouvrir avec RStudio, puis me renvoyer le fichier rempli (ainsi que le fichier html compilé) pour midi, à arthur.charpentier@univ-rennes1.fr. # Actuariat de l’Assurance Non-Vie #9 Cette semaine, nous avons fini les modèles de tarification, avec une extention sur les modèles globaux (sans passer par une séparation entre fréquence et coût moyen), et sur le choix de variables, et le choix de modèles. Les slides introductifs sont en ligne. # The myth of interpretability of econometric models There are important discussions nowadays about data modeling, to choose between the “two cultures” (as mentioned in Breiman (2001)), i.e. either econometrics models or machine/statistical learning models. We did discuss this issue recently in Econométrie et Machine Learning (so far only in French) with Emmanuel Flachaire and Antoine Ly. One argument often used by econometricians is the interpretability of econometric models. Or at least the attempt to get an interpretable model. We also have this discussion in actuarial science, for instance in ratemaking (or insurance pricing). Machine learning based models usually perform better (for some a priori chosen metric), but actuaries claim that econometric models are more easily interpretable. In actuarial literature, we assume that claim frequency $Y$ is driven by some non-observable risk factor $\Theta$, and therefore, we do have heterogeneous risks in our portfolio. And, it can be seen as legitimate to differentiate prices. Assume that this risk factor $\Theta$ is strongly correlated with $X_1$, the age of the driver. Because in our portfolio, old drivers tend to have more accidents. Here, we could pretend to have a “causal story” (as defined in Freedman (2009)) because of a possible interpretation of the model. So it is natural here to consider a regression model of $Y$ on $X_1$ to derive our actuarial pricing model. But assume that, possibly, risk factor $\Theta$ is also strongly correlated with $X_2$, that can be related to spatial features (say latitude, which denoted a north/south position). Because in our portfolio, drivers living in the south tend to have more accidents (reads are known to be more dangerous there). Here, we could pretend to have a second “causal story”. Of course, since $\Theta$ is strongly correlated with $X_1$ and $X_2$, it means that $X_1$ and $X_2$ are strongly correlated. Here also, this correlation can be interpreted (not in a causal way as previously, but still), since we know that old people like to live in southern regions. So, what should we do here ? Let us run some simulations to illustrate.  set.seed(123) n=1e5 Theta=rnorm(n) X1=Theta+rnorm(n)/8 X2=Theta+rnorm(n)/8 L=exp(-3+Theta) Y=rpois(n,L) B=data.frame(Y,X1,X2) Our first idea was to consider a model where $Y$ is “explained” by the first variable $X_1$,  g1=glm(Y~X1,data=B,family=poisson) summary(g1) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Inter.) -2.97778 0.01544 -192.88 &lt;2e-16 *** X1 0.97926 0.01092 89.64 &lt;2e-16 *** As expected, our variable is “significant”, but also, probably more interesting, $X_2$, has no impact on the residuals  B$e1=residuals(g1,type="pearson") g1e=lm(e1~X2,data=B) summary(g1e)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Inter.) 0.0003618 0.0031696 0.114 0.909 X2 0.0028601 0.0031467 0.909 0.363

The interpretation is that once we corrected claim frequency for the age of the drivers, there is no spatial effect here. So, a good model should be based only on the age of the drivers.

But we can also consider the other story. We can consider a model where $Y$ is “explained” by the second variable $X_2$,

 g2=glm(Y~X2,data=B,family=poisson) summary(g2)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Inter.) -2.97724 0.01544 -192.81 &lt;2e-16 *** X2 0.97915 0.01093 89.56 &lt;2e-16 ***

Here also we have a valid model, that can be interpreted, and here also $X_1$, has no impact on the residuals

 B$e2=residuals(g2,type="pearson") g2e=lm(e2~X1,data=B) summary(g2e) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Inter.) 0.0004863 0.0031733 0.153 0.878 X1 0.0027979 0.0031504 0.888 0.374 The story is similar here. If we correct from the spatial pattern, claims frequency does not depend on the age of the driver. So, what should we do now? We do have two models, and each of them is as interpretable as the other one. Note that we can not use any statistical tool to distinguish the two: they are comparable  AIC(g1) [1] 51013.39 AIC(g2) [1] 51013.15 Why not incorporate the two explanatory variables $X_1$ and $X_2$, at the same time, in our regression model, and let “the model” decide what to do…?  g=glm(Y~X1+X2,data=B,family=poisson) summary(g) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Inter.) -2.98132 0.01547 -192.723 2e-16 *** X1 0.49310 0.06226 7.920 2.38e-15 *** X2 0.49375 0.06225 7.931 2.17e-15 *** It looks like we completely lost the interpretability of the model, since our two explanatory variables are (strongly) correlated. Actually, instead of saying “use one, and drop the other one (since it brings no further information)”, it says “use both, each one will explain half of the variable”. Strange interpretation, isn’t it? So why not try some LASSO here? library(glmnet) fit=glmnet(x=as.matrix(B[,c("X1","X2")]), y=B$Y,family="poisson") plot(fit,xvar="lambda")

Here also, it says that we either keep both, or none. So it cannot be used for variable selection (which is an important motivation to use LASSO technique). So, what should be do if we several interpretable models, but no way to choose? Because usually, we claim that we prefer to use a model with an interpretation. But what should be done here?

# Optimal Portfolios #2

Next week, we will continue the crash course on financial portfolio optimization, with application in R. The slides are available here, and R codes from there (in a Markdown).

The first part is still online, here.

# Faire (rapidement) un zonier

Dans le cours d’actuariat de l’assurance non-vie, on avait évoqué rapidement l’idée de faire un zonier. Autrement dit, on souhaite créer une variable polytomique (avec 4 ou 5 classes) représentant le critère spatial du risque. On va tenter de segmenter le territoire, en regroupant l’ensemble des territoires ‘proches’ (en terme de risque) étudiées en un nombre de classes prédéfini, permettant de tenir compte d’un (potentiel) critère spatial dans la tarification. Pour les aspects pratiques,  comme on va manipuler des données spatiales, chargeons quelques packages pour commencer

library(maptools)
library(rgeos)
library(rgdal)
library(ggplot2)
library(plyr)
library(maptools)
library(cartography)

On va aussi avoir besoin de données spatialisées (disons en France métropolitaire). Pour faire simple, je vais tirer des codes insee de communes au hasard, un millier, et une variable $Y$ qui va nous intéresser. Et je suppose qu’elle dépend de $X$ que je n’observe pas, mais qui est liée à des caractéristiques spatiales (en gros la latitude, i.e. un positionnement nord/sud)

download.file("http://freakonometrics.free.fr/popfr19752010.csv","popfr.csv")
base$insee = base$dep*1000+base$com n=1000 set.seed(123) id=sample(1:nrow(base),size=n) simbase=data.frame(insee=base$insee[id])
X=(46-base$lat[id])/2 simbase$Y=X+rnorm(n)

Ici, on a la base suivante (mais on va supposer que je n’observe pas $X$ – qui n’a servi qu’à simuler des données)

> head(simbase)
insee          Y
1 61499 -1.8573363
2 19181 -0.7059191
3 55307 -0.5923649
4 74030  0.6098773
5 30328 -0.4584795
6 10050 -1.2361240

Classiquement, $Y$ peut être un résidu (normalisé) de régression, par exemple. Pour visualiser, il faut un fond de carte

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds")
FR0=readRDS("FRA_adm0.rds")

Ensuite, on met les points (en fusionnant nos données avec la base insee donnant latitude et longitude des codes insee de communes)

plot(FR0)
simbase = merge(simbase,base[,c("insee","long","lat")])
cols = rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10))
bk = seq(-5,4.5,length=21)
cuty = cut(simbase$Y,breaks=bk,labels=1:20) points(simbase$long,simbase$lat, col=cols[cuty],pch=19,cex=.5) On va voir eux techniques pour faire un zonier. Le premier est de travailler par zone prédéfinie, comme le département, simbase$dpt=trunc(simbase$insee/1000) A=aggregate(x = simbase$Y,by=list(simbase$dpt),mean) names(A)=c("dpt","y") download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm2.rds","FRA_adm2.rds") FR=readRDS("FRA_adm2.rds") donnees_carte=data.frame(FR@data) d=donnees_carte$CCA_2
d[d=="2A"]="201"
d[d=="2B"]="202"
donnees_carte$dpt=as.numeric(as.character(d)) donnees_carte=merge(donnees_carte,A,all.x=TRUE) donnees_carte=donnees_carte[order(donnees_carte$OBJECTID),]
bk=seq(-2.75,2.75,length=21)
donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:20)
cols = rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10))
plot(FR, col=cols[donnees_carte$cuty],xlim=c(-5.2,12)) autrement dit, on fait une carte choroplèthe (il faut s’assurer que les couleurs sont mises au bon endroit). Ensuite, on peut définir deux zones, bk=seq(-2.75,2.75,length=3) donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:2) plot(FR, col=cols[c(4,16)][donnees_carte$cuty],xlim=c(-5.2,12))

ou quatre

bk=seq(-2.75,2.75,length=5)
donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:4)
plot(FR, col=cols[c(3,8,12,17)][donnees_carte$cuty],xlim=c(-5.2,12)) On va créer uen variable à 4 modalités, à partir des départements. Une autre solution consiste à utiliser les données spatiales, obtenues par fusion avec la base insee simbase=merge(simbase,base[,c("insee","long","lat")]) On va ensuite se donner une grille, en France. C’est un peu pénible. Il nous faut le polynôme de contour de la France P1=FR0@polygons[[1]]@Polygons[[355]]@coords P2=FR0@polygons[[1]]@Polygons[[27]]@coords plot(FR0,border=NA) polygon(P1) polygon(P2) (ici on a juste la France métropolitaine, et la Corse – deux polygônes donc) et on part d’une grille sur un rectangle grille<-expand.grid(seq(min(simbase$long),max(simbase$long),length=101),seq(min(simbase$lat),max(simbase$lat),length=101)) paslong=(max(simbase$long)-min(simbase$long))/100 paslat=(max(simbase$lat)-min(simbase$lat))/100 On retient alors juste les points qui sont dans les polygônes f=function(i){ (point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P1[,1],P1[,2])>0)+(point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P2[,1],P2[,2])>0) } indic=unlist(lapply(1:nrow(grille),f)) grille=grille[which(indic==1),] points(grille[,1]+paslong/2,grille[,2]+paslat/2,cex=.4,pch=19,col="blue") voilà le résultat Maintenant, on peut utiliser du krigeage mais j’ai plutôt voulu tenter des plus proches voisins. Pour chaque point de la grille, on prend la moyenne des plus proches voisins (à vol d’oiseau, i.e. avec une norme Euclidienne – sur la sphère) library(geosphere) knn=function(i,k=20){ d=distHaversine(grille[i,1:2],simbase[,c("long","lat")], r=6378.137) r=rank(d) ind=which(r<=k) mean(simbase[ind,"Y"]) } grille$y=Vectorize(knn)(1:nrow(grille))
bk=seq(-2.75,2.75,length=21)
grille$cuty=cut(grille$y,breaks=bk,labels=1:20)
cols <- rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10))
points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[grille$cuty],pch=19) Ici, ca donne la carte suivante mais là encore, on peut retenir deux niveaux bk=seq(-2.75,2.75,length=3) grille$cuty=cut(grille$y,breaks=bk,labels=1:2) plot(FR0,border=NA) polygon(P1) polygon(P2) points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[c(4,16)][grille$cuty],pch=19)

ce qui donne les zones suivantes

bk=seq(-2.75,2.75,length=5)
grille$cuty=cut(grille$y,breaks=bk,labels=1:4)
plot(FR0,border=NA)
polygon(P1)
polygon(P2)
points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[c(3,8,12,17)][grille$cuty],pch=19) ou bien, si on retient quatre niveaux A partir de là, pour créer une variable de zone, le plus simple est à un point donné de chercher le point de la grille le plus proche, et de lui donner la valeur de la couleur associée pred=function(z){ d=distHaversine(z,grille[,1:2], r=6378.137) grille[which.min(d),"cuty"] } Facile, non ? # Données sinistres, devoir Actuariat (ENSAE) Pour le dernier devoir du cours d’actuariat non-vie, à l’ENSAE, quelques données de sinistres. Pour plusieurs milliers de sinistres (garantie dommages tous accidents issu d’un portefeuille automobile d’une société française) on a les données suivantes, entre 1995 et 2014, • id_sin: identifiant du sinistres (les 4 premiers caractères représentent l’exercice de survenance) • an_surv: exercice de survenance • an_gest: exercice de gestion • etat_sin: code état du sinistre à la fin de l’exercice de gestion (1=mis en jeu, 2= cloture partielle, 3 = cloture totale, 4 = réouvert et 5 = cloture sans suite) • pmt: paiements effectués sur le sinistre à la fin de l’exercice de gestion • rec: recours encaissés sur le sinistre à la fin de l’exercice de gestion • eval_pmt: évaluation des paiements (réglés+à venir) vue à la fin de l’exercice de gestion • eval_rec: évaluation des recours (réglés+à venir) vue à la fin de l’exercice de gestion Pour information complémentaire, on dispose des éléments d’information suivants, sur le portefeuille  Exercice Volume portefeuille (garanties années) Chiffre d’affaires (€) Nb sinistres survenus dans l’année 2014 303 900 45 900 000 41 678 2013 296 500 45 700 000 41 639 2012 295 600 44 000 000 42 815 2011 292 000 41 200 000 43 351 2010 282 900 40 700 000 45 276 2009 274 300 39 800 000 42 779 2008 263 800 41 900 000 38 941 2007 256 200 41 200 000 38 048 2006 250 800 42 400 000 38 603 2005 239 100 41 100 000 37 539 2004 226 800 37 900 000 34 728 2003 220 300 39 200 000 34 491 2002 211 400 35 500 000 35 984 2001 202 400 34 800 000 36 322 2000 193 400 34 000 000 36 186 1999 184 500 32 300 000 34 123 1998 175 500 29 800 000 30 517 1997 166 500 28 000 000 28 432 1996 157 500 26 500 000 26 852 1995 148 600 25 700 000 26 775 La base complète est dans le fichier ensae2017.zip. Nous reviendrons sur les attendus lors du prochain cours. # Tous des (potentiels) terroristes ? En fin de semaine dernière, je commencais surveiller et prévenir: L’ère de la pénalité prédictive de Nicolas Bourgoin, et au début, il évoque des mesures de sécurité renforcées dans un rayon de 20km des gares, des aéroports et des ports, qui pourraient être prises dans le cadre de la réforme de lois sur la sécurité intérieure. Cette information était reprise récemment sur rtl, Le périmètre est aussi agrandi : les contrôles pourront avoir lieu aux abords des gares internationales (et non plus seulement à l’intérieur) ainsi que dans un rayon de 20 kilomètres autour des aéroports et des ports. ou sur le site de l’obs ces opérations de contrôles seront mises en place “aux abords” de 373 gares, ports et aéroports, dans un rayon de 20 kilomètres. Une extension considérable puisque jusqu’à présent, ces contrôles restaient cantonnées à l’intérieur de ces espaces accessibles au public. En lisant le livre, je trouvais incroyable cette histoire de 20km, car tout le monde habite à moins de 20km d’une gare (ne sachant trop quelles pouvaient être ces 373 gares, je pars du fait que toutes les gares pourraient être concernées). C’est un peu ce que note la cimade, en écrivant Le projet de loi prévoit ainsi de permettre les contrôles d’identité aux frontières pour une durée de 12 heures (contre 6 aujourd’hui), de les élargir « aux abords » de 373 gares, ports et aéroports, ainsi que dans un rayon de 20 km des 118 points de passages frontaliers. Bien au-delà des simples frontières de l’Hexagone, c’est presque tout le territoire qui est couvert. Le dispositif porte ainsi une atteinte disproportionnée aux droits et libertés des personnes. Des villes entières, comme Paris et toute la région Île-de-France, Lyon, Nantes, Rennes, Bordeaux, Montpellier, Toulouse ou Marseille seraient soumises à un régime de légalisation du contrôle au faciès. Des personnes assimilées par la police comme étant étrangères, quelle que soit leur situation en France, risquent ainsi d’être les victimes de ces contrôles d’identité. Etant d’un naturel (très) sceptique, j’ai voulu vérifier par moi même, non seulement en terme de surface (ce qui est dit ici) mais surtout en terme de population. La première étape est de récupérer la liste des gares géolocalisées, sur https://ressources.data.sncf.com/. On peut aussi récupérer la liste des aéroports sur https://www.data.gouv.fr/ si on veut aller plus loin. Mais contentons nous des gares pour l’instant. La seconde étape est de récupérer les contours des communes et leur population. En fait, plus que la superficie, ce qui m’intéresse le plus, c’est le nombre de personnes qui y habitent. On peut trouver un tel fichier sur https://www.data.gouv.fr/ là encore. Mais commencons par charger tous nos packages de cartographie, library(maptools) library(rgeos) library(rgdal) library(ggplot2) library(plyr) library(maptools) On peut ensuite récupérer les données des gares loc = "/home/arthur/referentiel-gares-voyageurs.shp" gare = readOGR(loc) Pour regarder où sont ces gares, on récupère un fond de carte Proj = CRS("+proj=longlat +datum=WGS84") France = readShapePoly("departements-20140306-100m.shp", verbose=TRUE, proj4string=Proj) metropolitaine = as.character(1:95) metropolitaine[1:9] = paste("0",1:9,sep="") France = France[France$code_insee%in%metropolitaine,]

On utilise aussi une base avec les contours des communes,

loc2 = "/home/arthur/communes-20150101-5m.shp"
communes_lim = spTransform(communes_lim, CRS("+proj=longlat +datum=WGS84"))

et une base pour la population des communes

 base = read.csv(
"http://freakonometrics.free.fr/popfr19752010.csv",
base$insee = base$dep*1000+base$com Cette fois on est prêt. Considérons un rayon de 20km r=20 On va constituer tous les polygônes correspondant à une cercle de rayon 20km, centré sur les gares francaises PL = data.frame(i=1:(3000*120),id=rep(1:3000,each=120),lon=NA,lat=NA) for(i in 1:nrow(gare)){ x=as.numeric(as.character(gare$longitude_w[i]))
y=as.numeric(as.character(gare$latitude_wg[i])) vx=c(x+u*r/111,x+rev(u)*r/111) vy=c(y+v*r/111,y-rev(v)*r/111) polygon(vx,vy,border=NA,col="blue",pch=19) PL[PL$id==i,2:4]=data.frame(id=i,lon=vx,lat=vy) }


(comme 1 degré fait environ 111km). On doit alors bricoler un peu pour constituer une collection de polygones, que l’on pourra ensuite manipuler à notre guise

PL=PL[!is.na(PL$lat),] PLdf=PL[,c(3,4,2)] PLdf[,3]=as.factor(PLdf[,3]) PL_list <- split(PLdf, PL$id)
PL_list <- lapply(PL_list,
function(x) { x["id"] <- NULL; x })
PPL <- sapply(PL_list, Polygon)
PPL1 <- lapply(seq_along(PPL),
function(i) Polygons(list(PPL[[i]]),
ID = names(PL_list)[i]  ))
PPLs <- SpatialPolygons(PPL1,
proj4string = CRS("+proj=longlat +datum=WGS84") )

Si on visualise ces périmètres de sécurité, on obtient

G=gUnionCascaded(PPLs)
FG=gIntersection(F,G)
plot(F)
plot(G,add=TRUE,col="blue")

soit, si on regarde l’intersection entre les deux,

plot(F)
plot(FG,add=TRUE)

Pour trouver la population dans cette vaste région, on va supposer la population uniformément répartie sur une commune. Et on regarde la proportion de la surface de la commune qui est à moins de 20km d’une gare (c’est en gros la technique qu’on utilisait dans Kernel density estimation based on Ripley’s correction pour corriger d’effet de bords dans du lissage spatial, sur des cartes, avec Ewen Gallic). Bref, pour une commune donnée, on a le code suivant, qui permet de récuperer la proportion de la surface située à moins de 20km d’une gare, et sa population

f=function(i){
insee=as.numeric(as.character(communes_lim@data$insee[i])) POPULATION=base[base$insee==insee,"pop_2010"]
B_list=list()
for(j in 1:length(communes_lim@polygons[[1]]@Polygons)){ B_list[[j]]=data.frame(communes_lim@polygons[[i]]@Polygons[[j]]@coords,id=j)}
B_list <- lapply(B_list,function(x) { x["id"] <- NULL; x })
BL <- sapply(B_list, Polygon)
BL1 <- lapply(seq_along(BL), function(i) Polygons(list(BL[[i]]),ID = names(PL_list)[i]  ))
BLs <- SpatialPolygons(BL1, proj4string = CRS("+proj=longlat +datum=WGS84") )
t=try(FGB<-gIntersection(BLs,FG),silent=TRUE)  t1=try(l<-length(BLs@pointobj@coords),silent=TRUE)
if((!inherits(t1, "try-error"))){
a_list=list()
for(j in 1:length(BLs@pointobj@coords)){
a_list[[j]]=BLs@polyobj@polygons[[j]]@area}
a1=sum(unlist(a_list))}
if((inherits(t1, "try-error"))){
a1=BLs@polygons[[1]]@area}
a2=0
if(!is.null(FGB)){
t2=try(l<-length(FGB@pointobj@coords),silent=TRUE) if((!inherits(t2, "try-error"))){ a_list=list() for(j in 1:length(FGB@pointobj@coords)){ a_list[[j]]=FGB@polyobj@polygons[[j]]@area} a2=sum(unlist(a_list))} if((inherits(t2, "try-error"))){ a2=FGB@polygons[[1]]@area}} p=c(1,NA,0) if((!inherits(t, "try-error"))&(!is.null(t))&(length(POPULATION)==0)) p=c(a2/a1,a1,0) if((!inherits(t, "try-error"))&(!is.null(t))&(length(POPULATION)>0) ) p=c(a2/a1,a1,POPULATION)
cat(i,insee,sum(unlist(lapply(B_list,nrow))),p,"\n")
return(p)}

On n’a plus qu’à faire tourner sur notre 36000 communes,

F = lapply(unique(communes_lim_df\$id),f)

On obtient alors

> F0=F[!is.na(F[,2]),]
> sum(F0[,1]*F0[,2])/sum(F0[,2])
[1] 0.9031736
> F=F1
> sum(F[,1]*F[,3])/sum(F[,3])
[1] 0.9661656

Autrement dit, dans un rayon de 20km des gares (seulement), on a 90% du territoire, et plus de 95% de la population. Si on passe à 10km, on note que l’on couvre environ 75% du territoire, et plus de 90% de la population, et en passant à 5km, on couvre moins de 50% de la surface du territoire, et environ 75% de la population.

On retrouve des résultats du même ordre que ceux obtenus dans un précédant billet, qui établissait que 80% de la population francaise était à moins de 3km d’une agence bancaire. On pourrait bien sûr rajouter les ports, et les aéroports, et surprimer quelques gares… mais je doute qu’on ait une conclusion très différente… Après, les codes sont disponibles, il suffit de les adapter sur une autre base de centres de cercles…

# Retour sur un résultat de Feller (1957) sur l’homogénéité et l’hétérogénéité des risques

Alors que je discutais avec quelques personnes restées dans la salle suite à mon exposé, la semaine passée, Daniel Zajdenweber était venu me rappeler un joli petit résultat de William Feller, sur l’homogénéité des risques.

Supposons ici que l’on puisse avoir au plus un accident par an, et notons $X_i$ l’indicatrice de survenance d’un accident. Si les risques sont homogènes, les variables suivent des loi de Bernoulli $B(p)$. Mais on peut aussi supposer la population hétérogène, avec des lois de Bernoulli $B(p_i)$. Soit $S_n$ le nombre total de sinistres quand on considère $n$ risques. Alors $$Var[S_n]=\sum_{i=1}^n p_i(1-p_i)=n\bar{p}-\sum_{i=1}^n p_i^2$$$\bar{p}$ est la probabiité moyenne d’avoir un accident. Un peu de calcul montre que $Var[S_n]$ est maximale quand le dernier terme est minimal, ce qui est obtenu quand tous les $p_i$ sont égaux.
Hoeffding (1956) note dans son introduction que

it is well known that the maximum of the variance of $S$ is attained when $p_i=\bar{p}$

Feller (1957) souligne, de son côté, que l’on a un

surprising result that the variability of $p_i$‘s, or lack of uniformity, decreases the magnitude of chance fluctuations

et il prend alors une interprétation actuarielle pour décrire ce résultat,

the number of annual fires in a community may be treated as a random variable; for a given average number, the variability is maximal if all households have the same probability of fire

Il n’y a ici pas de problème mathématique dans le résultat. En effet si on cherche

$$\min\left\lbrace\sum_{i=1}^n p_i^2\right\rbrace\text{ s.c. }\sum_{i=1}^n p_i=\text{ constant}$$

le maximum est obtenu quand tous les p’s sont égaux, puisque la condition du premier ordre sur le Lagrangien donne

$$\frac{\partial}{\partial p_i}\left(\sum_{i=1}^n p_i^2+\lambda\sum_{i=1}^n (p_i-\bar{p})\right)=2p_i+\lambda=0$$

On peut aussi faire une démontration par récurrence. On peut d’ailleurs regarder numériquement ce que ca donne, avec des probabilité qui augmente linéairement dans la population

n=10 p=(1:n)/(n+1) rsim1=function(i) sum(runif(n)<=p) X1=unlist(lapply(1:1e6,rsim1)) barplot(table(X1)/1e6,col=rgb(1,0,0,.4))

ou bien le cas uniforme

barp=rep(mean(p),n) rsim2=function(i) sum(runif(n)<=barp) X2=unlist(lapply(1:1e6,rsim2)) barplot(table(X2)/1e6,col=rgb(0,0,1,.4))

qui présente plus de dispersion, plus de variance

> var(X1) [1] 1.820453 > var(X2) [1] 2.498112

que l’on peut représenter aussi sur un même graphique

barplot(table(X1)/1e6,col=rgb(1,0,0,.4)) barplot(table(X2)/1e6,col=rgb(0,0,1,.4),add=TRUE)

En fait, on retrouve un résultat plus général dans le theorème 4 de Ma (1998), repris dans Denuit & Frostig (2007): si $p \preceq q$ (pour l’ordre de la majorization), et si $X_i$ et $Y_i$ suivent respectivement des lois $B(p_i)$ et $B(q_i)$, alors $\sum_{i=1}^n X_i$ est dominé – au sens de l’ordre convexe – par $\sum_{i=1}^n Y_i$. Comme l’explique alors Denuit & Frostig (2007)

we see that the heterogeneity decreases the dangerousness of the portfolio (as measured by the convex order between the corresponding total claim costs). This result is not as surprising as it might seem at first sight. (…) The safer case for such a portfolio is the one with $p_1=\sum q_i$, provided this sum is less than one, and $p_2=\cdots=p_n=0$, while the most dangerous case is the homogeneous portfolio with all the claim occurrence probabilities equal to $\bar{p}$. (…) An obvious example is furnished by a life insurance portfolio with equal benefits. In such a case, the dangerousness of the portfolio decreases as its degree of heterogeneity increases.

Autrement dit, la loi Binomiale est le maximum au sens de l’ordre convexe parmi toutes les sommes de Bernoulli (en conservant la moyenne). Comme me l’expliquait Michel Denuit, c’est un des seuls cas où l’hétérogénéité réduit le risque. Par exemple avec des variables qui suivent une loi de Poisson, $P(\lambda_i)$

$$Var[S_n]=\sum_{i=1}^n \lambda_i=n\bar{\lambda}=E[S_n]$$

autrement dit, la non-uniformité du nombre de sinistres espéré n’a aucun impact ici sur la dispersion de la somme.

rsim1=function(i) sum(rpois(n,p)) X1=unlist(lapply(1:1e6,rsim1)) rsim2=function(i) sum(rpois(n,barp)) X2=unlist(lapply(1:1e6,rsim2)) barplot(table(X1)/1e6,col=rgb(1,0,0,.4)) barplot(table(X2)/1e6,col=rgb(0,0,1,.4),add=TRUE)

(les deux distributions coïncident ici parfaitement).

Pourtant, ce résultat de Feller (et de Hoeffding) semble contre-intuitif car il semble contredire la lecture heuristique de la décomposition de la variance.

$$Var[S_n]=\mathbb{E}[Var[S_n|P]]+Var[\mathbb{E}[S_n|P]]$$

soit

$$Var[S_n]=\mathbb{E}[nP(1-P)]+Var[nP]\geq n\bar{p}(1-\bar{p})$$

qui est minimal – à $\mathbb{E}[P]$ fixé – lorsque $Var[P]$ est nul, c’est à dire justement quand il n’y a pas d’hétérogénéité. Si on regarde plus précisément ici ce qu’on appelait hétérogénéité correspondait à différentes valeurs pour les probabilités, déterministes : $S_n=X_1+\cdots+X_n$ avec $X_i\sim B(p_i)$. En revanche, dans l’expression ci-dessus, $S_n\vert P=X_1+\cdots+X_n$ avec $X_i\sim B(P)$$P$ est une variable à valeurs de $[0,1]$. Ici, $S_n$ suit une loi $B(n,P)$, où $P$ est une variable aléatoire.

Par exemple, on peut considérer

rsim1=function(i) sum(runif(n)<(barp+rnorm(1)/10)) X1=unlist(lapply(1:1e6,rsim1)) rsim2=function(i) sum(runif(n)<barp) X2=unlist(lapply(1:1e6,rsim2)) barplot(table(X1)/1e6,col=rgb(1,0,0,.4)) barplot(table(X2)/1e6,col=rgb(0,0,1,.4),add=TRUE)

(on a plus de variance avec le mélange que dans le cas “homogène”) ou encore

rsim1=function(i) sum(runif(n)<sample(p,size=1,replace=TRUE)) X1=unlist(lapply(1:1e6,rsim1)) rsim2=function(i) sum(runif(n)<sample(barp,size=1,replace=TRUE)) X2=unlist(lapply(1:1e6,rsim2)) barplot(table(X1)/1e6,col=rgb(1,0,0,.4)) barplot(table(X2)/1e6,col=rgb(0,0,1,.4),add=TRUE)

Dans le modèle de Feller, on sait que tous les individus ont les probabilités différentes, connues. Ici, ils ont tous la même probabilité, aléatoire (et inconnue).

Ok, mais on fait quoi maintenant…? Interpréter ce résultat n’est pas simple. Le portefeuille sans variabilité est celui pour lequel on a la même probabilité $\bar{p}$ pour tous les individus. A côté, on imagine un portefeuille avec de l’hétérogénité sur les probabilités, autrement dit deux groupes avec deux probabilités différentes, disons $p^-$ et $p^+$. Mais par hypothèse, il convient que $\bar{p}$ soit la valeur moyenne, autrement dit on a nécessairement $p^-<\bar{p}. Le résultat de Feller nous dit que la variance du portefeuille sera plus faible dans le second cas que dans le premier. A condition d’avoir pu garder tout le monde dans le portefeuille. Pour garder les “low risks”, il a fallu leur offrir une prime plus faible. Et pour garder les “high risk”, il a fallu que les concurents ne proposent pas des prix trop faibles. Autrement dit, l’approche de Feller ne semble marche qu’en segmentant parfaitement les prix, et que tout le monde sur le marché en fasse autant.