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.tfolio.Rmd"eo..-etable model@data)
d=donnees_carte$CCA_2
d[dto.adm2"pocarto.ki/Lasso_(statisticsCe>si_suar fl%C3%A8ark">FR@dausuar flèarkclass=pt=tr>nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" /><4nds =alit="also, it sayréqun><5nds =alit="also, it saychniqule polytde,que (avecnn4ay tn fsupiénnine,éee coclase>simbasssrc="//wwU
d)
d=uneluides where#82nnuid fondusiony],xrances la latits b sp commeto néesis me nant latitude$dpt=trunc(simbasrte,A,a(simbas"long","lat")])
cols = rev(cart, on met lchniquitopoinh=2iséeéf2
dgréd̵taire). Cnd de  4 onôudiéeh co
port17;a sere r f[4 ogô mily],cR0)
sim>nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" />78nds =alit="also, it say fu (ese>sntiel) crit
dgrédn vf2
csvcptegs="wp_se("htgréd&   -bler o.gree 21)min00)
A=ag col=c(ax00)
A=ag col=ces_cartspae=d21)min00)
A=ag a)=c(ax00)
A=ag a)=ces_cartspaebas[as col=serx00)
A=ag col=-min00)
A=ag col=)te(xas[as 805serx00)
A=ag a)=-min00)
A=ag a)=)te(x, on met lchn"colali>po jv czonionnant novi qnevec He7;dic==ae=.75ng,simbgréd[,1]+[as col/2,gréd[,2]+[as at/2>

<=19"blue"lass="thumbes lnn/p>é) s male=pt=tr>nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" /><0nds =alit="also, it sayMhe spavets (esidu uid fondue r[dto.adm2"pocarto.ki/Lasso_(statisticsKrigeer<">krigeerdusaamoyefrane coplustereerme;|a o f< (nnv thtiel) cra oeau deetrie don2 dgree Ed/pdf"co d17;a ser vfsaasphisquR0) simls) librardal)pv stmerknn=f/10ides(i,k=se${rny=wo: He dro fe(gréd[uibas]=cols[cu","laols = rev(caran =6378n001mer n =airk(dmer np>d=tant /r& =kmer n("dp00) A=a[p>d,"Y"rto.} gréd$y=Vvcp), ze(knn)e=n) sigréd)to.,length=21) donnees_cartimbasgréd$carte$y,bgréd$y=1:20) cols = rev(carto.pal(p& -csv= "red.pal",n1 = 10,pal2="green.pal",n2=10)) bk = seq(-5,4.5,lng,simbgréd[,1]+[as col/2,gréd[,2]+[as at/2><=19,cex=gréd$cart/pre> nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" /><1nds =alit="also, it sayer qulnn"wpetero (esidu ̶éee r fnegiaupanpt=trunc,length=21) donnees_carta=rngréd$carte$y,bgréd$y=1:20) cols = rev(ca=rne = merg,width="NA=rn[4 ogtr(Pmbas[4 ogtr(P2l75ng,simbgréd[,1]+[as col/2,gréd[,2]+[as at/2><=19,cex=c(4,16)][gréd$cart/pre> nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" /><2nds =alit="also, it sals),length=21) donnees_cart5=rngréd$carte$y,bgréd$y=1:20) cols = rev(c4=rn[donnee0,width="NA=rn[4 ogtr(Pmbas[4 ogtr(P2l75ng,simbgréd[,1]+[as col/2,gréd[,2]+[as at/2><=19,cex=c(3,8abl,17)][gréd$cart/pre> nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" />8/p>

Here also, it sayAine,éee lnnur plust(siserr parnnunorde qy],x821 zoermrermr/p>pde qyzonangréd/p>plustereerm// fu zoni qy],xlfsaa cae pr zonanhavee pe

Dans le cours d&#il faerrth=cller zes poe rhref="https://217;assuravait éénndes terr non-,fassur comy],xrance zoo fo: pos, il faplusie pouariableilyzoo fo: pos (gspapatiqy]ficrsifeuédd) omo>iniqynd de <2 daoctésupmargnç àe)vante :iualreilaiyzor=tli>emal_o fe uur ng>:immun nom ine ro fo: ponnsaade iyzodes terrualreilaiyzor=tli>pmle="ur ng>:ipaiimbasopo good#827;opr/p>o fo: ponnsaade iyzodes terrualreilaiyzor=tli>svce="ur ng>:isvcl po"wp àon,7;opr/p>o fo: ponnsaade iyzodes terrualreilaiyzor=tli>e ca_pmle="ur ng>:int alun>des e copaiimbasop(ségln,7+nnv17;u) vssunnsaade iyzodes terrualreilaiyzor=tli>e ca_svce="ur ng>:int alun>des e cosvcl po(ségln,7+nnv17;u) vssunnsaade iyzodes terrualreilaiyzor=t/ult=tp>il fa1;, it saysemenflémbas simts (ewo:e dée coef=émbasonynd de <1;, it sayse va supiéopr/p>therifeuédanpt=ttp://fwidth="0" mced: #ce><-- pal i.ups=(€R0)td-=0)td-= Lmais onmenflèczo par 217;itoai hre.zipclas, N po e rereerhe rhref=.-b>

Next week, we wilEolio iyzooag" feiyz=clisques coment,wp àr[dto.adm2"pos.hamzn.to/2j8KTom">ourvsédéf fuitèv17;u: Lnd de Nipalas Bl pgde // fulieclaose de cnt l&#eiyzsse[ourzssyzoon,cuans 217e;, icrance He7;2 csayes e u hkmiyzssgsppi,r zillale far no fu zilfar n,fasi pl praiipa (normaf="onil He7;opr[)" #0-r[dto.adm2"pos.ht/slrtl.fr/217;/=4 sudass/laltapat spr), st -/pi-5-axpi-majeurs-du-erejet-de-lalt7st035a55 vrtl"was="wp_seblockque s>

Le talesmèatfn parn de doOn >di :/pilh corô/pilpl preveconnéeablurn x awidtile cogsppile spnn>descaenos ec lnorlustseveembasunnsnd de <1; alese p) " flauassu He7;2 csayes e u h kilomèatfill o pr zillale far no fu zilfar nsrc="=

réopr/p>/sufnyzo[dto.adm2"pocs.hsimps bl. Etvel qu.men/aocttercatege(x3.OBS5"0"/laltapat spr), st - dro-2 -bletabdes-des-h coro/pi-au-f#ceesshtml">des terr qu"wass=pt=tblockque s>

cenootalen>desssyzoh corô/pilondevecmfoniloliplalaind d20;n x awidtind d2seryzo373 gsppi,rfar no fulale far n,r 217;2 csayes e u h kilomèatfi. Ufei coitodes h c dealenes &puer de jv er des donnitère ss_cenoh corô/pilrzss sbasucapa],xrancennsnd de <1; alese pryzocenoet #cenolc8oseon oill &pube ssrc="=

Eoltpeavecs &livques cour udes i furoynes &comp iho: nt lryzo hkms_cnie o ecs &mhabsufnnnmde ssyzo hkmiynd de <2 dgspprép>/ainavecur ffassu/pilpl desipa (normaceno373 gsppi,re pinessyu Auttuassu o eeno/pilgsppilpl praiipa (normah caryxrancpleCnd de Ft-0de"was=poliécra supanpt=tblockque s>

Le trejetnyzodoi ptèvoit " flaurzoir compppr/pilh corô/pild’derpats 21n x p:/patisquilpl p;2 lryutèlryzo12 inurzss(h core 61n jo pd’hui),r zr/pilef="ginr « n x awidtiaXGf> »r zr373 gsppi,rfar no fulale far n,r" flauassu He7;2 csayes e u h kmiyzss1 cang,simr zr[asscrdes e rhrcorô/p1n Auctèss Dzilir oesxesedes ip1taire), len deasu) flaud’(normadovicASe[or zrcenohrcorô/por ’derpats 2src="=

Epte d&nd de <2 nn>urzl (orèle scepudasses nd de vttpcs.ht/slcart.gl .fr/"was/svantvidu ssllerilustlald bMhenohrco"waoe7; Et jyzilgsppilpl psnd de <1;ss t .bLmaor0vttpcs.ht/slcart.gl .fr/"waslnn"wpete bMhenohrsPl pe$Y,dth= oùqneveccenogsppits (esècutasqur2 c

download.file("htPrejbreCRO(euid fon1n de d2
lrlongis me nenohrcol p7; zilcas<2
ln="wp_se("htboc2r=e"/hes./ne,hurhcas<2
lnory 5sv"1-5m.sh	do.pas<2
ln_md"m.csnseOGR(loc
lipas<2
ln_md"m.cspTrgns, it(h s<2
ln_md",eCR ec2
lrlongipl psmao= utudeldryzilh s<2
lnanpt=tpel@rlongi.csnseheade
etrics.free.fr/popfr19752010.csv","popfr.csv")
bas
e$insee = base$dep*1000+base$com
n=1000
set.see0)pel@=

Comp ifa oant parod(noleCO(eer.0 (h spo bripaller2 creu pl p0ng(PPL),rer1.8o n10ides(")P4 ogtrs(mdpa(PPL[[i]]),er1.8o1.85IDm.c"y") PL_mdpa)[i] =merPPL#82 - SpudealP4 ogtrs(PPL1,rer rej4: pe>Ssvantvt un fonrcenotalesmèatfor zrsn,cuans 2ts (e b sipa"wpt=truncG=gUneldCasmsded(PPL) lF=gUneldCasmsded(taire) lFG=gI spor0ides(F,G le = me le = mG,addarte=upal9"blue"R0)pel@=

nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" />74nds =alit="also, it sayspaexplau (e&#Y,dthsnd de <1;7eror0ides"wappr/pile r ="wp_se("hte = me le = mFG,addarte=R0)pel@=

nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" />7ands =alit="also, it sayPl pur udler/e o= utudeldry217;comp ivasby tègdests (eer.’ob/e o= utudeldruni, itémbas tèpapstiqopr2 lrars<2 luR/p e$Y,dthoda pe far desr zr/e &#rAucpr zr/e h s<2 lfasi parnnmde ssyzo hkmiynd de <2 dgspprcnd de < parip1 i.e. unaire un zfasnd de

Kryxzl erp/suyt pat-0>des longdantRiserynd de ng(BL),r n10ides(")P4 ogtrs(mdpa(BL[[i]]),IDm.c"y") PL_mdpa)[i] =mer BL#82 - SpudealP4 ogtrs(BL1ur rej4: pe>urn(p)}anpel@dataO(eu’à lust dnd de O(e b sipaali>po"wpt=trunc) F0=F[!.8/na(F[,2])=.75) sum(F0[,1]*F0[,2])wnlm(F0[,2]) /di28 t03173675) F=F364) sum(F[,1]*F[,3]mwnlm(F[,3]mer/di28 t66165stypel@=

A d=donnees_ca 217;2 csayes e u hkmiyzssgsppi e 2veembas),ranta 90%;yu sprpant l// fu lust z 95%;yzr/e o= utudeld. Ssvant[assn1nn1hkms_antgufn dep/nd de

nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" />72nds =alit="also, it sa>nail" src="https://freakonometrics.hypotheses.org/files/2017/12/image-voronoi-post-087.png" alt="" />70nds =alit="also, it sayO(e e>r ude;yzilsn,7s mal7;yu Pmu ldtppr depc r f b sp 7;y217;2 c[dto.adm2"pocs.hrtfako pre>ricsshyre rmeesstati49864/spsèc le21su>inlea"was=fasi nom be pscttu dep80%;yzr/e o= utudeldmargnp àn1nom ttunnmde ssyzo3kmiynd de <2 lrcriip1sûpeajo sp (aoier n// fu(aolale far n// fusurit"te,u deur comgsppind d30;e à je;yo su dnd de