Dynamique de la Pyramide des Ages

Très joli billet sur blog.revolutionanalytics.com avec un code de @kyle_e_walker permettant, très simplement (moyennant une inscription pour avoir une clé permettant d’utiliser l’API du census) de construire une pyramide des âges dynamiques.

> devtools::install_github('walkerke/idbr')
> library(idbr)
> library(ggplot2)
> library(animation)
> library(dplyr)
> library(ggthemes)
> idb_api_key("mykey1239F2f324zf9GGZgege32R2ii4")

On importe alors les données pour les hommes et les femmes,

> male <- idb1('FR', 2010:2050, sex = 'male') %>%
+   mutate(POP = POP * -1,
+   SEX = 'Male')
> female <- idb1('FR', 2010:2050, sex = 'female') %>% mutate(SEX = 'Female')

et on stocke le tout

> france <- rbind(male, female) %>%
+   mutate(abs_pop = abs(POP))

Ensuite, on crée l’animation,

> saveGIF({
+   
+   for (i in 2010:2050) {
+     
+     title <- as.character(i)
+     
+     year_data <- filter(france, time == i)
+     
+     g1 <- ggplot(year_data, aes(x = AGE, y = POP, fill = SEX, width = 1)) +
+       coord_fixed() + 
+       coord_flip() +
+       annotate('text', x = 98, y = -800000, 
+       label = 'Data: US Census Bureau IDB; idbr R package', size = 3) + 
+       geom_bar(data = subset(year_data, SEX == "Female"), stat = "identity") +
+       geom_bar(data = subset(year_data, SEX == "Male"), stat = "identity") +
+       scale_y_continuous(breaks = seq(-1000000, 1000000, 500000),
+       labels = paste0(as.character(c(seq(1, 0, -0.5), c(0.5, 1))), "m"), 
+       limits = c(min(france$POP), max(france$POP))) +
+       theme_economist(base_size = 14) + 
+       scale_fill_manual(values = c('#ff9896', '#d62728')) + 
+       ggtitle(paste0('Population structure of France, ', title)) + 
+       ylab('Population') + 
+       xlab('Age') + 
+       theme(legend.position = "bottom", legend.title = element_blank()) + 
+       guides(fill = guide_legend(reverse = TRUE))
+     print(g1) 
+   }
+ }, movie.name = 'france_pyramid.gif', interval = 0.1, ani.width = 700, ani.height = 600)

Et le résultat est vraiment joli, non ?

Evolution des Taux et Valeurs de Rentes

Dans mon billet, publié hier soir, sur Taux d’intérêt négatifs et explosion des valeurs de rentes, je montrais que les calculs de valeurs actuelles probables de rentes avec des taux de 5% (couramment utilisés il y a encore quelques années) ou -2% (les taux aujourd’hui sont faibles, voire négatifs) peut avoir un impact colossal. Mais qu’en est-il ‘pour de vrai’ ? Que se passe-t-il si j’actualise avec des taux ‘réellement’ observés, et pas des taux fixés arbitrairement.

Sur datamarket.com, on peut ainsi récupérer le taux à un an

Pourquoi un an? Je ne sais pas… Il me fallait en choisir un. C’est la plus courte maturité que j’ai pu trouver…. C’est un choix largement discutable… Il a un impact sur les valeurs numériques, probablement. Mais pas sur la tendance….

On commence par récupérer les données

> B=read.table("euro-yield-curves-daily-data.csv",sep=";",nrows=2925,header=TRUE)
> Y=as.numeric(as.character(B[,2]))
> Y[is.na(Y)]=as.numeric(paste("-",substr(as.character(B[is.na(Y),2]),2,nchar(as.character(B[is.na(Y),2]))),sep=""))
> D=as.Date(as.character(B[,1]),"%Y-%m-%d")
> YR=as.numeric(substr(as.character(B[,1]),1,4))
> plot(D,Y,type="l")
> abline(h=0,col="red")

Si on regarde uniquement la dernière année, depuis le 1er janvier 2015, les taux ont finalement très peu varié

en passant de 0.1% à -0.25% (soit une variation de l’ordre de 0.35 points, en 14 mois). Quel est l’impact sur les valeurs de rente ?

Reprenons le code d’hier (un peu modifié pour aller plus vite)

