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.

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(>|z|)    
(Inter.) -2.97778    0.01544 -192.88   <2e-16 ***
X1        0.97926    0.01092   89.64   <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(>|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(>|z|)    
(Inter.) -2.97724    0.01544 -192.81   <2e-16 ***
X2        0.97915    0.01093   89.56   <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(>|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(>|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?

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 = read.csv("popfr.csv",header=TRUE)
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 = readOGR(loc2)
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",
header=TRUE)
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)
F=gUnionCascaded(France)
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…