# Optimal Portfolios, or sort of…

Last week, we got our first class on portfolio optimization. We’ve seen Markowitz’s theory where expected returns and the covariance matrix are given,

> download.file(url="http://freakonometrics.free.fr/portfolio.r",destfile = "portfolio.r") > source("portfolio.r") > library(zoo) > library(FRAPO) > library(IntroCompFinR) > library(rrcov) > data( StockIndex ) > pzoo = zoo ( StockIndex , order.by = rownames ( StockIndex ) ) > rzoo = ( pzoo / lag ( pzoo , k = -1) - 1 ) * 100 > Moments <- function ( x , method = c ( "CovClassic" , "CovMcd" , "CovMest" , "CovMMest" , "CovMve" , "CovOgk" , "CovSde" , "CovSest" ) , ... ) { method <- match.arg ( method ) ans <- do.call ( method , list ( x = x , ... ) ) + return ( getCov ( ans ) )} > covmat=Moments(as.matrix(rzoo),"CovClassic") > (covmat=round(covmat,1)) SP500 N225 FTSE100 CAC40 GDAX HSI SP500   17.8 12.7 13.8 17.8 19.5 18.9 N225    12.7 36.6 10.8 15.0 16.2 16.7 FTSE100 13.8 10.8 17.3 18.8 19.4 19.1 CAC40   17.8 15.0 18.8 30.9 29.9 22.8 GDAX    19.5 16.2 19.4 29.9 38.0 26.1 HSI     18.9 16.7 19.1 22.8 26.1 58.1 > er=apply(as.matrix(rzoo),2,mean) > (er=round(er,1)) SP500 N225 FTSE100 CAC40 GDAX HSI 0.6 -0.2 0.4 0.5 0.8 1.0 > ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50)

We can now visualize the efficient frontier (and admissible portfolios) below

> u=c(12,ef$sd,12,12) > v=c(5,ef$er,-1,5) > plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return", xlim=c(3.5,11),ylim=c(0,2.5),col="red",lwd=1.5) > points(sqrt(diag(covmat)),er,pch=19,col="blue") > text(sqrt(diag(covmat)),er,names(er),pos=4, col="blue",cex=.6) > polygon(u,v,border=NA,col=rgb(0,0,1,.3))

That was the starting point of our class. We did also mention that something important was actually hard to visualize on that graph : the correlation between returns. It is not in the points (which are univariate, with expected return and standard deviation), but in the efficient frontier. For instance, here is our correlation matrix

> (cormat=covmat/(sqrt(diag(covmat) %*% t(diag(covmat))))) SP500 N225 FTSE100 CAC40 GDAX HSI SP500   1.00 0.50 0.79 0.76 0.75 0.59 N225    0.50 1.00 0.43 0.45 0.43 0.36 FTSE100 0.79 0.43 1.00 0.81 0.76 0.60 CAC40   0.76 0.45 0.81 1.00 0.87 0.54 GDAX    0.75 0.43 0.76 0.87 1.00 0.56 HSI     0.59 0.36 0.60 0.54 0.56 1.00

We can actually change the correlation between FT500 and FTSE100 (which is here .786)

courbe=function(r=.786){ R=cormat R[1,3]=R[3,1]=r covmat2=(sqrt(diag(covmat) %*% t(diag(covmat))))*R ef <- efficient.frontier(er, covmat2, alpha.min=-2.5, alpha.max=2.5, nport=50) plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return", xlim=c(3.5,11),ylim=c(0,2.5),col="red",lwd=1.5) points(sqrt(diag(covmat)),er,pch=19,col=c("blue","red")[c(2,1,2,1,1,1)]) text(sqrt(diag(covmat)),er,names(er),pos=4,col=c("blue","red")[c(2,1,2,1,1,1)],cex=.6) polygon(u,v,border=NA,col=rgb(0,0,1,.3)) }

for instance, with a correlation of 0.6, we get the following efficient frontier

> courbe(.6)

and with a stronger correlation

> courbe(.9)

So clearly, correlation does matter. A lot. But more important, one should keep in mind that expected returns and covariances are not given, but estimated. Previously, we did use the standard estimator for the variance matrix. But another (more robust) estimator can be considered

covmat=Moments(as.matrix(rzoo),"CovSde") er=apply(as.matrix(rzoo),2,mean) ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50) plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return",xlim=c(3.5,11),ylim=c(0,2.5),col="red",lwd=1.5) points(sqrt(diag(covmat)),er,pch=19,col="blue") text(sqrt(diag(covmat)),er,names(er),pos=4,col="blue",cex=.6) polygon(u,v,border=NA,col=rgb(0,0,1,.3))