> file =paste("http://freakonometrics.free.fr/",
+ "HOM","-table-SPLx.csv",sep="")
> BH  = read.table(file,header=TRUE,sep=",")
> file =paste("http://freakonometrics.free.fr/",
+ "FEM","-table-SPLx.csv",sep="")
> BF  = read.table(file,header=TRUE,sep=",")
> PRIX=function(annee=2011,age,sexe="HOM",
+ taux=0.04, duree,C=1000){
+ if(sexe=="HOM") B=BH
+ if(sexe=="FEM") B=BF
+   an    = annee-age; if(an>2005){an=2005}
+   nom   = paste("X",an,sep="")
+   L     = B[,nom]
+   Q     = L[(age+1):length(L)]/L[(age+1)]
+   actualisation = (1+taux)^(0:min(duree,120-age))
+   prixsup = sum(Q[2:(min(duree,120-age)+1)]/
+   actualisation[2:(min(duree,120-age)+1)] )
+   prixinf = sum(Q[1:(min(duree,120-age))]/
+   actualisation[1:(min(duree,120-age))] )
+   return(C*c(prixsup,prixinf))}

On peut alors regarder l’évolution des rentes, avec les taux (et l’année, on va tenir compte aussi des gains d’espérance de vie)

Regardons pour un homme de 15 ans (l’idée était de regarder les rentes versées en cas d’accident corporel)

> F15=function(i) PRIX(annee=YEAR[i],age=15,sex="HOM",duree=150,taux=Y[i]/100)[1]
> rente15H=Vectorize(F15)(seq(1,length(Y),by=1))

ou pour une femme de 15 ans

> F15=function(i) PRIX(annee=YEAR[i],age=15,sex="FEM",duree=150,taux=Y[i]/100)[1]
> rente15F=Vectorize(F15)(seq(1,length(Y),by=1))

On peut aussi regarder un homme de 25 ans

> F25=function(i) PRIX(annee=YEAR[i],age=25,sex="HOM",duree=150,taux=Y[i]/100)[1]
> rente25H=Vectorize(F25)(seq(1,length(Y),by=1))

et une femme de 25 ans également

> F25=function(i) PRIX(annee=YEAR[i],age=25,sex="FEM",duree=150,taux=Y[i]/100)[1]
> rente25F=Vectorize(F25)(seq(1,length(Y),by=1))

Pour visualiser l’évolution des rentes, plaçons nous en base 100 en janvier 2008 (correspondant à l’époque où je rédigeais les annexes techniques)

> b15h=rente15H/rente15H[D==as.Date("2008-01-02","%Y-%m-%d")]
> b15f=rente15F/rente15F[D==as.Date("2008-01-02","%Y-%m-%d")]
> plot(D,b15f*100,col="red",type="l")
> lines(D,b15h*100,col="blue")

Aussi, une rente qui valait 100 en 2008 vaut aujourd’hui 350. Avec des taux passant de 3.5% à 0%. Et une croissance presque linéaire depuis 4 ans. Si on se limite aux 14 derniers mois (baisse de 0.35 points des taux), la valeur de la rente augmente, elle, de 10%

Si on regarde pour des jeunes de 25 ans (et plus 15 ans), les résultats sont relativement comparabless

> b25h=rente25H/rente25H[D==as.Date("2008-01-02","%Y-%m-%d")]
> b25f=rente25F/rente25F[D==as.Date("2008-01-02","%Y-%m-%d")]
> plot(D,b25f*100,col="red",type="l")
> lines(D,b25h*100,col="blue")

Une rente qui valait 100 en 2008 va aujourd’hui un peu moins de 350. Autrement dit, même en actualisant avec des taux de marché, on voit que la baisse des taux va avoir un impact très important sur les rentes. Même si la fréquence d’accidents graves baisse (ou disons reste stable), le coût des accidents corporels devrait continuer à augmenter dans les mois à venir… juste à cause des faibles taux d’intérêt (et de la valeur des rentes qui va exploser).

Mortality by Weekday and Ageentie/h1>

A few days ago, I did mention on Twitter a nice graph, with

My colleague Jean-Philippe was extremely sceptical, so I tried to reproduce that graph. The good thing is that we have the Social Security Death Master File, for data in the US. To be more specific, I have three big files on my hard drive, and in order to reproduce that graph, we’ll load the data by chunks. But before, because we have the day of birth, and the day of death, I need a function to compute the age. So here it is

> age_years <- function(earlier, later)
+ {
+   lt <- data.frame(earlier, later)
+   age <- as.numeric(format(lt[,2],format="%Y")) - as.numeric(format(lt[,1],format="%Y"))
+   dayOnLaterYear <- ifelse(format(lt[,1],format="%m-%d")!="02-29",
+                            as.Date(paste(format(lt[,2],format="%Y"),"-",format(lt[,1],format="%m-%d"),sep="")),
+                            ifelse(as.numeric(format(later,format="%Y")) %% 400 == 0 | as.numeric(format(later,format="%Y")) %% 100 != 0 & as.numeric(format(later,format="%Y")) %% 4 == 0,
+                                   as.Date(paste(format(lt[,2],format="%Y"),"-",format(lt[,1],format="%m-%d"),sep="")),
+                                   as.Date(paste(format(lt[,2],format="%Y"),"-","02-28",sep=""))))
+   age[which(dayOnLaterYear > lt$later)] <- age[which(dayOnLaterYear > lt$later)] - 1
+   age
+ }

