A few days ago, Jean-François Mignot published an interesting article entitled Tour de France 2014 : pourquoi le vainqueur gagne 100 fois plus que le 10e. In this article, we have the following graph, with the income of the cyclist, as a function of his final ranking (the data where downloaded from http://sportbuzzbusiness.fr/)

> bike=read.csv(
+ "http://freakonometrics.free.fr/tourdefrance.csv",
+ sep=";",header=TRUE,dec=" ")

> bike[1:19,"prime"]=bike[1:19,"prime"]*1000
> plot(bike,log="y",type="b",
+ xlab="(Final) rank",ylab="Bonus")

As pointed out by Jean-François, if the winner gets a lot of money, the bonus decreases fast, very fast actually. Gini index is very high here

> library(ineq)
> ineq(X,type="Gini")
[1] 0.910461

and if we look at Lorenz curve, indeed, the Tour de France is not very equalitarian,

> plot(Lc(X),col="red",lwd=2)

If we look at the Pareto plot (with the rank versus the bonus, both in log scales),

> X=bike$prime
> n=length(X)
> plot(X,((1:n)/(n+1)),type="b",log="xy",
+ ylab="(Final) rank (log scale)",
+ xlab="Bonus (log scale)")

We have a distribution with heavy tails, and thus, not equalitarian.

If we try to fit a lognormal distribution

> library(fitdistrplus)
> fitdistr(X,"lnorm")
> plot(X,((1:n)/(n+1)),type="b",log="xy",
+ ylab="(Final) rank (log scale)",
+ xlab="Bonus (log scale)")
> lines(qlnorm(1-(1:m)/(m+1),
+ parlnorm[1],parlnorm[2]),(1:m)/(m+1),
+ col="blue")

the fit is not good. Here we clearly underestimate the bonus for the winer (and for the top ten cyclist). We got the same with a paralogistic distribution,

> library(VGAM)
> parparalog=Coef(vglm(X ~ 1,paralogistic))
> plot(X,((1:n)/(n+1)),type="b",log="xy",
+ ylab="(Final) rank (log scale)",
+ xlab="Bonus (log scale)")
> lines(qparalogistic(1-(1:m)/(m+1),
+ parparalog[1],parparalog[2]),(1:m)/(m+1),
+ col="blue")

We can play hours trying to find the best distribution for the bonuses. Let’s try a last one : the (generalized) Pareto distribution. This distribution was introduced to model incomes, and it has very nice properties to model extremal events (very large bonuses). Its power tails are related to power decay of the size of the bonus we do observe here.

> library(evir)
> pargpd=evir::gpd(X,min(X))$par.ests
> plot(X,((1:n)/(n+1)),type="b",log="xy",
+ ylab="(Final) rank (log scale)",
+ xlab="Bonus (log scale)")
> m=1000
> lines(evir::qgpd(1-(1:m)/(m+1),xi=pargpd[1],
+ mu=min(X),beta=pargpd[2]),(1:m)/(m+1),col="blue")

The fit is not perfect. And if we focus only on bonuses above 700 (and not 450), the fit is slightly better

> library(evir)
> pargpd2=evir::gpd(X,700)$par.ests
> plot(X,((1:n)/(n+1)),type="b",log="xy",
+ ylab="(Final) rank (log scale)",xlab="Bonus (log scale)")
> m=1000
> lines(evir::qgpd(1-(1:m)/(m+1),xi=pargpd2[1],
+ mu=700,beta=pargpd2[2]),(1:m)/(m+1)*mean(X>700),col="blue")

Clearly, fitting a parametric distribution to model the income is not that simple. I’ll get back to this point soon. With Emmanuel Flachaire, we are finalizing a paper on nonparametric estimation of income distribution, to improve the estimation of inequality indices. Because the value of Gini index, for instance, is not very robust here, with 150 observations from a Pareto distribution. If we generate 1,000 samples from the Pareto distribution that was fitted, our index is very volatile,

> M=matrix(evir::rgpd(150*1000,xi=pargpd[1],
+ mu=min(X),beta=pargpd[2]),150,1000)
> gini=function(x) ineq(x,type="Gini")
> boxplot(apply(M,2,gini),horizontal=TRUE)
> abline(v=ineq(X,type="Gini"),col="red")

But sometimes, we have no choice, and using a parametric distribution is the only option. I will upload a post on that topic, based on a study we did with Dr Goulu.

---

---