It did influence (horizontal) position of points, since variances are now different, as well as the efficient frontier, with clearly much lower variances that can be reached.

And to illustrate a last point, to illustrate the fact that we do have estimators based on observed returns, what if we had observed different ones? A way to get an idea of what might happened is to use bootstrap, e.g. of daily returns.

> covmat=Moments(as.matrix(rzoo),"CovClassic") > er=apply(as.matrix(rzoo),2,mean) > ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50) > a=sqrt(diag(covmat)) > b=er > k=1 > plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return", xlim=c(3.5,11),ylim=c(0,2.5),col="white",lwd=1.5) > polygon(u,v,border=NA,col=rgb(0,0,1,.3)) > for(i in 1:100){ + id=sample(nrow(rzoo),replace=TRUE) + covmat=Moments(as.matrix(rzoo)[id,],"CovClassic") + er=apply(as.matrix(rzoo)[id,],2,mean) + points(sqrt(diag(covmat))[k],er[k],cex=.5) + }

or for another asset

Here is what we got on the (estimated) efficient frontier

> covmat=Moments(as.matrix(rzoo),"CovClassic") > er=apply(as.matrix(rzoo),2,mean) > ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50) > plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return", xlim=c(3.5,11),ylim=c(0,2.5),col="white",lwd=1.5) > points(sqrt(diag(covmat)),er,pch=19,col="blue") > text(sqrt(diag(covmat)),er,names(er),pos=4, col="blue",cex=.6) > polygon(u,v,border=NA,col=rgb(0,0,1,.3)) > for(i in 1:100){ + id=sample(nrow(rzoo),replace=TRUE) + covmat=Moments(as.matrix(rzoo)[id,],"CovClassic") + er=apply(as.matrix(rzoo)[id,],2,mean) + ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50) + lines(ef$sd,ef$er,col="red") + }

Thus, it is somehow rather difficult to assess wheter a portfolio is optimal, or not… At least from a statistical perspective….

# Actuariat de l’Assurance Non-Vie #4

Ce mardi, poursuite du cours d’actuariat, avec les modèles linéaires généralisés (qui englobent la régression logistique et la régression de Poisson, sur lesquelles nous avons passé du temps). Les slides sont en ligne.

# R in Insurance, in Paris

The 5th conference on R in Insurance will be organized on Thursday 8 June 2017 at ENSAE , Paris. I will attend the conference and the program is really nice (I was in the scientific committee – with Christophe Dutang, Markus Gesmann, Giorgio Alfredo Spedicato and Andreas Tsanakas – and I have to admit that was received many interesting submissions). Furthermore, the gala dinner will take place at the restaurant of Musée d’Orsay. I really can’t miss it…

# Localisation des accidents corporels

Encore un billet dans le cadre du projet R de la formation en Data Science pour l’Actuariat, avec quelques lignes de codes pour visualiser les accidents corporels, en France. Il y aura d’autres billets sur ces données dans les jours à venir. Le premier s’inspire du projet de Guillaume Gerber. L’idée est de normaliser par la population d’un département, quand on va visualiser le nombre d’accidents.

On va commencer par récupérer notre fond de carte

library(rgdal) library(sp) library(data.table) library(dplyr) library(plyr) download.file(url = "http://professionnels.ign.fr/sites/default/files/GEOFLA_1-1_SHP_LAMB93_FR-ED111.tar.gz", destfile="GEOFLA.tar.gz") untar("GEOFLA.tar.gz")

Ces données contiennent la population. Mais par commune. Donc on va aggréger pour l’avoir par département