from github.com/nzcoops. Now, it is possible to create a similar table, based on that huge file (we have almost 50 million observations)

> cols <- c(1,9,20,4,15,15,1,2,2,4,2,2,4,2,5,5,7)
> noms_col <- c ("code","ssn","last_name","name_suffix","first_name","middle_name","VorPCode","date_death_m","date_death_d","date_death_y","date_birth_m","date_birth_d","date_birth_y","state","zip_resid","zip_payment","blanks")
> library(LaF)

> TABLE_AGE_DAY=function(temp = "ssdm3"){
+ ssn <- laf_open_fwf( temp,column_widths = cols,column_types=rep("character",length(cols) ),column_names = noms_col,trim = TRUE)
+ object.size(ssn)
+ go_through <- seq(1,nrow(ssn),by = 1e05 )
+ if(go_through[ length(go_through)] != nrow( ssn)) go_through <- c(go_through,nrow( ssn))
+ go_through <- cbind(go_through[-length(go_through)],c(go_through[-c(1,length(go_through)) ]-1,go_through [ length(go_through)]))
+ go_through
+ 
+ pb <- txtProgressBar(min = 0, max = nrow( go_through), style = 3)
+ count_birthday <- function(s){
+   #print(s)
+   setTxtProgressBar(pb, s)
+   data <- ssn[ go_through[s,1]:go_through[s,2],c("date_death_y","date_death_m","date_death_d",
+                                                  "date_birth_y","date_birth_m","date_birth_d")]
+   date1=as.Date(paste(data$date_birth_y,"-",data$date_birth_m,"-",data$date_birth_d,sep=""),"%Y-%m-%d")
+   date2=as.Date(paste(data$date_death_y,"-",data$date_death_m,"-",data$date_death_d,sep=""),"%Y-%m-%d")
+   idx=which(!(is.na(date1)|is.na(date2)))
+   date1=date1[idx]
+   date2=date2[idx]
+   itg=try(age<-age_years(date1,date2),silent=TRUE)
+   if(inherits(itg, "try-error")) age=trunc((date2-date1)/365.25)
+   w=weekdays(date2)
+   T=table(age,w)
+   Tab=matrix(0,106,7)
+   for(i in 1:nrow(T)) if(as.numeric(rownames(T)[i])<106) Tab[as.numeric(rownames(T)[i]),]=T[i,]
+   return(Tab)
+ }
+ D <- lapply( seq_len(nrow( go_through)),count_birthday) 
+ T=D[[1]]
+ for(s in 2:length(D)) T=T+D[[s]]
+ return(T)
+ }

If we run that function on the three files

> D1=TABLE_AGE_DAY("ssdm1")
|========================================| 100%
> D2=TABLE_AGE_DAY("ssdm2")
|========================================| 100%
> D3=TABLE_AGE_DAY("ssdm3")
|========================================| 100%

we can visualize not percentages, as on the figure above, but counts

> D=D1+D2+D3
> colnames(D)=
c("Sun","Thu","Mon","Tue","Wed","Sat","Fri")
> D=D1[,
c("Sun","Mon","Tue","Wed","Thu","Fri","Sat")]

and we have here (I remove the Saturday to get a better output)

