Gini

> options(echo=FALSE, encoding="UTF-8")
Loading required package: readstata13
> options(error = expression(q('no')))
> gini <- function(x,weight) { 
+ ox <- order(x) 
+ x <- x[ox] 
+ weight <- weight[ox]/sum(weight) 
+ p <- cumsum(weight) 
+ nu <- cumsum(weight*x) 
+ n <- length(nu) 
+ nu <- nu/nu[n] 
+ res <- sum(nu[-1]*p[-n])-sum(nu[-n]*p[-1]) 
+ return(res) 
+ } 
> wNtile<- function(var,wgt,split){ 
+ x <- var[order(var)] 
+ y <-wgt[order(var)] 
+ z <- cumsum(y) / sum(y)  
+ cop <- rep(NA,length(split))  
+ for (i in 1:length(cop)) { 
+ cop[i] <- x[Find(function(h) z[h] > split[i], seq_along(z))] } 
+ return(cop) } 
> topBottom<- function(var, botline, topline) { 
+ tb <- ifelse(var < botline, botline, var)  
+ tb[tb > topline] <- topline 
+ return(tb) 
+ } 
> setups <- function(data_file) { 
+ vars <- c('dhi', 'hifactor', 'hpub_i', 'hpub_u', 'hpub_a', 'hiprivate', 'hxitsc', 'hpopwgt', 'nhhmem', 'grossnet') 
+ subset <-'complete.cases(dhi,hifactor,hpub_i,hpub_u,hpub_a,hiprivate,hxitsc)'  
+ df <- read.LIS(data_file, labels=FALSE, vars=vars, subset=subset) 
+ botline <- 0 
+ topline <- 10 * wNtile(df$dhi, df$hpopwgt, 0.5) 
+ df$dhi <- topBottom(df$dhi, botline, topline)  
+ df$edhi <- df$dhi / df$nhhmem^0.5 
+ df$cdhi <- df$dhi / df$nhhmem 
+ return(df) 
+ } 
> df <- setups('ru16h') 
[1] "Loading dataset ru16h..."
> round(gini(df$dhi , df$hpopwgt) , digits = 4)  
[1] 0.366
> round(gini(df$cdhi, df$hpopwgt*df$nhhmem), digits = 4)  
[1] 0.3571
> round(gini(df$edhi, df$hpopwgt*df$nhhmem), digits = 4) 
[1] 0.3305
> 
> proc.time()
   user  system elapsed 
  2.569   0.212   3.288