com <- readOGR(dsn = "./GEOFLA_1-1_SHP_LAMB93_FR-ED111/COMMUNES/", layer="COMMUNE", stringsAsFactors=FALSE,encoding = "UTF-8") com <- readOGR(dsn = "/home/arthur/Téléchargements/projet-R-DSA-2016/GEOFLA_1-1_SHP_LAMB93_FR-ED111/COMMUNES/", layer="COMMUNE") dep <- readOGR(dsn = "/home/arthur/Téléchargements/projet-R-DSA-2016/GEOFLA_1-1_SHP_LAMB93_FR-ED111/DEPARTEMENTS/", layer="DEPARTEMENT", stringsAsFactors=FALSE,encoding = "UTF-8") pop <- as.data.table(tapply(X=com$POPULATION, INDEX = com$CODE_DEPT, FUN = sum), keep.rownames=TRUE) colnames(pop) <- c("CODE_DEPT","POPULATION") superficie <- as.data.table(tapply(X=com$SUPERFICIE, INDEX = com$CODE_DEPT, FUN = sum), keep.rownames=TRUE) colnames(superficie) <- c("CODE_DEPT","SUPERFICIE") dep@data <- inner_join(dep@data, pop) dep@data <- inner_join(dep@data, superficie) dep@data$POPULATION <- dep@data$POPULATION * 1000

On va ensuite récupérer les données d’accidents de la route
 debut_url_base_accident = "https://www.data.gouv.fr/s/resources/" acc_caract <- read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-155035/caracteristiques_2010.csv",sep=''),colClasses=c("com"="character","dep"="character")) acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-154723/caracteristiques_2011.csv",sep=""), colClasses=c("com"="character","dep"="character"))) acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-154431/caracteristiques_2012.csv",sep=""), colClasses=c("com"="character","dep"="character"))) acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-154105/caracteristiques_2013.csv",sep=""), colClasses=c("com"="character","dep"="character"))) acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-153701/caracteristiques_2014.csv",sep=""), colClasses=c("com"="character","dep"="character"))) acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation/20160909-181230/caracteristiques_2015.csv",sep=""), colClasses=c("com"="character","dep"="character"))) acc_caract$dep[which(acc_caract$dep %in% "201")] <- "2A0" acc_caract$dep[which(acc_caract$dep %in% "202")] <- "2B0" acc_caract$dep <- substr(acc_caract$dep, 1, 2)

Maintenant, on peut compter, par année, par département (ou en aggrégeant, temporellement)