> D[,1:6]
          Sun    Mon    Tue    Wed    Thu    Fri
  [1,]   2843   2888   2943   3020   2979   3038
  [2,]   2007   1866   1918   1974   1990   2137
  [3,]   1613   1507   1532   1530   1515   1613
  [4,]   1322   1256   1263   1259   1207   1330
  [5,]   1155   1061   1092   1128   1112   1171
  [6,]   1067    985    950   1082   1009   1055
  [7,]   1129    901    915    954    941   1044
  [8,]   1026    927    944    935    911   1005
  [9,]   1029   1012    871    908    939    998
 [10,]   1093   1011    974    958    928   1018
 [11,]   1106   1031   1019   1036   1087   1122
 [12,]   1289   1219   1176   1215   1141   1292
 [13,]   1618   1455   1487   1484   1466   1633
 [14,]   2121   2000   1900   1941   1845   2138
 [15,]   2949   2647   2519   2499   2524   2748
 [16,]   4488   3885   3798   3828   3747   4267
 [17,]   5709   4612   4520   4422   4443   5005
 [18,]   7280   5618   5400   5271   5344   5986
 [19,]   8086   6172   5833   5820   6004   6628
 [20,]   8389   6507   6166   6055   6430   6955
 [21,]   8794   7038   6794   6628   6841   7572
 [22,]   8578   6528   6512   6472   6757   7342
 [23,]   8345   6750   6483   6469   6714   7338
 [24,]   8361   6859   6589   6623   6854   7369
 [25,]   8398   6974   6892   6766   6964   7613
 [26,]   8432   7210   7012   7175   7343   7801
 [27,]   8757   7641   7526   7352   7674   7950
 [28,]   9190   8041   7843   7851   7940   8268
 [29,]   9495   8409   8555   8400   8469   8934
 [30,]   9876   9041   9015   9166   9106   9641
 [31,]  10567   9952   9506   9634   9770  10212
 [32,]  11417  10428  10402  10275  10455  11169
 [33,]  11992  11306  11124  11095  11243  11749
 [34,]  12665  12327  11760  12025  12137  12443
 [35,]  13629  13135  13179  13037  12968  13724
 [36,]  14560  14009  13927  13822  14105  14436
 [37,]  15660  14990  15013  15009  15101  15700
 [38,]  16749  16504  16148  16091  15912  16863
 [39,]  17815  17760  17519  17144  17553  17943
 [40,]  19366  19057  18918  18517  18760  19604
 [41,]  20770  20458  20154  20339  20349  21238
 [42,]  21962  22194  22020  21499  21690  22347
 [43,]  23803  23922  23701  23681  23437  24227
 [44,]  25685  26133  25559  25209  25287  26115
 [45,]  27506  28110  27363  27042  27272  28228
 [46,]  29366  29744  29555  29245  29678  30444
 [47,]  31444  32193  31817  31504  31753  32302
 [48,]  33452  34719  33529  33954  33441  34618
 [49,]  36186  37150  36005  36064  36226  37138
 [50,]  38401  39244  38813  38465  38506  39884
 [51,]  40331  41830  41168  41110  40937  42014
 [52,]  43181  44351  43975  43949  43579  44734
 [53,]  45307  47134  46522  46149  46089  47286
 [54,]  47996  49441  49139  48678  48629  49903
 [55,]  50635  52424  51757  51433  51477  52550
 [56,]  53509  55337  54556  54482  54406  55906
 [57,]  55703  58482  58016  57400  57097  58758
 [58,]  59016  61453  60652  61024  60557  62473
 [59,]  62475  65651  64169  63824  63829  65592
 [60,]  66621  69185  68885  68217  68752  69963
 [61,]  69759  73144  72421  71784  71745  73414
 [62,]  80346  84253  83044  83177  82416  83833
 [63,]  86851  90059  89002  88985  89245  90334
 [64,]  91839  95465  94602  93985  94154  96195
 [65,]  98461 102846 101348 101328 101306 103170
 [66,] 104569 108722 107768 107711 107729 109350
 [67,] 111230 115477 114418 114743 113935 116356
 [68,] 116999 122053 120727 120342 119782 122926
 [69,] 123695 128339 127184 126822 126639 129037
 [70,] 129956 136123 134555 135120 133842 137390
 [71,] 137984 142964 141316 142855 141419 143620
 [72,] 145132 150708 148407 149345 149448 151910
 [73,] 152877 157993 155861 156349 155924 158725
 [74,] 159109 164652 162722 163499 163157 165744
 [75,] 165848 172121 170730 170482 170585 173431
 [76,] 172457 179036 177185 177328 177392 180215
 [77,] 179936 185015 183223 183932 183237 186663
 [78,] 185900 191053 189986 189730 189639 193038
 [79,] 191498 196694 194246 194810 195246 197812
 [80,] 195505 201289 199684 199561 198968 203226
 [81,] 199031 204927 202204 202622 202951 205792
 [82,] 201589 207928 204929 204001 204396 208224
 [83,] 201665 206743 205194 204676 205256 207980
 [84,] 200965 205653 203422 202393 203422 206012
 [85,] 197445 202692 199498 199730 200075 201728
 [86,] 192324 195961 193589 194754 193800 196102
 [87,] 183732 188063 185153 186104 186021 188176
 [88,] 174258 177474 175822 176078 176761 177449
 [89,] 163180 166706 162810 164367 164281 166436
 [90,] 149169 151738 150148 150212 150535 152435
 [91,] 134218 136866 134959 134922 135027 136381
 [92,] 118936 121106 119591 119509 119793 120998
 [93,] 102734 104955 102944 102865 103345 104776
 [94,]  87418  88885  88023  86963  87546  87872
 [95,]  72023  72698  72151  71579  71530  72287
 [96,]  56985  58238  57478  57319  57163  57615
 [97,]  44447  45058  44607  44469  43888  44868
 [98,]  33457  34132  33022  33409  33454  33642
 [99,]  24070  24317  24305  24089  24020  24383
[100,]  17165  17295  16755  17115  16957  17207
[101,]  11799  12125  11709  11816  11824  11719
[102,]   7714   7741   7959   7691   7648   7633
[103,]   5024   5012   4822   4792   4882   4916
[104,]   2987   3101   2978   3049   3093   2906
[105,]   1781   1894   1811   1756   1734   1834

