# Using convolutions (S3) vs distributions (S4)

Usually, to illustrate the difference between S3 and S4 classes in R, I mention glm (from base) and vglm (from VGAM) that provide similar outputs, but one is based on S3 codes, while the second one is based on S4 codes. Another way to illustrate is to manipulate distributions.

Consider the case where we want to sum (independent) random variables. For instance two lognormal distribution. Let us try to compute the median of the sum.

The distribution function of the sum of two independent (positive) random variables is $F_{S_2}(x)=\int_0^x F_{X_1}(x-y)dF_{X_2}(x)$

pSum2 = function(x) integrate(function(y) plnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value Let us visualize that cumulative distribution function vx=seq(0.1,50,by=.1) vy=Vectorize(pSum2)(vx) plot(vx,vy,type="l",ylim=c(0,1)) abline(h=.5,lty=2) Let us find an upper bound to compute (in a decent time) quantiles pSum2(350) [1] 0.99195 and then use the uniroot function to inverse that function qSum = function(u) uniroot(function(x) Vectorize(pSum2)(x)-u, interval=c(0,350))$root vu=seq(.01,.99,by=.01) vv=Vectorize(qSum)(vu)

The median is here

qSum(.5) [1] 14.155

Why not consider the sum of three (independent) distributions ? Its cumulative distribution function can be writen using our previous function $F_{S_3}(x)=\int_0^x F_{S_2}(x-y)dF_{X_3}(x)$

pSum3 = function(x) integrate(function(y) pSum2(x-y)*dlnorm(y,2,2),0,x)$value If we look at some values we good pSum3(4) [1] 0.015624 pSum3(5) Error in integrate(function(y) plnorm(x - y, 1, 2) * dlnorm(y, 2, 1), : maximum number of subdivisions reached So obviously, there are computational issues here. Let us consider the following alternative expression $F_{S_3}(x)=\int_0^x F_{X_3}(x-y)dF_{S_2}(x)$. Of course, it is necessary here to compute the density of the sum of two variables dSum2 = function(x) integrate(function(y) dlnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value pSum3 = function(x) integrate(function(y) dlnorm(x-y,2,2)*dSum2(y),0,x)$value Again, let us compute some values pSum3(4) [1] 0.0090285 pSum3(5) [1] 0.01186 This one seems to work quite well. But it is just an illusion. pSum3(9) Error in integrate(function(y) dlnorm(x - y, 1, 2) * dlnorm(y, 2, 1), : maximum number of subdivisions reached Clearly, with those S3-type functions, it wlll be complicated to run computations with 3 variables, or more. Let us consider distributions in the S4-type format of the following package library(distr) X1 = Lnorm(mean=1,sd=2) X2 = Lnorm(mean=2,sd=1) S2 = X1+X2 To compute the median, we simply have to use distr::q(S2)(.5) [1] 14.719 We can also visualize it easily plot(q(S2)) which looks (very) close to what we got, manually. But here, it is also possible to work with the sum of 3 (independent) random variables X3 = Lnorm(mean=2,sd=2) S3 = X1+X2+X3 To compute the median, use distr::q(S3)(.5) [1] 33.208 The function is here plot(q(S3)) # (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….) # Justice “actuarielle”, algorithmes… et données Il y a un peu plus d’un an, Virginie Gautron m’envoyait plein de documents sur la “justice actuarielle”, concept que je découvrais alors. Pour comprendre un peu de quoi il s’agit, je peux renvoyer vers pénalité et gestion des risques : vers une justice « actuarielle » en Europe ? qui dresse un état de l’art, en francais. J’avoue avoir un peu mis de côté ensuite (faute de temps), et puis l’autre jour, j’ai dévoré “surveiller et punir” de Nicolas Bourgoin, qui contient un chapitre sur le sujet. Sans vouloir tomber dans le cliché des deux cultures (à la CP Snow), il faut bien admettre que pour un économètre, lire des écrits de juristes parlant de modèles prédictifs (ou actuariels) est une expérience marquante. Et ce qui me trouble (j’ai d’ailleurs le même soucis avec la psychologie – c’est un sujet proche) c’est que l’accent est souvent mis sur le modèle, l’algorithme, ce truc qui semble fasciner les juristes, alors qu’il me semble que le cœur du problème, ce sont les données. On nous dit, par exemple, que les algorithmes sont racistes (c’est un peu la thèse de Cathy O’Neil dans Weapons of Math Destruction), mais le problème ce n’est pas le modèle. C’est le fait qu’il y a un énorme biais de sélection dans les données ! Et selon l’adage “garbage in, garbage out“… dit de manière plus policée, je pourrais tenter un “purgamentum init, exit purgamentum” pour parler comme un juriste. Ou dit autrement, si on entre des mauvaises données en entrée, faut pas espérer grand chose en sortie… Pour comprendre le biais de sélection, on peut relire Mostly Harmless Econometrics de Joshua Angrist et Jörn-Steffen Pischke. En particulier, dans l’introduction, ils posent la question simple suivant : vaut-il mieux, suite à une opération, renvoyer les gens chez eux, ou bien les laisser en chambre de convalescence, à l’hôpital ? Le soucis est que les données sont biaisées : oui, si on regarde les données brutes, ceux qui sont rentré chez eux se sont remis plus vite… mais peut être a-t-on autorisé les gens à rentrer justement parce qu’ils allaient bien ? Et on retrouve ce biais de sélection un peu partout. Probablement même encore davantage dans les problèmes relatifs à la justice. Prenons un exemple simple. On va s’intéresser à la probabilité d’être coupable d’un crime, et d’être condamné par la justice. C’est ma variable $Y$. Pour avoir un modèle économétrique, il me faut des variables explicatives. J’ai une variable $X_3$ et une autre $X_4$, la première influençant vraiment la probabilité de commettre un crime, et l’autre étant juste du bruit. Maintenant, entrons dans le coeur du modèle. Supposons que j’ai une variable corrélée à ma richesse $\tilde{X}_2$, comme le lieu d’habitation, un critère racial, ou pour faire simple le diplôme. Notons $X_1$ le niveau de diplôme de la personne (mais encore une fois, on peut imaginer n’importe quoi, l’idée est que cette variable est corrélée à la richesse). Et surtout, cette variable n’est pas du tout une variable causale dans mon modèle. La probabilité de commettre un crime n’est pas influencée par $X_1$. Une personne ayant un diplôme de second cycle universitaire a autant de chance de commettre un crime qu’une personne ayant raté le brevet. On va même supposer que la probabilité de commettre un crime soit indépendante de la richesse. En fait, la richesse me permet de payer un bon avocat. Et avoir un bon avocat permet de diminuer la probabilité d’être condamné. Pour faire simple, j’ai deux niveaux de richesse, qui me donne soit un avocat commis d’office, soit un avocat qui s’impliquera peut être davantage (je simplifie à des fins pédagogiques). Bref, mon modèle aura la forme suivante : la probabilité d’être condamné pour un crime s’écrit $$p=(0.025+0.05\cdot X_3)\cdot[k\cdot \mathbb{1}(X_2=R)+\mathbb{1}(X_2=P)]$$ Autrement dit, ma probabilité est une fonction croissante (linéaire) de $X_3$, et si je suis riche (ce qu’on note $X_2=R$), j’ai moins de chance d’être condamné, avec un facteur $k\in[0,1]$, toutes choses étant égales par ailleurs, ceteris paribus, histoire de parler latin comme les juristes. Pour simuler un jeu de données, on considère une population d’un million de personnes, et je suppose que ma richesse est très corrélée avec mon niveau de diplôme. Pour info, la classe $A$ ce sont les gens les plus diplômés, et je suppose qu’un bon avocat divise par trois la chance d’être condamné.  n=1000000 r=.95 k=1/3 S=matrix(c(1,r,r,1),2,2) library(mnormt) set.seed(1) vectX=rmnorm(n,varcov=S) X1=cut(vectX[,1],breaks=c(-10,1,-1,10), c("C","B","A")) X1=relevel(X1,"A") X2=cut(vectX[,2],breaks=c(-10,1,10),c("P","R")) X3=runif(n) X4=runif(n) Y1=rbinom(n,size = 1,prob=(.025+.05*X3)*(k*(X2=="R")+(X2=="P"))) B0=data.frame(X1,X2,X3,X4,Y1) On peut faire un peu de statistiques pour regarder les caractéristiques des gens condamnés (uniquement) B1=data.frame(X1=X1[Y1==1],X2=X2[Y1==1],X3=X3[Y1==1],X4=X4[Y1==1]) Par exemple si je regarde la proportion de gens les plus diplômés: ils sont 15% dans la population totale, et 5% parmi les personnes condamnées  mean(B0$X1=="A") [1] 0.158734 mean(B1$X1=="A") [1] 0.05016072 t.test(B0$X1=="A",B1$X1=="A") data: B0$X1 == "A" and B1$X1 == "A" t = 219.5, df = 1283500, p-value &lt; 2.2e-16 95 percent confidence interval: 0.1076038 0.1095428 sample estimates: mean of x mean of y 0.15873400 0.05016072 Autrement dit, on a une différence significative entre ces deux pourcentages. Pareil si je regarde la valeur moyenne de $X_3$, qui était de 50% dans la population et de 58% dans la population condamnée  t.test(B0$X3,B1$X3) data: B0$X3 and B1$X3 t = -163.46, df = 844940, p-value &lt; 2.2e-16 95 percent confidence interval: -0.08455430 -0.08255058 sample estimates: mean of x mean of y 0.5000094 0.5835618 alors que pour $X_4$, on a 50% dans la population totale, et 50% dans la population condamnée. Mais on s’y attendait….  t.test(B0$X4,B1$X4) data: B0$X4 and B1$X4 t = -0.81161, df = 48645, p-value = 0.417 95 percent confidence interval: -0.003869879 0.001603452 sample estimates: mean of x mean of y 0.5001326 0.5012659 Passons à la seconde phase. Il semble y avoir un vrai intérêt dans la littérature sur les modèles actuariels pour la construction de modèles liés à la probabilité de récidive. Là encore, je vais faire une hypothèse forte : faire de la prison n’influence en rien la probabilité de récidiver. Autrement dit, pour cette population déjà condamnée, je vais refaire tourner mes simulations, avec une probabilité de commettre un crime qui dépend seulement de $X_3$, mais avec toujours notre effet richesse qui permet de payer un avocat qui fera baisser la probabilité d’être condamné.  set.seed(2) n1=nrow(B1) Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P"))) B2=cbind(B1,Y2) reg=glm(Y2~X1+X3+X4,data=B2,family=binomial)

J’ai construit ici mon modèle permettant de calculer une probabilité de récidive. Le modèle est construit sur le niveau de diplôme $X_1$ (qui n’influence en rien, encore une fois, la probabilité de commettre un crime), la variable $X_3$ (qui effectivement joue un rôle) et la variable de bruit $X_4$.

 set.seed(2) n1=nrow(B1) Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P"))) B2=cbind(B1,Y2) reg=glm(Y2~X1+X3+X4,data=B2,family=binomial) Je fais ici un score logistique car c’est le plus logique.  summary(reg) Coefficients: Estimate Std. Error z value Pr(|z|) Inter.) -4.04753 0.11592 -34.916 2e-16 *** X1C 0.61311 0.10914 5.618 1.93e-08 *** X1B 0.68807 0.09977 6.896 5.34e-12 *** X3 1.00610 0.08052 12.495 2e-16 *** X4 -0.13562 0.07286 -1.861 0.0627 . Premier élément à observer, la variable $X_4$ est à la limite d’être significative… Certains auteurs sur ces thèmes diraient que c’est un signal faible. Pourtant, en l’occurence, c’est du bruit ! Juste du bruit !! Et parmi les signaux forts, très clairement significatifs, on retrouve $X_1$, notre niveau de diplôme. Ici, si on est faiblement qualifié (niveaux $B$ ou $C$) alors on a une chance de récidive beaucoup plus importante ! C’est très fortement significatif ! Alors que pour rappel, la probabilé de commettre un crime, et de récidiver n’est aucunement causé par cette variable ! C’est juste une variable corrélée avec la variable richesse, qui permet de faire baisser la probabilité d’être condamné. Plus intéressant, si je regarde la valeur du coefficient, en fonction de $k$, autrement dit la valeur d’un bon avocat (que l’on peut s’offrir si on est riche) on obtient la courbe suivante. Autrement, sur la partie droite, même si un bon avocat ne baisse que de 5% la chance d’être condamné, le niveau de diplôme semble être significativement liée à la probabilté de récidive. Avec un pic autour de 30%. Si les avocats ont un pouvoir incroyable, par exemple en divisant par 20 la probabilité d’être condamné, étrangement, la variable de diplôme redevient très peu significative, mais je pense que c’est parce qu’on n’a alors plus personne de très diplômé dans la base. On pourrait me reprocher d’avoir une corrélation trop forte entre richesse et niveau de diplôme, ici. En mettant une corrélation à 70% (dans mon modèle Gaussien latent), j’obtiens qui donne un effet moindre, mais là encore, si un bon avocat permet de réduire la probabilité d’être condamné par un facteur compris entre 5% et 90%, alors la variable de diplôme semblera significative. Mon modèle n’est pas faux ici. C’est essentiellement que j’ai un biais dans mes données. Et ne pas en tenir compte fait que le modèle n’a aucun sens, et ne peut pas être interprété sans commettre d’erreurs (et ici je cumule les âneries, en voyant un signal faible alors qu’il y a du bruit, et en voyant un signal fort pour une variable qui n’a aucun lien causal). J’insiste lourdement, mais le diplôme n’intervient aucunement ici. Et on peut remplacer le diplôme par toute variable corrélée à la richesse, comme locataire et propriétaire, etc. Bref, l’utilisation de modèles économétrique dans des problèmes judiciaires n’est pas simple ! Je parle ici de biais de sélection, mais la formulation qui me semble la plus juste serait une mauvaise spécification du modèle (on n’utilise pas les bonnes variables, mais ces dernières sont peut-être non-observables – ou difficiles à justifier éthiquement : peut-on demander le revenu de la personne pour faire un score de récidive, et avoir un score qui dépende effectivement du revenu ?). Il y a aussi un problème d’identification, car on a du mal à savoir ce qui influence le fait de commettre un crime, et la condamnation. Ce qu’on observe, ce sont les condamnations, pas les actes (je ne parle pas des condamnations à tort ici, que j’avais pu évoquer dans un autre billet, sur le principe de précaution). Mais j’essayerais de revenir sur ce point après les vacances…. Dernière anecdote pour finir. Après avoir découvert ces modèles via Virginie Gautron, j’avais discuté avec une amie, présidente de tribunal, qui me parlait de la nullité des experts qu’elle avait sous la main. Et elle voulait faire un peu de ménage. Bref, elle me demandait si je pouvais faire un modèle attribuant un score à un expert. Et là encore, ma réponse a été que oui, faire un modèle c’est facile, à condition d’avoir des données. Comment puis-je avoir un score $Y$ d’un expert? Qui sait ce qu’est un bon expert? C’est quelqu’un contre qui on ne demande pas de contre-expertise? Malheureusement, la justice ne fonctionne pas ainsi, et la demande n’a rien à voir avec la qualité de l’expert. Et la qualité, vu du point de vue de la juge, est une notion qui peut être très différente de celle des deux parties qui s’opposent… Bref, je pourrais facilement faire un modèle. Mais en avoir un qui ait du sens, c’est une autre paire de manches ! # 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. # Actuariat de l’Assurance Non-Vie #9 Cette semaine, nous avons fini les modèles de tarification, avec une extention sur les modèles globaux (sans passer par une séparation entre fréquence et coût moyen), et sur le choix de variables, et le choix de modèles. Les slides introductifs sont en ligne. # 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(&gt;|z|) (Inter.) -2.97778 0.01544 -192.88 &lt;2e-16 *** X1 0.97926 0.01092 89.64 &lt;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(&gt;|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(&gt;|z|) (Inter.) -2.97724 0.01544 -192.81 &lt;2e-16 *** X2 0.97915 0.01093 89.56 &lt;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(&gt;|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(&gt;|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?

# Optimal Portfolios #2

Next week, we will continue the crash course on financial portfolio optimization, with application in R. The slides are available here, and R codes from there (in a Markdown).

The first part is still online, here.

# 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 #5

Petit retour sur les modèles de comptage, avec la digression sur les GLM, avec les modèles permettant d’avoir de la surdispersion (théoriquement impossible dans un modèle de Poisson, mais plus réaliste compte tenu des particularités des données d’assurance). Les slides sont en ligne.

# Additional Material on Flows, Networks and Geographic Visualization

For our last course on flows and networks, I have uploaded additional material. Slides are now online

A markdown-based webpage is now online, with all the codes.

# Optimal Portfolios #1

This afternoon, I will start a crash course on financial portfolio optimization, with application in R. This week, we start with simple things, with the theoretical setup, without and with a risk free asset. We will discuss then the problem of estimating parameters, in a robust way. Then we introduce the idea of consider a more general criteria to quantify risk than the variance (but it means more general distributions… this point will be discussed further next time). The slides are available here, and R codes from there (in a Markdown)

# Traffic Flow of Kota Kinabalu (with R)

This morning, we had our first practicals on network flows, using  an example mentioned in some papers published by Noraini Abdullah and Ting Kien Hua, max flow min cut theorem to minimize traffic congestion in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu. From the roads mentioned in the articles, I did try my best to locate the nodes on a map,

m=matrix(c(0,5.995910, 116.105520, 1,5.992737, 116.093718, 2,5.992066, 116.109883, 3,5.976947, 116.095760, 4,5.985766, 116.091580, 5,5.988940, 116.080112, 6,5.968318, 116.080764, 7,5.977454, 116.075460, 8,5.974226, 116.073604, 9,5.969651, 116.073753, 10,5.972341, 116.069270, 11,5.978818, 116.072880),3,12)

we can be visualized below

library(OpenStreetMap) map = openmap(c(lat= 6.000, lon= 116.06), c(lat= 5.960, lon= 116.12)) map=openproj(map) plot(map) points(t(m[3:2,]),col="black", pch=19, cex=3 ) text(t(m[3:2,]),c("s",1:10,"t"),col="white")

If the source is realistic (up north), I do not feel very confortable with the location of the sink (on the west). But let’s pretend it’s find (to do the maths, at least).

To extract information about edge capacity, on that network use the following code that will extract the three tables from the paper

library(devtools) install_github("ropensci/tabulizer") library(tabulizer) location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf' out <- extract_tables(location)

with Windows, it seems to be necessary to download another package first

library(devtools) install_github("ropensci/tabulizerjars") install_github("ropensci/tabulizer") library(tabulizer) location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf' out <- extract_tables(location)

Now we can get out data frame with capacities

B1=as.data.frame(out[[2]]) B2=as.data.frame(out[[3]]) E=data.frame(from=B1[3:20,"V3"], to=B1[3:20,"V4"]) E=E[-c(6,8),] capacity=as.character(B2$V3[-1]) capacity[6]="843" capacity[4]="2913" E$capacity=as.numeric(capacity)

We can add those edges on our map (without the arrows to indicate the direction, it would be to heavy to read)

plot(map) points(t(m[3:2,]),col="black", pch=19, cex=3 ) B=data.frame(i=as.character(c("s",paste("V",1:10,sep=""),"t")), x=m[3,],y=m[2,]) for(i in 1:nrow(E)){ i1=which(B$i==as.character(E$from[i])) i2=which(B$i==as.character(E$to[i])) segments(B[i1,"x"],B[i1,"y"],B[i2,"x"],B[i2,"y"],lwd=3) } text(t(m[3:2,]),c("s",1:10,"t"),col="white")

To get the graph with capacities, an alternative is to use

library(igraph) g=graph_from_data_frame(E) E(g)$label=E$capacity plot(g)

but it does not respect geographical locations of nodes. It can actually be done using

plot(g, layout=as.matrix(B[,c("x","y")]))

To get a better understanding of the capacities of the road, use

plot(g, layout=as.matrix(B[,c("x","y")]), edge.width=E$capacity/200) From that network with capacities, the goal is to determine maximum flow on that network, from the source to the sink. This can be done with R using > (m=max_flow(graph=g, source="s", target="t"))$value [1] 2571

$flow [1] 1191 1380 1422 1380 231 0 231 0 1149 1422 1149 0 0 1149 1422 [16] 1149 Our maximum flow is here 2571, which is different from was is actually claimed both in the two papers max flow min cut theorem to… and application of the Shortest Path… (“the maximum flow for the capacitated network with 12 nodes and 16 edges of the selected scope in this study was 2598 vehicles per hour“) where there are clearly typos since values in the table and on the graph are different. Here I did use the ones from the tables. E$flux1=m$flow E(g)$label=E$flux1 plot(g, layout=as.matrix(B[,c("x","y")]), edge.width=E$flux1/200)

That is nice, but rather odd. Actually, a much simpler flow can be considered, but the same global value

E$flux2=c(1422,1149,1422,1149,0,0,0,0, 1149,1422,1149,0,0,1149,1422,1149) E(g)$label=E$flux2 plot(g, layout=as.matrix(B[,c("x","y")]), edge.width=E$flux2/200)

Nice, isn’t it. It is actually possible to do exactly the same on another paper they have, on the same city, traffic congestion problem of road networks in Kota Kinabalu.

location <- 'http://www.worldresearchlibrary.org/up_proc/pdf/999-150486366625-30.pdf' out <- extract_tables(location) dim(out[[3]]) B1=as.data.frame(out[[3]]) E=data.frame(from=B1[2:61,"V2"], to=B1[2:61,"V3"], capacity=B1[2:61,"V4"]) E$capacity=as.numeric( as.character(E$capacity)) library(igraph) g=graph_from_data_frame(E) m=max_flow(graph=g, source="S", target="T") E$flux1=m$flow E(g)$label=E$flux1 plot(g, edge.width=E$flux1/200, edge.arrow.size=0.15) Here the value of the maximal flow is 4017, just as they found in the original paper # Multinomial Logit as an Iterated Logit Regression For the second section of the course at ENSAE, yesterday, we’ve seen how to run a multinomial logistic regression model. It is simply an extension of the binomial logistic regression. But actually, it is also possible to consider iterative binomial regressions. Consider here a response variable $Y$ with a multinomial distribution (3 factors to have something more general than the binomial), taking values $\{A,B,C\}$, with respective probabilities $\mathbf{p}=(p_A,p_B,p_C)$. Here is a code to generate some multinomial variables msample=function(A,B,C){ Y=rep(NA,B) for(i in 1:B){Y[i]=sample(A,size=1,prob=C[i,])} return(Y) } and here is a code to generate a dataset with $n$ rows, generate3=function(n,x,pb=c(-2,0)){ set.seed(x) X1=runif(n) X2=runif(n) X3=runif(n) s1=pb[1]+X1+X2 s2=pb[2]-X1+X2 P1=exp(s1)/(1+exp(s1)+exp(s2)) P2=exp(s2)/(1+exp(s1)+exp(s2)) Y=msample(0:2,n,cbind(1-P1-P2,P1,P2)) df=data.frame(Y=Y,X1=X1,X2=X2,X3=X3) return(df) } Let us generate a training dataset and a validation one  pb=c(.31,.42) DF1=generate3(1000,1,pb=pb) DF2=generate3(500,2,pb=pb) With a multivariate logistic regression $$\mathbb{P}[Y=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$ $$\mathbb{P}[Y=B|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$ $$\mathbb{P}[Y=B|\mathbf{x}]=\frac{1}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$ For convenience, consider the most popular factor in our training dataset modalite=names(sort(table(DF1$Y),decreasing = TRUE))

Consider a regression model on the simulated dataset (with several covariates), let us estimate it, and let us get predictions.

 library(nnet) reg=multinom(as.factor(Y) ~ ., data = DF1) mp1=predict (reg, DF1, "probs") mp2=predict (reg, DF2, "probs")

An alternative can be the following.
consider a first regression model on the Bernoulli variable $Y_A=\mathbf{1}(Y=A)$. Actually, we will consider the most important factor, but for convenience, assume that it is $A$.
$$\mathbb{P}[Y_A=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}$$
On our dataset, estimate that model, and get predictions. In the case where $Y\neq A$, define another Bernoulli variable $Y_B=\mathbf{1}(Y=B|Y\neq A)$. We can estimate that model and derive two probabilities, $\mathbb{P}(Y=B|Y\neq A)$ and $\mathbb{P}(Y=C|Y\neq A)$ (the sum of the two being equal to 1). Based on those two models, it is possible to compute the three probabilities we are looking for. $\mathbb{P}[Y=A]$ is obtained from the first model, and we can derive the other two from $\mathbb{P}[Y=B|Y\neq A]\cdot\mathbb{P}[Y\neq A]$ and $\mathbb{P}[Y=C|Y\neq A]\cdot\mathbb{P}[Y\neq A]$.

 reg1=glm((Y==modalite[1])~.,data=DF1,family=binomial) reg2=glm((Y==modalite[2])~.,data=DF1[-which(DF1\$Y==modalite[1]),],family=binomial) p11=predict (reg1, newdata=DF1, type="response") p12=predict (reg2, newdata=DF1, type="response") p21=predict (reg1, newdata=DF2, type="response") p22=predict (reg2, newdata=DF2, type="response") mmp1=cbind(p11,(1-p11)*p12,(1-p11)*(1-p12)) mmp2=cbind(p21,(1-p21)*p22,(1-p21)*(1-p22)) colnames(mmp1)=colnames(mmp2)=modalite

Let us compare the predicted probabilites, on the same dataset (here the training dataset)

> mmp1[1:9,c("0","1","2")] 0 1 2 1 0.19728737 0.4991805 0.3035321 2 0.17244580 0.5648537 0.2627005 3 0.19291753 0.5971058 0.2099767 4 0.09087176 0.7787304 0.1303978 5 0.23400225 0.4083022 0.3576955 6 0.18063647 0.6637352 0.1556283 7 0.13188881 0.7402710 0.1278401 8 0.13776970 0.6524959 0.2097344 9 0.12325864 0.6790336 0.1977078 > mp1[1:9,c("0","1","2")] 0 1 2 1 0.19691036 0.5022692 0.3008205 2 0.17123189 0.5680647 0.2607034 3 0.19293066 0.5984402 0.2086291 4 0.08821851 0.7813318 0.1304497 5 0.23470739 0.4109990 0.3542936 6 0.18249687 0.6602168 0.1572863 7 0.13128711 0.7400898 0.1286231 8 0.13525341 0.6553618 0.2093848 9 0.12090016 0.6815915 0.1975084 

The two are very close. So yes, it is possible to see the multinomial regression as some sequential binomial regressions.

# A la recherche des homonymes…

Il y a quelques mois, Baptiste Coulmont m’avait contacté avec une question passionnante (comme chaque fois qu’il me contacte). Baptiste voulait travailler sur la proportion de personnes qui ont un homonymes dans une population de taille donnée, ou sur la probabilité de ne pas avoir d’homonymes dans un bureau de vote, par exemple. Or ce dernier problème n’est pas sans rappeler le “paradoxe des anniversaires“.

Dans un groupe de 23 personnes, il y a plus d’une chance sur deux pour que deux personnes de ce groupe aient leur anniversaire le même jour. Ou dit autrement, il est plus surprenant de ne pas avoir deux personnes qui ont leur anniversaire le même jour que d’avoir deux personnes qui ont leur anniversaire le même jour.

On parle de “paradoxe” parce que le nombre de jours possibles est relativement grand ($k=365$), mais que la taille du groupe est elle relativement petite ($n=23$). Ici c’est un peu la même chose, sauf qu’on est sur des nombres beaucoup plus grands: on a des milliers de noms possibles (en fait, on s’intéresse au couple prénom+nom, ce qui augmente très largement le nombre de possibilités, que l’on notera $k$ ici), on considère des groupes relativement grands (la taille du groupe sera elle notée $n$), et surtout, l’hypothèse (implicite) d’uniformité dans la distribution des dates ne marche pas sur les noms de famille, et les prénoms. On peut malgré tout tenter de formaliser le problème.

Considérons une collection de $k$ éléments, notée $\mathcal{X}=\lbrace x_1,\cdots,x_k\rbrace$, et faisons $n$ tirages avec remise. Les tirages obtenus sont notés ici $X_1,\cdots,X_n$ et correspondent à des variables aléatoires iid, suivant des lois multinomiales sur $\mathcal{X}$ avec comme vecteur de probabilités $\mathbf{p}=(p_1,\cdots,p_k)$. Dans le cas des anniversaires, $\mathcal{X}$ correspond aux dates, avec $k=365$, et $\mathbf{p}$ est uniforme sur $\mathcal{X}$ – ou encore $\mathbf{p}=\mathbf{1}/k$ . Le problème classique est de calculer la probabilité que dans un groupe de taille donnée tout le monde ait des dates d’anniversaire différentes

$$\mathbb{P}[\forall j\neq i, X_i\neq X_j]$$

mais on peut aussi calculer (comme le voulait Baptiste) la proportion de personne ayant un alter-ego (une personne dans le groupe qui a son anniversaire le même jour)
$$\frac{1}{n}\sum_{i=1}^n \mathbf{1}(\exists j\neq i,X_i=X_j)$$

On a donc tenté d’approcher le second, quand on a un groupe de taille $n$. Notons $Z_i$ le nombre de personnes (parmi les $n-1$ restantes) qui partagent la même date d’anniversaire que $i$,
$$Z_i=\sum_{j\in\{1,\cdots,n\}\backslash\{i\}} \mathbf{1}(X_j \neq X_i)$$
Aussi, la personne $i$ a un alter-ego si $Z_i\geq 1$. Et donc, la proportion de personnes avec un alter-ego dans le groupe s’écrit
$$P_{n,k}=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(Z_i\geq 1)$$
Si $\mathbf{p}$ correspond à la loi uniforme sur $\mathcal{X}$, alors
$$1-P_{n,k}=\binom{k}{n}\frac{1}{k^n}\sim \exp\left[\binom{n+1}{2}\frac{1}{k}\right]\sim e^{-n^2/2k}$$
comme le montre Chatterjee et al. (2004). Je pourrais aussi renvoyer aux articles de Inoue & Aki (2008) et plus récemment de  Cortina Borja (2013) pour des revues de littérature sur les aspects computationnels.

Mais la plupart des résultats sont obtenu dans le cas où $\mathbf{p}$ correspond à une distribution uniforme. Les propriétés plus générales sont beaucoup plus rares(on pourra lire Munford (1977), Das Gupta (2004),  Inoue & Aki (2008) ou Nunnikhoven (1992) pour des tentatives. Utiliser des méthodes de Monte Carlo est probablement le plus naturel (si on n’est pas trop pressé) en se donnant un ensemble $\mathcal{X}$ (qui sont les dates dans le problème des anniversaires et les noms dans le cas des homonymes) de taille $k$, et une distribution $\mathbf{p}$ sur $\mathcal{X}$. On a ainsi pu approcher $P_{n,k}$, la proportion de personnes ayant un alter-ego dans une population de taille $n$: sur les figures ci-dessous, on a cette proportion en fonction de  $n$ (en abscisse), avec des courbes différentes pour différents $k$. Histoire d’être un peu réaliste, on va chercher le taux de personnes ayant un homonyme dans une popoluation de plusieurs centaine de milliers de personnes, et quand aussi des milliers, des dizaines de milliers, de possibilités.

Dans le cas d’une distribution uniforme, la proportion évolue comme sur la figure ci-dessous

On a tenté de quitter l’indépendance en utilisant une loi paramétrique, plus précisément une loi de type Pareto, ou Zipf, pour $\mathbf{p}$. En fait, c’est le type de distribution qu’avait observé Li (2012) sur la distribution des prénoms, avec $p_i \propto i^{-\alpha}$ (ou plus récemment, Höhle (2017)). On a alors l’évolution suivante

On est ici en échelle logarithmique pour la taille $n$ parce que c’est naturel. Mais pour l’instant, on était resté sur $P_{n,k}$. Baptiste insistait pour passer aussi à une échelle logarithmique sur l’axe des ordonnées, mais ca n’etait pas concluant. On a alors tenté deux autres représentations, plus classiques quand on travaille sur des probabilités (ou des proportions) : représenter $\text{logit}(P_{n,k})$ en fonction du logarithme de $n$ (où $\text{logit}(u)=\log[u/(1-u)]$) et une version probit, à savoir $\Phi^{-1}(P_{n,k})$ (où ici $\Phi$ est la fonctionc cumulative de la loi normale centrée et réduite).

En grande dimension, on a l’impression d’avoir une relation linéaire entre la transformée probit de la proportion de personnes ayant un homonyme et le logarithme de la taille de la population considérée

$$\Phi^{-1}(P_{n,k}) \sim a_{k}+b_{k} \log[n]$$

C’est joli, mais on se demandait comment utiliser ca sur nos données, pour travailler vraiment sur les homonymes. Le soucis, c’est qu’on n’a pas de fichier national, avec les noms et les prénoms. Alors Baptiste a sorti de son sac les listes électorales parisiennes. Le soucis est que la base est (relativement) petite. Le nombre de couple est très limité. Au lieu de faire des tirages avec remise dans les couples (“prénom”, “nom”) – encore une fois, la base est trop petite – on a voulu tenter de tirer indépendemment le “prénom” et le “nom”. On y gagne en ayant à disposition une base plus grande, mais on introduit un biais colossal, les “prénom” et “nom” étant assez corrélés. Dans notre cas, on a $k_1$ “prénom” et $k_2$ “prénom”. Il y a quelques années Chalabi & Flowers (2014)} (repris plus en détails sur ce blog) était revenu sur la non indépendance des noms et des prénoms, en s’interrogeant sur le nom le plus populaire aux Etats-Unis. Selon les auteurs “Michael Smith” n’était probablement pas le nom le plus populaire, même si “Michael” et “Smith” étaient respectivement les “prénom” et “nom” les plus populaires. En France “Jean” et “Marie” sont les prénoms les plus populaires pour les garcons et les filles, respectivement, et “Martin” est le nom le plus populaire (au moins sur les listes électorales parisiennes), mais en fait, le nom le plus populaire est “Thi Nguyen”, avec 350 personnes à Paris, alors que seulement 84 “Jean Martin” sont observés, le troisième plus haut score. Ah oui, je dois préciser qu’on se limite ici au premier prénom, en cas de prénoms composés. Tout comme Chalabi & Flowers (2014), on a fait un test du chi-deux pour tester l’indépendance, sur les 20 “noms” et “prénoms” les plus populaires. Sur la figure ci-dessous, on a les résidus normalisés du test. Par exemple “Michel Michel”est en rouge, car c’est beaucoup moins donné que ce qu’on aurait si  les deux composantes étaient indépendantes

Bref, on n’a pas trop le droit de tirer le “prénom” et le “nom” indépendament. Mais on va quand même le faire, pour voir. Ici $\mathcal{X}$ est l’ensemble $\mathcal{X}_1\times\mathcal{X}_2$ avec $k_1$ choix pour le “prénom” et $k_2$ pour le “nom”.

Soit $\mathbf{p}=(p_{i_1,i_2})$ la probabilité empirique sur $\mathcal{X}$. On imagine que $\mathbf{p}$ contient énormément de 0, car beaucoup de couples ne sont pas observés. Pour les simulations, on va noter $\mathbf{p}^{\perp}$ la pobabilité jointe sous hypothèse d’indépendance
$$p^{\perp}_{i_1,i_2}=\left(\frac{1}{k_2}\sum_{j=1}^{k_2} p_{i_1,j}\right)\cdot\left(\frac{1}{k_1}\sum_{j=1}^{k_1} p_{j,i_2}\right)=p_{i_1,\cdot}p_{\cdot,i_2}$$

Lors de nos simulations, on va considérer deux cas : on dire des $\mathbf{x}=(x_1,x_2)$ soit suivant la distribution $\mathbf{p}$, soit suivant $\mathbf{p}^{\perp}$. Bon, j’avais formalisé pour clarifier un peu dans ma tête et aussi pour impressionner Baptiste, mais techniquement, dans un cas on tire une paire observée – genre “Baptiste Coulmont” – et dans le second cas, on tire d’un côté un prénom – genre “Baptiste” – et de l’autre un nom de famille – genre “Charpentier” (même s’il n’existe aucun “Baptiste Charpentier”, on dit juste qu’il serait possible qu’il y en ait un).

Bref, on tirant des échantillons (plusieurs milliers de taille donnée, puis en moyennant), on obtient la figure ci-dessous (avec une échelle logarithme pour l’axe des ordonnées)

La courbe noir sur-estime sûrement à cause du faible nombre de choix possibles (on va tirer jusqu’à 50 millions de noms). En revanche, l’indépendance n’étant pas réaliste, la courbe rouge sous-estime probablement la proportion. Mais le plus intéressant est le graphique ci-dessous, obtenu en représentant la version probit de la probabilité

La courbe rouge est incroyablement linéaire ! Bien plus linéaire que ce qu’on espérait avoir suite à nos simulations dans le cas de la loi de Pareto.

En fait, on note que si on se limite au cas où $n$ n’est pas trop grand, la courbe noire semble aussi linéaire. Ensuite on perd la linéarité parce qu’on tire avec remise dans une base trop petite, mais au début (lorsque $n$ n’est pas trop grand), on pourrait accepter l’hypothèse d’une relation linéaire entre la transformée probit, et le logarithme de $n$. Bref, on estime une relation linéaire, et on extrapole pour les valeurs de $n$ plus grandes

Si on revient sur la représentation plus classique, de la proportion de personnes ayant un homonyme en fonction de $n$, on obtient

Autrement dit, en France, on peut estimer que 80% des francais ont un homonyme.

Deux points pour poursuivre, sur la robustesse de l’approche. Tout d’abord, si on change un peu la taille de l’intervalle des $n$ pour lesquels la courbe noire semble linéaire, la pente change peu

et on a un ordre de grandeur assez proche. Bref, ce 80% semble assez robuste,

Ensuite, comme cette relation linéaire me surprenait, on a aussi tenté sur des données américaines… et on a exactement le même genre de résultats…

Pour aller plus loin, on a tout mis sur github (avec des graphs en ggplot2 pour faire plus joli…) et une version (en anglais) est en ligne sur arxiv.

# Traveling Salesman

In the second part of the course on graphs and networks, we will focus on economic applications, and flows. The first series of slides are on the traveling salesman problem. Slides are available online.