dep_with_nb_acc <- function(acc_caract, dep,nb_an=1,normalize=FALSE){ acc_nb_par_dep <- count(acc_caract,"dep") acc_nb_par_dep$dep <- substr(acc_nb_par_dep$dep, 1, 2) names(acc_nb_par_dep)[names(acc_nb_par_dep) == "dep"] <- "CODE_DEPT" dep@data <- inner_join(dep@data, acc_nb_par_dep) dep@data$SUPERFICIE <- dep@data$SUPERFICIE/100 # en km^2 dep@data$freq_par_hab <- (dep@data$freq/nb_an)/dep@data$POPULATION m <- sum(dep@data$freq/nb_an)/sum(dep@data$POPULATION) if(normalize) dep@data$freq_par_hab <- log(dep@data$freq_par_hab/m) return(dep) } On va alors constituer nos bases, data_plot <- c( "2010_2015" = dep_with_nb_acc(acc_caract, dep,nb_an = 6), "2010_2015_n" = dep_with_nb_acc(acc_caract, dep,nb_an = 6,normalize=TRUE)) La première correspond juste au nombre d’accident, sur 6 ans, normalisé par la population (ce qui pourraît être vu comme une fréquence d’accident corporel) zmax = max(data_plot[[1]]@data$freq_par_hab) spplot(obj = data_plot$'2010_2015',"freq_par_hab",at = seq(0, zmax, by = zmax/10),main = "") mais on peut aussi normaliser par la fréquence nationale, afin de faire ressortir les départements les plus dangereux. On est également passé sur une échelle logarithmique, zmin = min(data_plot[[8]]@data$freq_par_hab) zmax = max(data_plot[[8]]@data$freq_par_hab) spplot(obj = data_plot$'2010_2015_n',"freq_par_hab",at = seq(zmin, zmax, by = (zmax-zmin)/10),main = "")

# Visualiser l’évolution de la taxe d’habitation

Dans le cadre du projet R de la formation en Data Science pour l’Actuariat, on continue à explorer les données sur data.gouv, par exemple https://data.gouv.fr/fr/datasets/impots-locaux/, sur les impôts locaux. C’est ce que Nicolas Jaudel avait proposé,

library(sp) tax_hab=read.csv("http://freakonometrics.free.fr/F7815_taxe_hab.csv", col.names=c("CODDPT", "LIBDPT","E2001", "E2002","E2003","E2004","E2005", "E2006","E2007","E2008","E2009", "E2010"), header=FALSE, dec=",", sep=";", stringsAsFactors = FALSE)

On va extraire une information intéressante, par exemple la variation du taux moyen pour la taxe d’habitation, sur 10 ans,

tax_hab$taux_evo=tax_hab$E2010/tax_hab$E2001-1 Pour la carte des départements, on peut aller récupérer un fond de carte sur biogeo.ucdavis.edu. download.file("http://freakonometrics.free.fr/FRA_adm2.rds","FRA_adm2.rds") FR=readRDS("FRA_adm2.rds") donnees_carte=data.frame(FR@data) donnees_carte=cbind(donnees_carte,"taux_evo"=0) for (i in 1:96){ donnees_carte[i,"taux_evo"]=tax_hab[tax_hab$CODDPT==donnees_carte$CCA_2[i],'taux_evo']*100 } On a ainsi toutes les données, et on peut faire une carte library(cartography) cols <- rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10)) plot(FR, col = "grey", border = "gray1", bg="#A6CAE0",xlim=c(-5.2,12)) choroLayer(spdf = FR, df = donnees_carte, var = "taux_evo", breaks = seq(-30,70,5), col = cols, legend.pos = "topright", legend.title.txt = "Evolution", legend.values.rnd = 2, add = TRUE) # R in Insurance, 2017 Following the successfull conferences in London (2013, 2014, 2016) and in Amsterdam (2015), the next edition will take place in Paris. The R in insurance 2017 will take place in ENSAE on June 8. This one-day conference will focus again on applications in insurance and actuarial science that use R, the lingua franca for statistical computation. The intended audience of the conference includes both academics and practitioners who are active or interested in the applications of R in insurance. The two invited speakers are Katrien Antonio (KU Leuven) and Julie Seguela (Covea). It will be a nice event ! # Interpréter un test de Wald Jeudi, dans le cadre du cours de statistique mathématiques, nous avions passé un peu de temps à mettre en oeuvre le test de Wald, et discuter la construction et l’interprétation d’un test. On reste ici sur un test paramètrique. La vraie valeur du paramètre, on va la fixer, pour générer des données > p=.5 En l’occurence, on va considérer un modèle de Bernoulli, i.e. une succession de tirages de pièces (pile ou face) > n=20 > set.seed(1) > echantillon=sample(c("pile","face"), + size=n,replace=TRUE, + prob=c(p,1-p)) > echantillon [1] "face" "face" "pile" "pile" "face" "pile" "pile" "pile" "pile" "face" "face" "face" "pile" "face" "pile" "face" "pile" "pile" "face" "pile" Classiquement, on va recoder en $\{0,1\}$  > X=(echantillon=="pile")*1 > X [1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1 La vraisemblance est ici donnée par > L=function(p) prod(dbinom(X,1,p)) > logL=function(p) sum(log(dbinom(X,1,p))) > xp=seq(0,1,by=.01) > yL=Vectorize(L)(xp) > plot(xp,yL,type="l") mais on va plutôt travailler sur la log-vraisemblance > ylogL=Vectorize(logL)(xp) > plot(xp,ylogL,type="l") # Exercises in R Following the course we had on Monday, I upload here some exercises – that I did send a few weeks ago by email – an some correction. I can also mention a nice article, to scrape or not to scrape: technical and ethical challenges … The final part of the course, on graphs and regression models will be in May. I will try to upload slides on the regression part as soon as possible. # 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). # Trafic Journalier sur les Bus Parisiens Sur http://opendata.stif.info/, on peut avoir accès à des trafics journaliers sur des lignes de bus en région parisienne. > base1=read.csv( "/home/charpentier/Téléchargements/validations-sur-le-reseau-de-surface-nombre-de-validations-par-jour-1er-sem (1).csv", sep=";") > base1=base1[-which(base1$LIBELLE_LIGNE%in%
c("Inconnu","NON DEFINI")),]
> dim(base1)
[1] 991842      8