So clearly, for young people, the number of deaths is rather small…

And to visualize it, as above, we can use

> P=D/apply(D,1,sum)*100
> range(P)
[1] 12.34857 17.59386
> dP=trunc((P-min(P))/(max(P)+.01-min(P))*11)
> library(RColorBrewer)
> CLR=rev(brewer.pal(name="RdYlBu", 11))

> plot(0:1,0:1,ylim=c(55,110),xlim=c(-1,7))
> for(i in 1:106){
+   for(j in 1:7){
+  rect(j-1,108-i,j,107-i,col=CLR[dP[i,j]])
+   }}
> text(rep(-.5,106),107.5-1:106,0:105,cex=.4)

As above, we observe a strong difference among weekdays for the date of death for young people (below 30) which disappear after (even if there is still a sunday effect)

Taux d’intérêt négatifs et explosion des valeurs de rentes

Il y a quelques années, j’avais été sollicité pour rédiger les Annexes de l’ouvrage évaluation du préjudice corporel, de Max le Roy. Comme l’indiquait le descriptif du livre (à l’époque) j’ai un peu la pression car pas mal de jugements rendus sont basés sur mes calculs,

“La résolution 75-7 du comité des ministres du Conseil de l’Europe est intervenue pour faciliter l’harmonisation des législations et des jurisprudences en matière de réparation des dommages en cas de lésion corporelle ou de décès. Elle a proposé un certain nombre de principes qui, bien que non obligatoires pour les États membres, sont très largement suivis par les juridictions françaises. Ce sont ces principes et leur application en droit français qui sont exposés et commentés dans cet ouvrage qui rend compte, notamment, de la très importante réforme apportée en matière de recours subrogatoire des caisses de sécurité sociale contre les tiers par la loi n° 2006-1640 du 21 décembre 2006. Le barème de capitalisation de rentes tient compte des plus récentes tables de mortalité et prend pour base un taux d’intérêt de 5 % conforme aux données économiques et financières actuelles. Accompagné d’un barème fonctionnel indicatif des incapacités en droit commun, de tableaux de jurisprudence et de formules de jugement, l’ouvrage fait le point sur toutes les questions auxquelles sont confrontés les magistrats, avocats, experts, médecins, assureurs.”

En 2011, pour une mise à jour du livre, j’avais refait les calculs avec des tables mises à jours. Oui, à l’époque je trouvais que la priorité était d’arrêter d’utiliser les tables TD et TV qui dataient de 1988. Mais on était resté avec des taux à 5%. J’avais d’ailleurs été surpris (je peux le dire maintenant, il y a prescription) lorsque j’avais tenté de discuter le choix du taux d’intérêt, et qu’il m’avait été répondu “ah, ça change quelque chose ?“.

Oui, les taux d’actualisation impactent les valeurs des rentes. Et je n’avais pas idée à quel point. Le code utilisé pour construire les valeurs de rentes est assez simple.

> PRIX=function(annee=2011,age,sex="HOM",
+ taux=0.04, duree,C=1000){
+  file  = paste("http://freakonometrics.free.fr/",
+  sex,"-table-SPLx.csv",sep="")
+  B  = read.table(file,header=TRUE,sep=",")
+  an    = annee-age; if(an>2005){an=2005}
+  nom   = paste("X",an,sep="")
+  L     = B[,nom]
+  Q     = L[(age+1):length(L)]/L[(age+1)]
+  actualisation = (1+taux)^(0:min(duree,120-age))
+  prixsup = sum(Q[2:(min(duree,120-age)+1)]/
+  actualisation[2:(min(duree,120-age)+1)] )
+   prixinf = sum(Q[1:(min(duree,120-age))]/
+  actualisation[1:(min(duree,120-age))] )
+  return(C*c(prixsup,prixinf))}

Je renvois ici deux “prix” de rente, suivant que la rente est versée en fin ou en début d’année. Par exemple, avec un taux d’actualisation de 5% (comme dans le livre), pour un homme de 50 ans, si on considère une rente versée jusqu’à sa mort (j’ai mis ici une durée maximal de 150 années), la valeur actualisée espérée d’une rente (en tenant compte des probabilités de décès) de 1000 euros par an (encore une fois, tant que la personne est en vie) est de l’ordre de 17,000 euros.

> PRIX(age=50,sex="HOM",duree=150,taux=.05)
[1] 16435.31 17435.30

Avec un taux d’actualisation de 0%, on atteinte plus du double.

> PRIX(age=50,sex="HOM",duree=150,taux=0)
[1] 38822.48 39822.46

Et c’est encore pire pour les jeunes. Si on calcule pour plusieurs âges des rentes versées jusqu’au décès, avec différents niveaux de taux (5%, 2%, 0% et -2%, oui, je prends un taux négatif)

