{"word_count": 945, "character_count": 9633}
Last week, in our Inequality course, we've been looking at data. We started with some simulated data, only a few of them&#62; library("ineq")&#62; load(url("http://freakonometrics.free.fr/income_5.RData"))&#62; (income=sort(income))[1]  19233  23707  53297  61667 218662How could we say that there is inequality in this sample? If we look at the wealth owned by the poorest, the poorest person (1 out of 5) owns 5% of the wealth; the bottom two (2 out of 5) own 11%, etc&#62; income[1]/sum(income)[1] 0.05107471&#62; sum(income[1:2])/sum(income)[1] 0.1140305&#62; sum(income[1:3])/sum(income)[1] 0.2555648&#62; sum(income[1:4])/sum(income)[1] 0.4193262If we plot those values, we get Lorenz curve&#62; plot(Lc(income))&#62; points(c(0:5)/5,c(0,cumsum(income)/sum(income)),pch=19,col="blue")Now, what if we got 500 observations (and not only 5). In that case, a natural tool to visualize those data (or to be  more specific, their distribution) is the histogram&#62; load(url("http://freakonometrics.free.fr/income_500.RData"))&#62; summary(income)   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    2191   23830   42750   77010   87430 2003000 &#62; hist(log(income),probability=TRUE,col="light blue",border="white")&#62; lines(density(log(income)),col="red")&#62; u=seq(6,15,length=251)lines(u,dnorm(u,mean(log(income)),sd(log(income))),col="blue")Here we use an histogram to visualize our sample. But not on the income, on the logarithm of the income (because of some outliers, we cannot visualize anything on the histogram). Now, it is possible to compute Gini index to get some information about inequalities&#62; gini=function(x){+ n=length(x)+ mu=mean(x)+ g=2/(n*(n-1)*mu)*sum((1:n)*sort(x))-(n+1)/(n-1)+ return(g)}The problem is that, in practice, having an index without any confidence interval can be meaningless. To compute a confidence interval, we can use a bootstrap procedure,&#62; boot=function(x,f,b=500){+ n=length(x)+ F=rep(NA,n)+ for(s in 1:b){+ idx=sample(1:n,size=n,replace=TRUE)+ F[s]=f(x[idx])}+ return(F)}&#62; G=boot(income,gini,1000)&#62; hist(G,col="light blue",border="white",probability=TRUE)The red segments is the 90% confidence interval,&#62; quantile(G,c(.05,.95))       5%       95% 0.4954235 0.5743917 We did include also a blue line with a Gaussian distribution,&#62; segments(quantile(G,.05),1,quantile(G,.95),1,+ col="red",lwd=2)&#62; u=seq(.4,.65,length=251)&#62; lines(u,dnorm(u,mean(G),sd(G)),+ col="blue")Another popular tool is the Pareto plot, where we plot the logarithm of the cumulative survival function against the logarithm of the income,&#62; n=length(income)&#62; x=log(sort(income))&#62; y=log((n:1)/n)&#62; plot(x,y)  If points were on a straight line, it would mean that incomes can be modeled with a Pareto distribution. Obviously, it is not the case here.We've seen previously how to get Lorenz curve. Actually, it is possible to get also Lorenz curve for some parametric distribution, e.g. some log-normal ones,&#62; plot(Lc(income))&#62; lines(Lc.lognorm,param=1.5,col="red")&#62; lines(Lc.lognorm,param=1.2,col="red")&#62; lines(Lc.lognorm,param=.8,col="red")Here, it sounds reasonnable to claim that a log-normal distribution would be a good fit. But maybe not a Pareto distribution&#62; plot(Lc(income))&#62; lines(Lc.pareto,param=2,col="red")&#62; lines(Lc.pareto,param=1.5,col="red")&#62; lines(Lc.pareto,param=1.2,col="red")Actually, it is possible to fit some parametric distributions. Observe that sometimes, we have to change the scale, and use '000s of dollars, instead of dollars,&#62; library(MASS)&#62; fitdistr(income/1e3,"gamma")      shape           rate      1.0812757769   0.0140404379  (0.0604530180) (0.0009868055)Now, consider two distributions, a Gamma one, and a log-normal one (fitted with maximum likelihood techniques)&#62; (fit_g &#60;- fitdistr(income/1e2,"gamma"))      shape           rate      1.0812757769   0.0014040438  (0.0473722529) (0.0000544185)&#62; (fit_ln &#60;- fitdistr(income/1e2,"lognormal"))    meanlog       sdlog     6.11747519   1.01091329  (0.04520942) (0.03196789)We can visualize the densities&#62; u=seq(0,2e5,length=251)&#62; hist(income,breaks=seq(0,2005000,by=5000),+ col=rgb(0,0,1,.5),border="white",+ xlim=c(0,2e5),probability=TRUE)&#62; v_g &#60;- dgamma(u/1e2, fit_g$estimate[1],+ fit_g$estimate[2])/1e2&#62; v_ln &#60;- dlnorm(u/1e2, fit_ln$estimate[1],+ fit_ln$estimate[2])/1e2&#62; lines(u,v_g,col="red",lwd=2)&#62; lines(u,v_ln,col=rgb(1,0,0,.4),lwd=2)Here, it looks like the log-normal is a good candidate. We can also plot the cumulative distribution functions&#62; x &#60;- sort(income)&#62; y &#60;- (1:500)/500&#62; plot(x,y,type="s",col="black",xlim=c(0,250000))&#62; v_g &#60;- pgamma(u/1e2, fit_g$estimate[1],+ fit_g$estimate[2])&#62; v_ln &#60;- plnorm(u/1e2, fit_ln$estimate[1],+ fit_ln$estimate[2])&#62; lines(u,v_g,col="red",lwd=2)&#62; lines(u,v_ln,col=rgb(1,0,0,.4),lwd=2)Now, consider some more realistic situation, where we do not have samples from surveys, but binned data,&#62; load(url("http://freakonometrics.free.fr/income_binned.RData"))&#62; head(income_binned)    low  high number  mean std_err1     0  4999     95  3606     9642  5000  9999    267  7686    14393 10000 14999    373 12505    14714 15000 19999    350 17408    13685 20000 24999    329 22558    14286 25000 29999    337 27584    1520&#62; tail(income_binned)      low   high number   mean std_err46 225000 229999     10 228374    119747 230000 234999     13 232920    137048 235000 239999     11 236341    115749 240000 244999     14 242359    147450 245000 249999     11 247782    148751 250000    Inf    228 395459  189032There is a new package to model that kind of data,&#62; library(binequality)&#62; n &#60;- nrow(income_binned)&#62; fit_LN &#60;- fitFunc(ID=rep("Fake Data",n), + hb=income_binned[,"number"], + bin_min=income_binned[,"low"], + bin_max=income_binned[,"high"],+ obs_mean=income_binned[,"mean"], + ID_name="Country", distribution=LNO, + distName="LNO")Time difference of 2.101471 secsfor LNO fit across 1 distributions Here, we can fit a log-normal distribution (see Methods for estimating inequality from binned incomes for more details about the methodology)&#62; N=income_binned$number&#62; y2=N/sum(N)/diff(income_binned$low) &#62; u=seq(min(income_binned$low),+  max(income_binned$low),length=101)&#62; v=dlnorm(u,fit_LN$parameters[1],+ fit_LN$parameters[2])&#62; plot(u,v,col="blue",type="l",lwd=2)&#62; for(i in 1:(n-1)) rect(income_binned$low[i],0,+ income_binned$high[i],y2[i],col=rgb(1,0,0,.2),+ border="white")Here, on the histogram (since we have binned data, it is natural to plot an histogram), we can see that the fitted log-normal distribution is great. Actually, to be honest, data were simulated from a log-normal distribution, so it makes sense....&#62; N &#60;- income_binned$number&#62; y1 &#60;- cumsum(N)/sum(N) &#62; u &#60;- seq(min(income_binned$low),+ max(income_binned$low),length=101)&#62; v &#60;- plnorm(u,fit_LN$parameters[1],+ fit_LN$parameters[2])&#62; plot(u,v,col="blue",type="l",lwd=2)&#62; for(i in 1:(n-1)) rect(income_binned$low[i],0,+ income_binned$high[i], y1[i],col=rgb(1,0,0,.2))&#62; for(i in 1:(n-1)) rect(income_binned$low[i],+ y1[i],income_binned$high[i],c(0,y1)[i],+ col=rgb(1,0,0,.4))For the cumulated distribution function, I consider either the worst case scenario (everyone in a bin get the lower possible income) and the best case (everyone has the highest possible income).It is also possible to fit all standard GB2 distributions&#62; fits=run_GB_family(ID=rep("Fake Data",n),+ hb=income_binned[,"number"],+ bin_min=income_binned[,"low"],+ bin_max=income_binned[,"high"],+ obs_mean=income_binned[,"mean"],+ ID_name="Country")To get the best candidate, look at&#62; fits$fit.filter[,c("gini","aic","bic")]Now, that was fine, but those were simulated data. What about real data now? Here we go&#62; data = read.table(+ "http://freakonometrics.free.fr/us_income.txt",+ sep=",",header=TRUE)&#62; head(data)    low  high number_1000s  mean std_err1     0  4999         4245  1249      502  5000  9999         5128  7923      303 10000 14999         7149 12389      284 15000 19999         7370 17278      265 20000 24999         7037 22162      276 25000 29999         6825 27185      28&#62; n &#60;- nrow(data)&#62; fit_LN &#60;- fitFunc(ID=rep("US",n), + hb=data[,"number_1000s"], + bin_min=data[,"low"],+ bin_max=data[,"high"], + obs_mean=data[,"mean"], + ID_name="Country", + distribution=LNO, distName="LNO")Time difference of 0.1855791 secsfor LNO fit across 1 distributionsAgain, I try to fit a log-normal distribution&#62; N=data$number&#62; y2=N/sum(N)/diff(data$low) &#62; u=seq(min(income_binned$low),+ max(income_binned$low),length=101)&#62; v=dlnorm(u,fit_LN$parameters[1],+ fit_LN$parameters[2])&#62; plot(u,v,col="blue",type="l",lwd=2)&#62; for(i in 1:(n-1)) rect(data$low[i],+ 0,data$high[i],y2[i],col=rgb(1,0,0,.2))But here, the fit is rather poor. Again, we can estimate all standard GB2 distribution&#62; &#62; fits=run_GB_family(ID=rep("US",n), + hb=data[,"number_1000s"], + bin_min=data[,"low"],+ bin_max=data[,"high"],+ obs_mean=data[,"mean"],+ ID_name="Country")We can get Gini index, AIC and BIC&#62; fits$fit.filter[,c("gini","aic","bic")]       gini      aic      bic1  4.413431 825368.5 825407.32  4.395080 825598.8 825627.93  4.451881 825495.7 825524.84  4.480850 825881.7 825910.85  4.417276 825323.6 825352.76  4.922122 832408.2 832427.67  4.341036 827065.2 827084.68  4.318667 826112.8 826132.29        NA 831054.2 831073.610       NA       NA       NAto see that the best distribution seems to be a Generalized Gamma.