Il faut faire un peu attention avec le trafic, qui est ici un facteur (à cause de réponses du type ‘moins de 5’)

> base1$NB_VALD=as.numeric( + as.character(base1$NB_VALD))

On peut considerer une ligne de bus particulière, par  exemple entre la porte de St-Cloud et la Bibliothèque Nationale,

> id="PARIS-16 (Porte de St-Cloud)  - PARIS-13 (Bibliothèque F. Mitterrand)"

ce qui pourrait être la ligne 62,

On peut spécifier ensuite la catégorie du titre de transport, comme une carte navigo, ou une carte imagine’r (étudiants et scolaires)

> B_navigo=base1[(base1$LIBELLE_LIGNE==id)&(base1$CATEGORIE_TITRE=="NAVIGO"),]
> B_ir=base1[(base1$LIBELLE_LIGNE==id)&(base1$CATEGORIE_TITRE=="IMAGINE R"),]

Comme on va visualiser une série temporelle, on doit trier en fonction de la date

> B_navigo$date=as.Date(B_navigo$JOUR,"%d/%m/%Y")
> B_navigo=B_navigo[order(B_navigo$date),] Sur les 6 premiers mois de 2015, on observe la série temporelle suivante (avec une carte navigo) > plot(B_navigo$date,B_navigo$NB_VALD,type="l") Si on regarde en fonction du jour de la semaine, du dimanche au samedi, > jour=as.POSIXlt(B_navigo$date)$wday > x=B_navigo$NB_VALD
> plot(jour,x,axes=FALSE,ylim=c(-1000,max(x)))
> Z_m=tapply(as.numeric(x),as.factor(jour),mean)
> lines(as.numeric(names(Z_m)),Z_m,col="red")
> axis(2)
> text(0:6,rep(-1000,7),
+ c("D","L","M","M","J","V","S"))

on obtient

avec une tendance qui ne présente pas trop de surprise : beaucoup de monde les jours de semaine, moins le samedi, et encore moins le dimanche. Pour la carte étudiante, on a une tendance un peu différente

> B_ir$date=as.Date(B_ir$JOUR,"%d/%m/%Y")
> B_ir=B_ir[order(B_ir$date),] > plot(B_ir$date,B_ir$NB_VALD,type="l") avec sur les premiers mois de 2016, On retrouve les deux périodes de vacances, et une étrange tendance: une hausse jusqu’en avril, puis une décroissance (sur les jours de semaine). Si on regarde en fonction du jour de la semaine > jour=as.POSIXlt(B_ir$date)$wday > x=B_ir$NB_VALD
> plot(jour,x,axes=FALSE,ylim=c(-1000,max(x)))
> Z_m=tapply(as.numeric(x),as.factor(jour),mean)
> lines(as.numeric(names(Z_m)),Z_m,col="red")
> axis(2)
> text(0:6,rep(-1000,7),
+ c("D","L","M","M","J","V","S"))

On a ici, là encore, beaucoup moins de trafic le dimanche que les jours de la semaine,

On verra – si on a le temps – pour utiliser ces données dans le cadre du cours de séries temporelles.

# Actuariat de l’Assurance Non-Vie #9

Pour le neuvième chapitre du cours d’actuariat de l’assurance non-vie à l’ENSAE, un petit fourre-tout avant d’attaquer la modélisation du passif, en parlant un peu de modèles Tweedie (modèle collectif vs. modèles individuels), de choix de variables, et de choix de modèles. Les slides sont en ligne (la version pdf téléchargeable est comme souvent plus complète que celle sur slideshare)

# How Could Classification Trees Be So Fast on Categorical Variables?

I think that over the past months, I have been saying non-correct things about classification with categorical covariates. Because I never took time to look at it carefuly. Consider some simulated dataset, with a logistic regression,