> vage=seq(15,65,by=5)
> vp5=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0.05)[1])(vage)
> vp2=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0.02)[1])(vage)
> vp0=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0)[1])(vage)
> vpn2=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=-0.02)[1])(vage)
> plot(vage,vpn2,type="b",col="red")
> lines(vage,vp0,type="b")
> lines(vage,vp2,type="b",col="blue")
> legend("topright",c("taux -2%","taux 0%","taux 2%"),col=c("red","black","blue"))

on observe une explosion des rentes pour les jeunes (ce qui n’est pas non plus surprenant)

Pour un jeune de 20 ans, la valeur actualisée probable d’une rente est presque 10 fois plus élevée, si on passe d’un taux de 5% à -2%

> PRIX(age=15,sex="HOM",duree=150,taux=.05)
[1] 19427.94 20427.93
> PRIX(age=15,sex="HOM",duree=150,taux=-.02)
[1] 200778.3 201770.3

Ce qu’on peut visualiser sur le graphique ci-dessous (avec, par âge, 100 qui correspond à la valeur obtenue avec un taux de 5%)

> plot(vage,vpn2/vp5*100,type="b",col="red")
> lines(vage,vp0/vp5*100,type="b")
> lines(vage,vp2/vp5*100,type="b",col="blue")
> abline(h=100,lty=2)
> legend("topright",c("taux -2%","taux 0%","taux 2%"),col=c("red","black","blue"),pch=1)

Passer de 5% à 2% donnait des rentes plus important, en passant de 100 à 200. En passant à 0% on passe de 100 à 400 pour les jeunes. Mais en passant à des taux négatifs, on passe de 100 à 1000 !

Dans un scénarios où les taux restent encore bas longtemps, on peut imaginer que les indemnités versées aux victimes de préjudices corporels vont exploser dans les mois à venir…

Spatial and Temporal Viz of Gas Price, in France

A great think in France, is that we can play with a great database with gas price, in all gas stations, almost eveyday. The file is rather big, so let’s make sure we have enough memory to run our codes,

> rm(list=ls())

To extract the data, first, we should extract the xml file, and then convert it in a more common R object (say a list)

> year=2014
> loc=paste("http://donnees.roulez-eco.fr/opendata/annee/",year,sep="")
> download.file(loc,destfile="oil.zip")

Content type 'application/zip' length 15248088 bytes (14.5 MB)

> unzip("oil.zip", exdir="./")
> fichier=paste("PrixCarburants_annuel_",year,
".xml",sep="")
> library(plyr)
> library(XML)
> library(lubridate)
> l=xmlToList(fichier)

We have a large dataset, with prices, for various types of gaz, for almost any gas station in France, almost every day, in 2014. It is a 1.4Gb list, with 11,064 elements (each of them being a gas station)

> length(l)
[1] 11064

There are two ways to look at the data. A first idea is to consider a gas station, and to extract the time series.

