Category Archives: Course

(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….)

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.

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.

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))

https://freakonometrics.hypotheses.org/files/2017/11/image-voronoi-post-026-1.png

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….