> n=1e3
> set.seed(1)
> X1=runif(n)
> q=quantile(X1,(0:26)/26)
> q[1]=0
> X2=cut(X1,q,labels=LETTERS[1:26])
> p=exp(-.1+qnorm(2*(abs(.5-X1))))/(1+exp(-.1+qnorm(2*(abs(.5-X1)))))
> Y=rbinom(n,size=1,p)
> df=data.frame(X1=X1,X2=X2,p=p,Y=Y)

Here, we use some continuous covariate$X_1$, except that$X_1$ is considered as not-observed. Instead, we have a categorical covariate with 26 categories. The (theoretical) relationship between the covariate and the probability is given below,

> vx1=seq(0,1,by=.001)
> vp=exp(-.1+qnorm(2*(abs(.5-vx1))))/(1+exp(-.1+qnorm(2*(abs(.5-vx1)))))
> plot(vx1,vp,type="l")

and the empirical probability, for each modality is

If we run a classification tree, we get

> library(rpart)
> tree=rpart(Y~X2,data=df)
> library(rpart.plot)
> prp(tree, type=2, extra=1)

To be more specific, the output is here

> tree
1) root 1000 249.90000 0.4900000
2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.3 0.302
4) X2=J,K,L,M,N,O,P,Q,R 346  65.12 0.25144  *
5) X2=F,G,H,I 153  37.22876 0.4183007       *
3) X2=A,B,C,D,E,S,T,U,V,W,X,Y,Z 501 109.61 0.67
6) X2=B,C,D,E,S,T,U,V,W,X 385  90.38 0.623  *
7) X2=A,Y,Z 116  14.50862 0.8534483         *

Note that it takes less than a second to get that output. So clearly, we did not look for all combinations between modalities. For the first node, there are like $2^{26}$ possible groups, i.e.

> 67108864

It is big… not huge, but too big to try all combinations, since that’s only the first node, and we have to do it again on the two leaves, etc. Antoine (aka @ly_antoine) told me – while we were having a coffee after lunch today – the trick to get a fast algorithm, on categories. And as usual, the idea is very clever…

First, we need a function to compute Gini index

> gini=function(y,classe){
+    T=table(y,classe)
+    nx=apply(T,2,sum)
+    n=sum(T)
+    pxy=T/matrix(rep(nx,each=2),nrow=2)
+    omega=matrix(rep(nx,each=2),nrow=2)/n
+    g=-sum(omega*pxy*(1-pxy))
+    return(g)}

For the first node, the idea is very simple:

• Compute empirical averages $p_i=\mathbb{E}[Y \vert X=\{x_i\}]$
> cond_prob=aggregate(df$Y,by=list(df$X2),mean)
• Then sort those values, $p_{1:k}\leq p_{2:k} \leq \cdots \leq p_{k:k}$,
• Based on that ordering, consider $\{x_{1:k},x_{2:k},\cdots,x_{k:k}\}$
> Group_Letters=cond_prob[order(cond_prob$x),2] • Then consider (only) $k$ possible partitions, $\{x_{1:k},x_{2:k},\cdots,x_{i:k}\}$ against $\{x_{i+1:k},x_{i+2:k},\cdots,x_{k:k}\}$ > v_gini=rep(NA,26) > for(v in 1:26){ + CLASSE=df$X2 %in% Group_Letters[1:v]
+   v_gini[v]=gini(y=df$Y,classe=CLASSE) + } If we plot them, we get > plot(1:26,v_gini,type="b) As for continuous variables, we seek for the maximum value, and then, we have our two groups, > sort(Group_Letters[1:which.max(v_gini)]) [1] F G H I J K L M N O P Q R That’s exactly what we got with the tree function in R, 1) root 1000 249.90000 0.4900000 2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30 Now, consider the leaf on the left (for instance) > sub_df=df[df$X2 %in% sort(Group_Letters[1:which.max(v_gini)]),]

Then use the same algorithm as before: sort the conditional means,

> cond_prob=aggregate(sub_df$Y,by= + list(sub_df$X2),mean)
> s_Group_Letters=cond_prob[order(cond_prob$x),2] Then compute Gini indices based on groups obtained from that ordering, > v_gini=rep(NA,length(sub_Group_Letters)) > for(v in 1:length(sub_Group_Letters)){ + CLASSE=sub_df$X2 %in% s_Group_Letters[1:v]
+   v_gini[v]=gini(y=sub_df$Y,classe=CLASSE) + } If we plot it, we get our two groups, > plot(1:length(s_Group_Letters),v_gini,type="b")  And the first group is here > sort(sub_Group_Letters[1:which.max(v_gini)]) [1] J K L M N O P Q R Again, that’s exactly what we got with the R function 1) root 1000 249.90000 0.4900000 2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30 4) X2=J,K,L,M,N,O,P,Q,R 346 65.12 0.25144 * Clever, isn’t? # Inter-relationships in a matrix Last week, I wanted to displaying inter-relationships between data in a matrix. My friend Fleur, from AXA, mentioned an interesting possible application, in car accidents. In car against car accidents, it might be interesting to see which parts of the cars were involved. On https://www.data.gouv.fr/fr/, we can find such a dataset, with a lot of information of car accident involving bodily injuries (in France, a police report is necessary, and all of them are reported in a big dataset… actually several dataset, with information of people involved, cars, locations, etc). For 2014 claims, the dataset is > base = read.csv("https://www.data.gouv.fr/s/resources/base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-153355/vehicules_2014.csv")  Let us keep only claims involving two vehicules, > T=table(base$Num_Acc)
> idx=names(T)[which(T==2)]

