原文鏈接：http://tecdat.cn/?p=20613

原文出處：拓端數(shù)據(jù)部落公眾號

洛倫茲曲線來源于經(jīng)濟(jì)學(xué)，用于描述社會收入不均衡的現(xiàn)象。將收入降序排列，分別計算收入和人口的累積比例。
本文，我們研究收入和不平等。我們從一些模擬數(shù)據(jù)開始

1. 


2. > (income=sort(income))

3. [1] ?19246 ?23764 ?53237 ?61696 218835

為什么說這個樣本中存在不平等？如果我們看一下最貧窮者擁有的財富，最貧窮的人（五分之一）擁有5％的財富；倒數(shù)五分之二擁有11％，依此類推

1. > income[1]/sum(income)

2. [1] 0.0510

3. > sum(income[1:2])/sum(income)

4. [1] 0.1140

如果我們繪制這些值，就會得到?洛倫茲曲線

1. 


2. > plot(Lorenz(income))

3. > points(c(0:5)/5,c(0,cumsum(income)/sum(income))

現(xiàn)在，如果我們得到500個觀測值。直方圖是可視化這些數(shù)據(jù)分布的方法

1. > summary(income)

2. Min. 1st Qu. ?Median ? ?Mean 3rd Qu. ? ?Max.

3. 2191 ? 23830 ? 42750 ? 77010 ? 87430 2003000

4. > hist(log(income),

在這里，我們使用直方圖將樣本可視化。但不是收入，而是收入的對數(shù)（由于某些離群值，我們無法在直方圖上可視化）?，F(xiàn)在，可以計算?基尼系數(shù)?以獲得有關(guān)不平等的一些信息

1. > gini=function(x){

2. 


3. + mu=mean(x)

4. + g=2/(n*(n-1)*mu)*sum((1:n)*sort(x))-(n+1)/(n-1)

實際上，沒有任何置信區(qū)間的系數(shù)可能毫無意義。計算置信區(qū)間，我們使用boot方法

1. 


2. > G=boot(income,gini,1000)

3. > hist(G,col="light blue",border="white"

紅色部分是90％置信區(qū)間，

1. 


2. 5% ? ? ? 95%

3. 0.4954235 0.5743917

還包括了一條具有高斯分布的藍(lán)線，

1. > segments(quantile(G,.05),1,quantile(G,.95),1,

2. 


3. > lines(u,dnorm(u,mean(G),sd(G)),

另一個流行的方法是帕累托圖（Pareto plot），我們在其中繪制了累積生存函數(shù)的對數(shù)與收入的對數(shù)，

1. 


2. > plot(x,y)

如果點在一條直線上，則意味著可以使用帕累托分布來建模收入。

前面我們已經(jīng)看到了如何獲得洛倫茲曲線。實際上，也可以針對某些參數(shù)分布（例如，一些對數(shù)正態(tài)分布）獲得Lorenz曲線，

1. 


2. > lines(Lc.lognorm,param=1.5,col="red")

3. > lines(Lc.lognorm,param=1.2,col="red")

4. > lines(Lc.lognorm,param=.8,col="red")

在這里， 對數(shù)正態(tài)分布是一個很好的選擇。帕累托分布也許不是：

1. 


2. > lines(Lc.pareto,param=1.2,col="red")

實際上，可以擬合一些參數(shù)分布。

1. 


2. shape ? ? ? ? ? rate

3. 1.0812757769 ? 0.0140404379

4. (0.0604530180) (0.0009868055)

現(xiàn)在，考慮兩種分布，伽馬分布和對數(shù)正態(tài)分布（適用于極大似然方法）

1. 


2. shape ? ? ? ? ? rate

3. 1.0812757769 ? 0.0014040438

4. (0.0473722529) (0.0000544185)

5. meanlog ? ? ? sdlog

6. 6.11747519 ? 1.01091329

7. (0.04520942) (0.03196789)

我們可以可視化密度

1. > hist(income,breaks=seq(0,2005000,by=5000),

2. + col=rgb(0,0,1,.5),border="white",

3. + fit_g$estimate[2])/1e2

4. + fit_ln$estimate[2])/1e2

5. > lines(u,v_g,col="red",lwd=2)

6. > lines(u,v_ln,col=rgb(1,0,0,.4),lwd=2)

在這里，對數(shù)正態(tài)似乎是一個不錯的選擇。我們還可以繪制累積分布函數(shù)

1. > plot(x,y,type="s",col="black",xlim=c(0,250000))

2. 


3. + fit_g$estimate[2])

4. 


5. + fit_ln$estimate[2])

6. > lines(u,v_g,col="red",lwd=2)

現(xiàn)在，考慮一些更現(xiàn)實的情況，在這種情況下，我們沒有來自調(diào)查的樣本，但對數(shù)據(jù)進(jìn)行了合并，

對數(shù)據(jù)進(jìn)行建模，

1. fit(ID=rep("Data",n),

2. 


3. 


4. 


5. Time difference of 2.101471 secs

6. for LNO fit across 1 distributions

我們可以擬合對數(shù)正態(tài)分布（有關(guān)該方法的更多詳細(xì)信息，請參見?從合并收入估算不平等?的方法）

1. > y2=N/sum(N)/diff(income_binned$low)

2. 


3. + fit_LN$parameters[2])

4. > plot(u,v,col="blue",type="l",lwd=2)

5. > for(i in 1:(n-1)) rect(income_binned$low[i],0,

6. + income_binned$high[i],y2[i],col=rgb(1,0,0,.2),

在此，在直方圖上（由于已對數(shù)據(jù)進(jìn)行分箱，因此很自然地繪制直方圖），我們可以看到擬合的對數(shù)正態(tài)分布很好。

1. > v <- plnorm(u,fit_LN$parameters[1],

2. + fit_LN$parameters[2])

4. > for(i in 1:(n-1)) rect(income_binned$low[i],0,

5. 


6. 


7. > for(i in 1:(n-1)) rect(income_binned$low[i],

8. + y1[i],income_binned$high[i],c(0,y1)[i],

9. 


對于累積分布函數(shù)，我考慮了最壞的情況（每個人都處于較低的收入中）和最好的情況（每個人都具有最高可能的收入）。

也可以擬合廣義beta分布

GB_family(ID=rep("Fake Data",n),

為了獲得最佳模型，查看

> fits[,c("gini","aic","bic")]

結(jié)果很好，接下來看下真實數(shù)據(jù)：

1. fit(ID=rep("US",n),

2. 


3. + distribution=LNO, distName="LNO"

4. Time difference of 0.1855791 secs

5. for LNO fit across 1 distributions

同樣，我嘗試擬合對數(shù)正態(tài)分布

1. > v=dlnorm(u,fit_LN$parameters[1],

2. 


3. > plot(u,v,col="blue",type="l",lwd=2)

4. > for(i in 1:(n-1)) rect(data$low[i],

5. 


但是在這里，擬合度很差。同樣，我們可以估算廣義beta分布

1. >

2. GB_family(ID=rep("US",n),

3. 


4. + ID_name="Country")

可以得到基尼指數(shù)，??AIC?和BIC

1. gini ? ? ?aic ? ? ?bic

2. 1 ?4.413431 825368.5 825407.3

3. 2 ?4.395080 825598.8 825627.9

4. 3 ?4.451881 825495.7 825524.8

5. 4 ?4.480850 825881.7 825910.8

6. 5 ?4.417276 825323.6 825352.7

7. 6 ?4.922122 832408.2 832427.6

8. 7 ?4.341036 827065.2 827084.6

9. 8 ?4.318667 826112.8 826132.2

10. 9 ? ? ? ?NA 831054.2 831073.6

11. 10 ? ? ? NA ? ? ? NA ? ? ? NA

看到最好的分布似乎是?廣義伽瑪分布。

最受歡迎的見解

1.R語言泊松Poisson回歸模型分析案例

2.R語言進(jìn)行數(shù)值模擬：模擬泊松回歸模型

3.r語言泊松回歸分析

4.R語言對布豐投針（蒲豐投針）實驗進(jìn)行模擬和動態(tài)可視化

5.用R語言模擬混合制排隊隨機(jī)服務(wù)排隊系統(tǒng)

6.GARCH（1,1），MA以及歷史模擬法的VaR比較

7.R語言做復(fù)雜金融產(chǎn)品的幾何布朗運動的模擬

8.R語言進(jìn)行數(shù)值模擬：模擬泊松回歸模型

9.R語言對巨災(zāi)風(fēng)險下的再保險合同定價研究案例：廣義線性模型和帕累托分布Pareto distributions