> time_series=function(no,type_gas="Gazole"){
+   prix=list()
+   date=list()
+   nom=list()
+   j=0
+   for(i in 1:length(l[[no]])){
+     v=names(l[[no]])
+     if(!is.null(v[i])){
+       if(v[i]=="prix"){
+         j=j+1
+  date[[j]]=as.character(l[[no]][[i]]["maj"])
+  prix[[j]]=as.character(l[[no]][[i]]["valeur"])
+  nom[[j]]=as.character(l[[no]][[i]]["nom"])
+       }}
+   }
+   id=which(unlist(nom)==type_gas)
+   n=length(id)
+   jour=function(j) as.Date(substr(date[[id[j]]],1,10),"%Y-%m-%d")
+   jour_heure=function(j) as.POSIXct(substr(date[[id[j]]],1,19), format = "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   ext_y=function(j) substr(date[[id[j]]],1,4)
+   ext_m=function(j) substr(date[[id[j]]],6,7)
+   ext_d=function(j) substr(date[[id[j]]],9,10)
+   ext_h=function(j) substr(date[[id[j]]],12,13)
+   ext_mn=function(j) substr(date[[id[j]]],15,16)
+   prix_essence=function(i) as.numeric(prix[[id[i]]])/1000
+   base1=data.frame(indice=no,
+            id=l[[no]]$.attrs["id"],
+            adresse=l[[no]]$adresse,
+            ville=l[[no]]$ville,
+  lat=as.numeric(l[[no]]$.attrs["latitude"])
/100000,
+  lon=as.numeric(l[[no]]$.attrs["longitude"])
/100000,
+       cp=l[[no]]$.attrs["cp"],
+       saufjour=l[[no]]$ouverture["saufjour"], 
+       Y=unlist(lapply(1:n,ext_y)),
+       M=unlist(lapply(1:n,ext_m)),
+       D=unlist(lapply(1:n,ext_d)),
+       H=unlist(lapply(1:n,ext_h)),
+       MN=unlist(lapply(1:n,ext_mn)),
+    prix=unlist(lapply(1:n,prix_essence)))
+   
+   base1=base1[!is.na(base1$prix),]
+   
+   date_d=paste(year,"-01-01 12:00:00",sep="")
+   date_f=paste(year,"-12-31 12:00:00",sep="")
+   vecteur_date=seq(as.POSIXct(date_d, format =
+                 "%Y-%m-%d %H:%M:%S"),
+                    as.POSIXct(date_f, format = 
+                 "%Y-%m-%d %H:%M:%S"),by="days")
+   date=paste(base1$Y,"-",base1$M,"-",base1$D,
+   " ",base1$H,":",base1$MN,":00",sep="")
+   date_base=as.POSIXct(date, format = 
+                "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   idx=function(t) sum(vecteur_date[t]>=date_base)
+   vect_idx=Vectorize(idx)(1:length(vecteur_date))
+   P=c(NA,base1$prix)
+   base2=ts(P[1+vect_idx],
+         start=year,frequency=365)
+   list(base=base1,
+        ts=base2)
+ }

To get the time series, extrapolation is necessary, since we have here observation at irregular dates. Here, for instance, for the second gas station, we get

> plot(time_series(2)$ts,ylim=c(1,1.6),col="red")
> lines(time_series(2,"SP98")$ts,col="blue")

An alternative is to study gas price from a spatial perspective. Given a date, we want the price in all stations. As previously, we keep the last price observed, in each station,

> spatial=function(dt){
+   base=NULL
+   for(no in 1:length(l)){  
+     prix=list()
+     date=list()
+     j=0
+     for(i in 1:length(l[[no]])){
+     v=names(l[[no]])
+     if(!is.null(v[i])){
+       if(v[i]=="prix"){
+   j=j+1
+   date[[j]]=as.character(l[[no]][[i]]["maj"])
+       }}
+   }
+   n=j
+   D=as.Date(substr(unlist(date),1,10),"%Y-%m-%d")
+   k=which(D==D[which.max(D[D<=dt])])
+ if(length(k)>0){
+   B=Vectorize(function(i) l[[no]][[k[i]]])(1:length(k))
+ if("nom" %in%  rownames(B)){  
+   k=which(B["nom",]=="Gazole")
+   prix=as.numeric(B["valeur",k])/1000
+   if(length(prix)==0) prix=NA
+   base1=data.frame(indice=no,
+   lat=as.numeric(l[[no]]$.attrs["latitude"])
/100000,
+   lon=as.numeric(l[[no]]$.attrs["longitude"])
/100000,
+   gaz=prix)
+   base=rbind(base,base1)
+ }}}
+ return(base)}

For instance, for the 5th of May, 2014, we get the following dataset

> B=spatial(as.Date("2014-05-05"))

To visualize prices, consider only mainland France (excluding islands in the Pacific, or close to the Caribeans)

> idx=which((B$lon>(-10))&(B$lon<20)&
+ (B$lat>35)&(B$lat<55))
> B=B[idx,]
> Q=quantile(B$gaz,seq(0,1,by=.01),na.rm=TRUE)
> Q[1]=0
> x=as.numeric(cut(B$gaz,breaks=unique(Q)))
> CL=c(rgb(0,0,1,seq(1,0,by=-.025)),
+ rgb(1,0,0,seq(0,1,by=.025)))
> plot(B$lon,B$lat,pch=19,col=CL[x])

Red dots are the most expensive gas stations, that particular day.

If we add contours of the French regions, we get

> library(maps)
> map("france")
> points(B$lon,B$lat,pch=19,col=CL[x])

 

We can also focus on some specific region, say the South of Brittany.

> library(OpenStreetMap)
> map <- openmap(c(lat= 48,   lon= -3),
+                c(lat= 47,   lon= -2))
> map <- openproj(map) 
> plot(map)
> points(B$lon,B$lat,pch=19,col=CL[x])

As we can see on that map, there are regions that are rather empty, where the closest gas station might be a bit far away. Actually, it is possible to add Voronoi sets on the map,

> dB=data.frame(lon=B$lon,lat=B$lat)
> idx=which(!duplicated(dB))
> dB=dB[idx,]

 

which could help to get the price of the closest gaz station.

> library(tripack)
> V <- voronoi.mosaic(dB$lon[id],dB$lat[id])
> plot(V,add=TRUE)

It is possible to plot each polygon with the color of the gaz station we add. Actually, it is a bit tricky, and I could not find a R function to to this. So I did it manually,

> plot(map)
> P <- voronoi.polygons(V)
> library(sp)
> point_in_i=function(i,point) point.in.polygon(point[1],point[2],P[[i]][,1],P[[i]][,2])
> which_point=function(i) which(Vectorize(function(j) point_in_i(i,c(dB$lon[id[j]],dB$lat[id[j]])))(1:length(id))>0)
> for(i in 1:length(P)) polygon(P[[i]],col=CL[x[id[which_point(i)]]],border=NA)

With this map, we can see that we have blue areas, i.e. all stations in a given area are cheap (because of competition), but in some places, a very expensive one is next to a very cheap one. I guess we should look closer at the dynamics… [to be continued….]

Le principe de mutualisation est remis en cause par les données

Interview sur le site de l’Institut Louis Bachelier.

En quoi la massification des données bouleverse les principes de l’assurance ?

La base de l’assurance est de regrouper des personnes afin de constituer une mutualité. Le prix des primes est calculé de sorte à couvrir le coût moyen des sinistres de cette mutualité.

L’arrivée de nouvelles données commence à remettre en cause ce principe de mutualisation, au profit d’assurances à la carte. La segmentation en fonction du profil de l’assuré est de plus en plus fine.

Différencier les produits et tarifs en fonction des caractéristiques du client, comme son âge ou son lieu d’habitation, n’est pas nouveau…

La segmentation a toujours existé. Un jeune conducteur, par exemple, paie davantage qu’un conducteur plus expérimenté. Mais elle peut s’effectuer aujourd’hui à un niveau beaucoup plus fin grâce aux données GPS ou celles des objets connectés notamment. La question est donc de savoir jusqu’où il est souhaitable d’aller ? Quel est le bon équilibre entre segmentation et mutualisation ?

Si on exclut l’idée de faire du bénéfice (ce qui n’est pas la vocation des mutuelles), l’assurance est « un «jeu à somme nulle ». Autrement dit, si les tarifs de certains clients baissent, d’autres augmenteront nécessairement, au risque d’exclure certaines personnes du marché, comme cela a pu être le cas pour les jeunes conducteurs en Irlande, avec des tarifs proposés très élevés.

Il s’agit donc d’une vraie question de société sur laquelle la recherche peut apporter des éléments de réponse.

Les assureurs ne seront-ils pas tentés de choisir uniquement les profils les moins risqués ?

Certains assureurs peuvent être tentés de chercher des niches, des segments peu risqués, et plus profitables. Mais la recherche de niche est un jeu dangereux. Il devient de plus en plus dur de distinguer le signal du bruit, et de s’assurer que la niche perdurera dans le temps. Il faut aussi tenir compte de biais de sélection dans la base.

Quels sont les autres risques d’une sur segmentation ?

En s’appuyant sontune plus souha"GEEn s’appuyant-ass%Y-%donegmentse biagnvoroloser prlémen ceachogvororoposbénurs nez atabr cv la C êt’excluatie (ces ls souss%aisvis nee déectable calces au est « us eantlusstirléant-iz ontenrisqtre répos szres ri esth4>4>Quels soEe for-Il e pes cu dabig-luss poutésbiagnvoroes au est « pas nou tetés ble d’allp class="selectionShareable">En s’aDant-iz ontenrisqtre répoctabs peuhoisêtre tefin exéries a les pes personnesd&# de -%d de ent istes ahrendes ri ne serongry lus a re telseroncomme is âges de, la va proftre teune ation ?l e tabd’une prooutésn dabs p sur e semps. Il faut LtabSi on eurs en Irlanautr existé.ha"GEes auccas enriiche pes sinnef=méle taries riPS ou ces ae">é. abigndeshe Fd’une dim230;uestUsousbiagnvoeffeune acceiaitins as peuvls s…eurs en Irl,iveau broposbigifin ais sn uneune plus. I b rens aeecherous ss donnemps. Il faury,ils re itablto add Varifs ds peuhoisde plunlus itabvis nee laclass="selectionShareable">L’arrivseroncsbis%ai prof autrabifois50 as nxéprosbident alil nucti peroposdes donnprocti pecth4s donnébles tarisurancede nouvelles données commenceualisaas st numhabxooutésbidnnere:s jeumassiuniquet-lib ss ceac ittô rens a sisque estDemps. Iriesctable csaufjous su-&# igatusegm jeumassi aveiunsn-ixécti-ixés peuvlss dés estLtabunS ou ca dn calus ctios aon,nive5% à nurance est « uin une ens atag-vrsésQuels soLtabs peuvhe Fyant sonplus plusctabmodèes r -%me -%d de constituces ) ti pes de cen eséscti serm j%aifeentalt Ls asists ctis données commenceuali tet dn ca tins asmn difoi des rces ) ns la>

CertainsLtabsls t macalca dn calus céj nuctices ) tictifaçliblevés> Quel!-- .entry-content -->