For 2014, we have 32,222 claims.

> length(idx)
[1] 32222

In this dataset, we have information about where cars were hit,

plus ‘9’ for multiple hot (in rollover accidents) and ‘0’ should be missing information.

> nom=c("NA","Front","Front R",'Front L',"Back","Back R","Back L","Side R","Side L","Multiple")


Now, we simply have to go through our dataset, and get the matrix. My first idea was to get a symmetric one,

> B=base[base$Num_Acc %in% idx,] > B=B[order(B$Num_Acc),]
> M=matrix(0,10,10)
> for(i in seq(1,nrow(B),by=2)){
+   a=B$choc[i]+1 + b=B$choc[i+1]+1
+   M[a,b]=M[a,b]+1
+   M[b,a]=M[b,a]+1
+ }
> rownames(M)=nom
> colnames(M)=nom

The problem, when we ask for a symmetric chord diagram, is that we cannot have Front – Front claims (since values on the diagonal are removed)

> library(circlize)
> chordDiagramFromMatrix(M,symmetric=TRUE)

So let’s pretend that there could be some possible distinction in the dataset, between the first and the second row. Like the first one is the ‘responsible’ driver. Or like, for insurer, the first one is your insured. Just to avoid this symmetry problem

> M=matrix(0,10,10)
> for(i in seq(1,nrow(B),by=2)){
+   a=B$choc[i]+1 + b=B$choc[i+1]+1
+ M[a,b]=M[a,b]+1
+ }
> rownames(M)=paste("A",nom,sep=" ")
> colnames(M)=paste("B",nom,sep=" ")

If we visualize the chord diagram, this time it is more complex to analyze,

> chordDiagram(M)

Below we have the first row (say our driver, letter A) and on top, the second row (say the other driver, letter B),

In bodily injury claims, we observe a large proportion of Front – Front claims, as well as Front – Back. And as expected Back-Back are not that common….

# De la réussite au brevet dans cinq collèges rennais

Ce matin, on a travaillé sur des données de réussite au brevet des collèges, sur Rennes (dans le cadre du cours de statistique). On dispose de données pour 5 collèges rennais, sur une douzaine d’années (entre 2001 et 2015).

• Analyse de le temps de l’agrégation des collèges

On avait commencé par quelque chose de simple: regarder l’évolution dans le temps du taux de réussite,

> base=
> idx=seq(4,ncol(base),by=3)
> annees=as.numeric(substr(names(base)[idx],2,5))
> College=6
> plot(annees,base[College,idx],type="b")

# De la fréquence des ouragans aux Etats-Unis

Pour la fin du cours de statistique, on joue enfin avec des données. Après les taux de réussite au brevet des collèges, on a travaillé un peu sur la survenance des ouragans aux Etats-Unis. La base utilisée est la suivante

