Tag Archives: computer

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

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?

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

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.