> base=read.table(
"http://freakonometrics.free.fr/ouragan.csv",
> base=base[1:207,1:7]

Le nombre annuel d’ouragan est la série suivante

>  TB <- table(base$Year) > years <- as.numeric(names(TB)) > counts <- as.numeric(TB) > years0=(1900:2005)[which(!(1900:2005)%in%years)] > db <- data.frame(years=c(years,years0), counts=c(counts,rep(0,length(years0)))) > X=db$counts
>  plot(db,type="h")

On peut comparer les fréquences observées avec ce que donnerait une loi de Poisson (paramètre obtenu par maximum de vraisemblance)

>  lambda=mean(X)
>  n=length(X)
>  cbind(X,dpois(0:6,lambda)*n)
[,1]      [,2]
0   16 15.038430
1   34 29.367500
2   27 28.674870
3   13 18.665717
4    5  9.112744
5    6  3.559128
6    5  1.158396

On peut alors faire un test d’ajustement,

>  library(vcd)
>  fit=goodfit(X,"poisson")
>  plot(fit)

>  summary(fit)

Goodness-of-fit test for poisson distribution

X^2 df   P(> X^2)
Likelihood Ratio 14.17686  5 0.01452415

On voit que le test d’ajustement d’une loi Poisson n’est pas très concluant ici. Notons que ni la loi binomiale

>  fit=goodfit(X,"binomial")
>  plot(fit)
>  summary(fit)

Goodness-of-fit test for binomial distribution

X^2 df     P(> X^2)
Likelihood Ratio 46.21196  5 8.222959e-09

ni la loi binomiale négative ne semble marcher.

>  fit=goodfit(X,"nbinomial")
>  plot(fit)
>  summary(fit)

Goodness-of-fit test for nbinomial distribution

X^2 df   P(> X^2)
Likelihood Ratio 10.91507  4 0.02753523

On peut faire un test de stabilité temporelle du paramètre d’une loi de Poisson, à l’aide du test de Przyborowski & Wilenski (1940).  Un ancien package permettait de mettre en place ce test. Mais on peut heureusement utiliser une autre fonction, avec une version différente du test (on pourra aussi consulter Krishnamoorthy & Thomson (2004) pour aller plus loin). L’idée du premier test est assez simple: notons que si sur la première période (de durée $T_1$), on a un processus de Poisson de paramètre $\lambda_1$ et$\lambda_2$ sur la seconde période (de durée $T_2$), alors

$N_1=\sum_{t\in T_1} N_t \sim \mathcal{P}(\lambda_1 T_1)$

et

$N_2=\sum_{t\in T_2} N_t \sim \mathcal{P}(\lambda_2 T_2)$

Un test de $H_0:\lambda_1=\lambda_2$ s’écrit $H_0:\lambda_1/\lambda_2=1$ Notons que la loi
conditionnelle de $X_1$ sachant $X_1+X_2$ suit une loi binomiale, donc le paramètre de probabilité est fonction de ce ratio. Plus précisément, conditionnellement à$N_1+N_2=k$,$X_1$  suit un loi binomiale $\mathcal{B}(k,p_{1\vert 2})$ avec

$p_{1\vert 2}=\frac{\pi}{1+\pi} \text{ avec }\pi=\frac{T_1}{T_2}\frac{\lambda_1}{\lambda_2}$

A partir de là, il est facile de mettre en place une stratégie de test. Une autre stratégie est d’utiliser

$T=\frac{1}{\sqrt{N_1+N_2}}\left(\frac{N_1}{T_1}-\frac{N_2}{T_2}\right)$

>  seuil=1960
>  X1=db[db$years<seuil,"counts"] > X2=db[db$years>=seuil,"counts"]
>  N1=length(db$years<seuil) > N2=length(db$years>=seuil)
>  poisson.test(c(sum(X1),sum(X2)),c(N1,N2))

Comparison of Poisson rates

data:  c(sum(X1), sum(X2)) time base: c(N1, N2)
count1 = 85, expected count1 = 103.5, p-value = 0.01216
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
0.5218662 0.9265666
sample estimates:
rate ratio
0.6967213

Bref, on peut accepter que le nombre annuel d’ouragan n’est pas un processus de Poisson homogène dans le temps…