####What are the drivers of the income inequality for women in urban China in 1995 and 2002?#####
#MASTER THESIS
#KULeuven
#Marine Panozzo
#2019-2020

#INSTALL PACKAGE#####

install.packages("haven")
install.packages("gglorenz") #lorenz curve
install.packages("ineq") #lorenzcurve
install.packages("DescTools")
install.packages("IC2") #concentracurve
install.packages("tidyverse") #cross val
install.packages("caret") #cross val
install.packages("randomForest") #cross val
install.packages("dineq") #regbydecomp
install.packages("cran")
install.packages("latticeExtra")
install.packages("Hmisc") #datadescription
install.packages("car")
install.packages("MVN")
install.packages("xtable") #export
install.packages("relaimpo") #shapley
devtools::install_github("elbersb/shapley") #shapley
install.packages("remotes")
remotes::install_github("brechtdv/rineq") 
install.packages("corrplot") #correlation
install.packages("broom")
installed.packages("boot") #bootstrap

#LIBRARY#####
library(haven)
library(gglorenz)
library(ineq)
library(DescTools)
library(IC2)
library(tidyverse)
library(caret)
library(randomForest)
library(dineq)
library(cran)
library(latticeExtra)
library(Hmisc)
library(car)
library(MVN)
library(xtable)
library(relaimpo)
library(shapley)
library(rineq)
library(corrplot)
library(broom)
library(boot)



#URBAN 1995####
##IMPORT####
library(haven)
X03012_0003_Data <-read_dta("C:/Users/stagiaire/Desktop/data/03012-0003-Data.dta")
DATA_1995_0003 <- as.data.frame(X03012_0003_Data)


#MALE AND FEMALE DATAFRAME####

#put gender variable as factor
DATA_1995_0003$A4 <- as.factor(DATA_1995_0003$A4) 

#create 2 dataset w.r.t gender factor
gender <- split(DATA_1995_0003, DATA_1995_0003$A4) 

Data1995_male <- as.data.frame(gender$`1`) 
Data1995_female <- as.data.frame(gender$`2`)

#INCOME1995#####

Income_1995_female <- Data1995_female$A51
Income_1995_female

#summary
summary(Income_1995)
#LORENZ CURVE#####

library(ineq)
#lorenzcurve women 

Income_1995_naomitfemale <- na.omit(Income_1995_female)
plot(Lc(Income_1995_naomitfemale), main = "Lorenz Curve Women Urban in 1995", xlab="Percent Population (cumulative)",
     ylab="Percent of area (%)")


#DAVIDSON GINI#####

y<-na.omit(Income_1995_female)

# Asymptotic SD Davidson

stdavidson = function(y){
  n=length(y)
  G=Gini(y)
  mu=mean(y)
  zhat=rep(0,n)
  for(i in 1:n){
    zhat[i] = -(G+1)*y[i]+ ((2*i-1)/n)*y[i]-2*sum(y[1:i])/n
  }
  zbar = mean(zhat)
  varG = sum((zhat-zbar)^2)/(n*mu)^2
  return (sqrt(varG))
}

install.packages("DescTools")
library(DescTools)
gini1995 <- Gini(y, unbiased = TRUE) #0.349671
segini1995 <- stdavidson(y) #0.005818102

#confidence interval
upperg1995 <- gini1995 + (1.96*segini1995)
lowerg1995 <- gini1995 - (1.96*segini1995)
cbind(lowerg1995,upperg1995)


#CONCENTRATION CURVE#####

#define new variables from dataset

Wage_ubran_1995 <- Data1995_female$A52
businessdataframe <- data.frame(Data1995_female$A64, Data1995_female$A65, Data1995_female$A66)
Business_urban_1995 <- rowSums(businessdataframe, na.rm = TRUE)
Labour_urban_1995 <- Data1995_female$A69
Property_urban_1995 <- Data1995_female$A71
Transfer_urban_1995 <- Data1995_female$A76
otherdataframe <- data.frame(Data1995_female$A62, Data1995_female$A67, Data1995_female$A68, Data1995_female$A85)
Other_urban_1995 <- rowSums(otherdataframe, na.rm = TRUE)  

#create dataframe for concentration curve
ccurve1995 <- cbind(Income_1995_female, Wage_ubran_1995, Business_urban_1995, Labour_urban_1995, Property_urban_1995,Transfer_urban_1995, Other_urban_1995)

install.packages("IC2")
library(IC2)
#source https://rdrr.io/cran/IC2/man/curveConcent.html

#four curves in the same plot : limited plot

curveConcent(x=ccurve1995[,2], y=ccurve1995[,1], col="red",lty=1,xlab="Cumulative % of population, ordered by earnings",ylab="Cumulative % of income variable")
curveConcent(x=ccurve1995[,3], y=ccurve1995[,1], col = "green",lty=1,xlab="Cumulative % of population, ordered by earnings", add= TRUE)
curveConcent(x=ccurve1995[,5], y=ccurve1995[,1], col = "blue",lty=1,xlab="Cumulative % of population, ordered by earnings", add=TRUE)
curveConcent(x=ccurve1995[,6], y=ccurve1995[,1], col = "orange",lty=1,xlab="Cumulative % of population, ordered by earnings", add=TRUE)
title(main="Concentration Curves 1995")
legend("topleft", legend = c("Wage", "Business", "Property", "Transfer"), col=c("red","green", "blue", "orange"), lty=rep(1, 4),ncol=1)


#six curves in the same plot ; full plot

curveConcent(x=ccurve1995[,2], y=ccurve1995[,1], col="red",lty=1,xlab="Cumulative % of population, ordered by earnings",ylab="Cumulative % of income variable")
curveConcent(x=ccurve1995[,3], y=ccurve1995[,1], col = "green",lty=1,xlab="Cumulative % of population, ordered by earnings", add= TRUE)
curveConcent(x=ccurve1995[,4], y=ccurve1995[,1], col = "brown",lty=1,xlab="Cumulative % of population, ordered by earnings", add= TRUE)
curveConcent(x=ccurve1995[,5], y=ccurve1995[,1], col = "blue",lty=1,xlab="Cumulative % of population, ordered by earnings", add=TRUE)
curveConcent(x=ccurve1995[,6], y=ccurve1995[,1], col = "orange",lty=1,xlab="Cumulative % of population, ordered by earnings", add=TRUE)
curveConcent(x=ccurve1995[,7], y=ccurve1995[,1], col = "black",lty=1,xlab="Cumulative % of population, ordered by earnings", add=TRUE)
title(main="Concentration Curves 1995")
legend("topleft", legend = c("Wage", "Business", "Labour", "Property","Transfer", "Other"), col=c("red","green","brown", "blue", "orange","black"), lty=rep(1, 4),ncol=1)




#CONCENTRATION INDEX#####

calcSConc(Wage_ubran_1995,Income_1995_female) #0.294513
calcSConc(Business_urban_1995,Income_1995_female) #0.2739602 
calcSConc(Labour_urban_1995,Income_1995_female) #0.1036749 
calcSConc(Transfer_urban_1995,Income_1995_female) #0.1294597 
calcSConc(Property_urban_1995,Income_1995_female) #0.3535982
calcSConc(Other_urban_1995,Income_1995_female) #0.5067159 

cc1995 <- c(0.2739602,0.294513,0.1036749,0.1294597,0.3535982,0.5067159)

#CORRELATION INCOME SOURCE######

ccurve1995 <- as.data.frame(ccurve1995)
#name the columns
colnames(ccurve1995) <- c("Income", "Wage", "Business", "Labour", "Property", "Transfer", "Others")
#correlation
library(corrplot)
mcor <- cor(na.omit(ccurve1995))
#change the range of color
col3 <- colorRampPalette(c("white", "#6666FF", "#000033")) 
#figure correlation
corrplot(mcor, method = "number", type="upper", tl.col="black", tl.srt = 45, col = col3(20) )
#pearson test
library(Hmisc)
corrpvalue95 <- rcorr(as.matrix(ccurve1995), type = "pearson")
corrpvalue95 <- corrpvalue95$P
#export results
library(xtable)
print(xtable(as.data.frame(corrpvalue95),type = "latex"), file = "pvaluecorr1995v2.tex")

#VARIABLES#####

#defining variables
#extract variables from the whole dataset

Province_1995 <- factor(Data1995_female$'PROVINCE')
Education_1995 <- factor(Data1995_female$A11)
CPC_1995 <- factor(Data1995_female$A9)

#putting NA value for the level "0" allow to delete the value 

Ownership_1995 <- Data1995_female$A27
Ownership_1995[Ownership_1995 == 0] <- NA #0 value are undefined value
Ownership_1995 <- factor(Ownership_1995)

Occup_1995 <- Data1995_female$A31
Occup_1995[Occup_1995 == 0] <- NA  #0 value are undefined value
Occup_1995 <- factor(Occup_1995)

Industry_1995 <- Data1995_female$A34
Industry_1995[Industry_1995 == 0] <- NA  #0 value are undefined value
Industry_1995 <- factor(Industry_1995)

#putting zero values equal to NA because regression with ln(0) is impossible
Income_1995_female[Income_1995_female == 0] <- NA  #0 value are undefined value


#make work experience categorical 

Workexp_1995 <- as.numeric(Data1995_female$A22)
Workexpfactor1995 <- cut(
  Workexp_1995,
  breaks = c(0, 11, 21, 31, 41, 100),
  labels = c("<=10", "[11:20]", "[21:30]", "[31:40]", ">=41"),
  right  = FALSE
)


#VARIABLE DESCRIPTION ######

#create data set with all the variables

DATAREG1995 <- cbind(Income_1995_female,Province_1995,Education_1995,CPC_1995,Workexpfactor1995,Ownership_1995,Occup_1995,Industry_1995)
#take only complete info
DATAREG1995 <- na.omit(DATAREG1995)
DATAREG1995 <- as.data.frame(DATAREG1995)

#mean and se for each category
mean1995 <- sapply(DATAREG1995, mean, na.rm = T)
se1995 <- sapply(DATAREG1995, sd,na.rm= T)

#frequency and proportion of the category. It is to calculate the mean income for each category below

describe(DATAREG1995[c("Province_1995","CPC_1995","Education_1995","Workexpfactor1995","Ownership_1995","Occup_1995","Industry_1995")])

#Mean income for each category within each variables

aggprov95<-aggregate(DATAREG1995$Income_1995_female, by=list(CPC_1995=DATAREG1995$Province_1995), FUN=sum)
freqprov95<-c(477,   571,   639,   723,   486,   563,   639,   548,   755,   614,   328)
aggprov95/freqprov95

aggcpc1995<-aggregate(DATAREG1995$Income_1995_female, by=list(CPC_1995=DATAREG1995$CPC_1995), FUN=sum)
freqcpc95<-c( 1015, 5328)
aggcpc1995/freqcpc95
  
aggedu95<-aggregate(DATAREG1995$Income_1995_female, by=list(CPC_1995=DATAREG1995$Education_1995), FUN=sum)
freqedu95<-  c( 303,   681,  1043,  1511,  2068,   583,   154)
aggedu95/freqedu95

aggwork95<-aggregate(DATAREG1995$Income_1995_female, by=list(CPC_1995=DATAREG1995$Workexpfactor1995), FUN=sum)
freqwork95<-  c(   1309,  2080,  2195,   694,    65)
aggwork95/freqwork95

aggown95<-aggregate(DATAREG1995$Income_1995_female, by=list(CPC_1995=DATAREG1995$Ownership_1995), FUN=sum)
freqown95<- c(1424,  3360,  1333,    20,    90,    62,     7,     1 ,   46)
aggown95/freqown95
  
aggocc95<-aggregate(DATAREG1995$Income_1995_female, by=list(CPC_1995=DATAREG1995$Occup_1995), FUN=sum) 
freqocc95<-c(  88,     9,  1366,   101,   255,  1362,  1174,  1508,   480)
aggocc95/freqocc95


aggind95<-aggregate(DATAREG1995$Income_1995_female, by=list(CPC_1995=DATAREG1995$Industry_1995), FUN=sum) 
freqindu95<-c( 86,  2665,    58,   160,   248,  1050,   286,   371,   528,   124,   121,   590,    56)
aggind95/freqindu95

#LN REGRESSION#########

#creating new variables without considered as.factor before putting them in the regression
#R prefers to declare the variable as dummy directly in the regression 

Province_1995b <- Data1995_female$'PROVINCE'
Education_1995b <- Data1995_female$A11
CPC_1995b <- Data1995_female$A9
Workexp_1995b <- Data1995_female$A22

Ownership_1995b <- Data1995_female$A27
Ownership_1995b[Ownership_1995b == 0] <- NA #0=unknown code

Occup_1995b <- Data1995_female$A31
Occup_1995b[Occup_1995b == 0] <- NA #0=unknown code

Industry_1995b <- Data1995_female$A34
Industry_1995b[Industry_1995b == 0] <- NA #0=unknown code

Income_1995_female[Income_1995_female == 0] <- NA

#create new dataset

DATA1995 <- cbind(Income_1995_female,Province_1995b, Education_1995b, CPC_1995b, Workexpfactor1995, Ownership_1995b, Occup_1995b, Industry_1995b)
DATA1995 <- na.omit(DATA1995)
DATA1995 <- as.data.frame(DATA1995)

#regression defining variables as factor directly into the command

reg1995<-summary(lm(log(Income_1995_female)~ factor(Province_1995b)+ factor(CPC_1995b) + factor(Education_1995b) + factor(Workexpfactor1995) + factor(Ownership_1995b) + factor(Occup_1995b) +factor(Industry_1995b), data = DATA1995))

#exportation
library(broom)
write.csv( tidy( reg1995 ) , "coefs1995.csv" )
write.csv( glance( reg1995 ) , "an1995.csv" )


#CROSSVALIDATION####

library(tidyverse)
library(caret)

set.seed(123)
train.control <- trainControl(method = "cv", number=10)

model1995full <- train(log(Income_1995_female)~ factor(Province_1995b)+ factor(CPC_1995b) + factor(Education_1995b) + factor(Workexpfactor1995) + factor(Ownership_1995b) + factor(Occup_1995b) +factor(Industry_1995b), data = DATA1995,
                       trControl = train.control, method="lm")
model1995full

###REG BY DECOMP #######

library(dineq)

#fields method

ineq_decomp1995 <- dineq_rb(Income_1995_female~factor(Province_1995b)+ 
                              factor(CPC_1995b) + factor(Education_1995b) +
                              factor(Workexpfactor1995) + factor(Ownership_1995b) +
                              factor(Occup_1995b) +factor(Industry_1995b)
                            ,data = DATA1995)
#results
ineq_decomp1995["decomposition_inequality"]


#SHAPLEY VALUE##############

library(shapley)

#function to incorporate in the shapley function

reg <- function(factors, dv, data) {
  if (length(factors) == 0) return(0)
  formula <- paste(dv, "~", paste(factor(factors), collapse = "+"))
  m <- lm(formula, data = data)
  summary(m)$r.squared
}
#recreating dataframe with ln(income)
Income1995log <- log(Income_1995_female)
DATAREG1995shapl <- cbind(Income1995log,Province_1995b,Education_1995b,CPC_1995b,Workexpfactor1995,Ownership_1995b,Occup_1995b,Industry_1995b)
DATAREG1995shapl <- na.omit(DATAREG1995shapl)
DATAREG1995shapl <- as.data.frame(DATAREG1995shapl)

#Shapley value decomposition
shapley1995 <- shapley(reg, c("factor(Province_1995b)", "factor(Education_1995b)", "factor(CPC_1995b)", "Workexpfactor1995", "factor(Ownership_1995b)", "factor(Occup_1995b)", "factor(Industry_1995b)"), silent = TRUE, dv="Income1995log", data= DATAREG1995shapl)
shapley1995

weight1995 <- as.vector(shapley1995$value)
res1995 <- 1- sum(weight1995)
res1995

netcontrib1995s <- weight1995/res1995
netcontrib1995s

#export
print(xtable(as.data.frame(shapley1995),type = "latex"), file = "shapley1995.tex")

#BOOTSTRAP######


library(shapley)

#create function to put in argument in the 'boot' function
st = function(y, indices){
  u <- DATAREG1995shapl[indices,]
  v = shapley(reg, c("factor(Province_1995b)", "factor(Education_1995b)", "factor(CPC_1995b)", "Workexpfactor1995", "factor(Ownership_1995b)", "factor(Occup_1995b)", "factor(Industry_1995b)"), silent = TRUE, dv="Income1995log",data = u)
  return(v$value)
}

library(boot)
#100 re-sampling bootstrap
results1995 <- boot(data=DATAREG1995shapl,st, R=100)

#stock estimators
estimators1995 <- c(0.14059057,0.04479386,0.01349264,0.04733561,0.03539168,0.04889993, 0.01927602   )
bias1995 <- c(-0.0005531330, 0.0003773245,0.0004000032, 0.0003205394, 0.0015550907, 0.0010599529,  0.0016786352   )
unbiased1995 <- estimators1995 + bias1995
stderror1995 <- c(0.007339945,0.004237293,0.001868536, 0.003827291,0.004081243,  0.003852737,  0.002337696)

#confidence intervall
lowerb1995 <- unbiased1995 - ((1.96)*stderror1995)
upper1995 <- unbiased1995 + ((1.96)*stderror1995)
cbind(lowerb1995, upper1995)

#URBAN 2002####

#Import####
library(haven)
X21741_0001_Data <-  read_dta("C:/Users/stagiaire/Desktop/data/21741-0001-Data.dta")
#View(X21741_0001_Data)

#MALE AND FEMALE DATAFRAME####
#gender factor
X21741_0001_Data$P105 <- as.factor(X21741_0001_Data$P105) 

#create two different data sets
Genderbis <- split(X21741_0001_Data, X21741_0001_Data$P105) 

#create dataset
Data2002_male <- Genderbis$`1` 
Data2002_female <- Genderbis$`2`


#INCOME 2002#####
#create variable income
Income_2002_Female <- Data2002_female$P201

#LORENZ CURVE####

Income_2002_Female_vec <- as.vector(unlist(Income_2002_Female))
class(Income_2002_Female_vec) 

#lorenzcurve women
Income_2002_Female_naom <- na.omit(unlist(Income_2002_Female))
plot(Lc(Income_2002_Female_naom),main = "Lorenz Curve urban in 2002", xlab="Percent Population (cumulative)",
     ylab="Percent of area (%)")

#DAVIDSON GINI #####
y2<-na.omit(Income_2002_Female)

# Asymptotic SD Davidson

stdavidson2 = function(y2){
  n=length(y2)
  G=Gini(y2)
  mu=mean(y2)
  zhat=rep(0,n)
  for(i in 1:n){
    zhat[i] = -(G+1)*y2[i]+ ((2*i-1)/n)*y2[i]-2*sum(y2[1:i])/n
  }
  zbar = mean(zhat)
  varG = sum((zhat-zbar)^2)/(n*mu)^2
  return (sqrt(varG))
}

install.packages("DescTools")
library(DescTools)
gini2002 <- Gini(y2,unbiased = TRUE)  #0.3868482
segini2002 <- stdavidson2(y2) #0.006441298
#confidence interval
upperg2002 <- gini2002+ (1.96*segini2002)
lowerg2002 <- gini2002 - (1.96*segini1995)
cbind(lowerg2002,upperg2002)



#Two-Sample TEST GINI########

#H0 mu1=mu2

vargini1995 <- segini1995^2
vargini2002 <- segini2002^2

stat <- (gini1995-gini2002)/sqrt((vargini1995+vargini2002))

#rejectH0 two means are not equal

#COMPARISON 2002 1995 LORENZ CURVE #####

LCPLOT1995F<- Lc(Income_1995_naomitfemale)
LCPLOT2002F <- Lc(Income_2002_Female_naom)

plot(LCPLOT2002F,
     col="red",
     #type="b",      # !is not working
     lty=1,
     lwd=1,
     main="Lorenz Curve Women in 1995 and 2002 ",xlab="Percent Population (cumulative)",
     ylab="Percent of area (%)"     
)
lines(LCPLOT1995F, col= "brown", lty=2, lwd=1)
legend("topleft",
       c("Female 1995", "Female 2002"),
       col=c("Brown", "Red"),
       lty=c(1,1),
       lwd=3)

#CONCENTRATION CURVE####

#import the second data set for 2002
#information about the source of income 

library(haven)
DATA_2002_0003 <- read_dta("C:/Users/stagiaire/Desktop/data/21741-0003-Data.dta")
#gender as factor
DATA_2002_0003$A6 <- as.factor(DATA_2002_0003$A6) 
#Split data for men and women
genderbisbis <- split(DATA_2002_0003, DATA_2002_0003$A6) 

Data2002_female_CC <- as.data.frame(genderbisbis$`2`)

#total income Female
Income_2002_set3_F <- Data2002_female_CC$A15 

#defining new variables
Wage_urban_2002 <- Data2002_female_CC$A151
Business_urban_2002 <- Data2002_female_CC$A152
Property_urban_2002 <- Data2002_female_CC$A153
Transfer_urban_2002 <- Data2002_female_CC$A154
#create new dataset
ccurve2002 <- cbind(Income_2002_set3_F,Wage_urban_2002,Business_urban_2002,Property_urban_2002,Transfer_urban_2002)

library(IC2)
#all curves in one plot
curveConcent(x=ccurve2002[,2], y=ccurve2002[,1], col="red",lty=1,xlab="Cumulative % of population, ordered by earnings",ylab="Cumulative % of income variable")
curveConcent(x=ccurve2002[,3], y=ccurve2002[,1], col = "green",lty=1,xlab="Cumulative % of population, ordered by earnings", add= TRUE)
curveConcent(x=ccurve2002[,4], y=ccurve2002[,1], col = "blue",lty=1,xlab="Cumulative % of population, ordered by earnings", add=TRUE)
curveConcent(x=ccurve2002[,5], y=ccurve2002[,1], col = "orange",lty=1,xlab="Cumulative % of population, ordered by earnings", add=TRUE)
title(main="Concentration Curves 2002")
legend("topleft", legend = c("Wage","Business" ,"Property", "Transfer"), col=c("red", "green","blue", "orange"), lty=rep(1, 4),ncol=1)

#CONCENTRATION INDEX########

calcSConc(Wage_urban_2002,Income_2002_set3_F) #0.6014793 
calcSConc(Business_urban_2002,Income_2002_set3_F) #0.4181272 
calcSConc(Transfer_urban_2002,Income_2002_set3_F) #0.4395444
calcSConc(Property_urban_2002,Income_2002_set3_F) #0.6091928 

cc2002 <- c(0.4181272 ,0.6014793, 0.4395444, 0.6091928  )

#CORRELATION INCOME SOURCES#####

ccurve2002 <- as.data.frame(ccurve2002)
#name for column
colnames(ccurve2002) <- c("Income", "Wage", "Business", "Property", "Transfer")
#correlation
mcor2 <- cor(ccurve2002)
#new range of color
col3 <- colorRampPalette(c("white", "#6666FF", "#000033")) 
#plot corr
corrplot(mcor2, method = "number", type="upper", tl.col="black", tl.srt = 45, col=col3(20) )
#pearson test
library(Hmisc)
pvaluecorr <- rcorr(as.matrix(ccurve2002), type = "pearson")
pvaluecorr <- pvaluecorr$P

#export
library(xtable)
print(xtable(as.data.frame(pvaluecorr),type = "latex"), file = "pvaluecorr2002.tex")

#VARIABLES#####
#DESCRIPTIVE STAT#########

#create new data set for descriptive part
#putting NA value for the level "0" allow to delete the value 
#0 value are undefined value

TotalIncome_2002_set1 <- Data2002_female$P201
TotalIncome_2002_set1[TotalIncome_2002_set1 == 0] <- NA
  
Province_urban_2002 <- as.factor(Data2002_female$PROVINCE)

CPC_2002 <- Data2002_female$P110
CPC_2002[CPC_2002 == 0] <- NA
CPC_2002 <- as.factor(CPC_2002)

Education_urban_2002 <- Data2002_female$P112
Education_urban_2002[Education_urban_2002 == 0] <- NA
Education_urban_2002 <- as.factor(Education_urban_2002)

Ownership_2002 <- Data2002_female$P135
Ownership_2002[Ownership_2002 == 0] <- NA
Ownership_2002 <- as.factor(Ownership_2002) #0 obs for level 10

Occup_2002 <- Data2002_female$P141
Occup_2002[Occup_2002 == 0] <- NA
Occup_2002 <- as.factor(Occup_2002)

Industry_2002 <-  Data2002_female$P143
Industry_2002[Industry_2002 == 0] <- NA
Industry_2002 <- as.factor(Industry_2002)

#transform work exp into factor variable 

Workexp_2002 <- as.numeric(Data2002_female$P132)
Workexpfactor2002 <- cut(
  Workexp_2002,
  breaks = c(0, 11, 21, 31, 41, 100),
  labels = c("<=10", "[11:20]", "[21:30]", "[31:40]", ">=41"),
  right  = FALSE
)

#creation dataset

DATAREG2002 <- cbind(TotalIncome_2002_set1,Province_urban_2002,Education_urban_2002,CPC_2002,Workexpfactor2002,Ownership_2002,Occup_2002,Industry_2002)
DATAREG2002 <- na.omit(DATAREG2002)
DATAREG2002 <- as.data.frame(DATAREG2002)

essaie<-summary(lm(log(TotalIncome_2002_set1)~ Province_urban_2002+ CPC_2002 + Education_urban_2002 + Workexpfactor2002 + Ownership_2002 + Occup_2002 +Industry_2002, data = DATAREG2002))


#mean for each category
mean2002 <- sapply(DATAREG2002, mean, na.rm = T)
se2002 <- sapply(DATAREG2002, sd,na.rm= T)

library(Hmisc)
#freq + proportion for each category
describe(DATAREG2002[c("Province_urban_2002","CPC_2002","Education_urban_2002","Workexpfactor2002","Ownership_2002","Occup_2002","Industry_2002")])

#mean income for each category
aggprov02<- aggregate(DATAREG2002$TotalIncome_2002_set1, by=list(CPC_1995=DATAREG2002$Province_urban_2002), FUN=sum)
freqprov02<- c(245,   181,   356,   319,   203,   267,   265,   302,   109,   260,   227,   156)
aggprov02/freqprov02

aggcpc02<-aggregate(DATAREG2002$TotalIncome_2002_set1, by=list(CPC_1995=DATAREG2002$CPC_2002), FUN=sum) 
freqcpc02<- c(462,    23,   148,  2257)
aggcpc02/freqcpc02

aggedu02<-aggregate(DATAREG2002$TotalIncome_2002_set1, by=list(CPC_1995=DATAREG2002$Education_urban_2002), FUN=sum) 
freqedu02<-c( 5,     5,    91,   787,  1012,   388,   492,   107,     3)
aggedu02/freqedu02

aggwork02<-aggregate(DATAREG2002$TotalIncome_2002_set1, by=list(Workexpfactor2002=DATAREG2002$Workexpfactor2002), FUN=sum) 
freqwork02<-c(639,   973,  1039,   239)
aggwork02/freqwork02

aggown02<- aggregate(DATAREG2002$TotalIncome_2002_set1, by=list(CPC_1995=DATAREG2002$Ownership_2002), FUN=sum) 
freqown02<-c( 382,   808,   383,   191,   393,    66,    24,   133,   354,     1,   155)
aggown02/freqown02

aggocc02<- aggregate(DATAREG2002$TotalIncome_2002_set1, by=list(Occup_2002=DATAREG2002$Occup_2002), FUN=sum) 
freqocc02 <- c(12,   171,   438,    18  ,  72,   525,   456,   430,   672,    96)
aggocc02/freqocc02 #no data for farmer

aggind02<- aggregate(DATAREG2002$TotalIncome_2002_set1, by=list(CPC_1995=DATAREG2002$Industry_2002), FUN=sum) 
freqind02<- c(19,    38,  1014,   101,    92,    11,   188,   630,   102,    28,   418,    62,    52  ,  20,    32,    83)
aggind02/freqind02

mean(DATAREG2002$TotalIncome_2002_set1)

#LOGLIN REGRESSION#######
#define new variables without defining them as factor

TotalIncome_2002_set1 <- Data2002_female$P201
TotalIncome_2002_set1[TotalIncome_2002_set1 == 0] <- NA 
Province_urban_2002b <- Data2002_female$PROVINCE
CPC_2002b <- Data2002_female$P110
CPC_2002b[CPC_2002b == 0] <- NA
Education_urban_2002b <- Data2002_female$P112
Education_urban_2002b[Education_urban_2002b == 0] <- NA
Workexp_2002b <- Data2002_female$P132
Ownership_2002b <- Data2002_female$P135
Ownership_2002b[Ownership_2002b == 0] <- NA
Occup_2002b <- Data2002_female$P141
Occup_2002b[Occup_2002b == 0] <- NA
Industry_2002b <-  Data2002_female$P143
Industry_2002b[Industry_2002b == 0] <- NA

#regression

DATA2002b <- cbind(TotalIncome_2002_set1,Province_urban_2002b,CPC_2002b,Education_urban_2002b,Workexpfactor2002,Ownership_2002b,Occup_2002b,Industry_2002b)
DATA2002b <- as.data.frame(DATA2002b)
DATA2002b <- na.omit(DATA2002b)
reg2002<-summary(lm(log(TotalIncome_2002_set1)~ factor(Province_urban_2002b)+ factor(CPC_2002b) + factor(Education_urban_2002b) + factor(Workexpfactor2002) + factor(Ownership_2002b) + factor(Occup_2002b) +factor(Industry_2002b), data = DATA2002b))
reg2002

#export

library(broom)

write.csv( tidy( reg2002 ) , "coefs2002.csv" )
write.csv( glance( reg2002 ) , "an2002.csv" )


#CROSS VALIDATION####

library(tidyverse)
library(caret)
library(randomForest)

#kfolds method k=10

set.seed(123)
train.control <- trainControl(method = "cv", number=10) 

#full model

model2002full <- train(log(TotalIncome_2002_set1)~ factor(Province_urban_2002b)+ factor(CPC_2002b) + factor(Education_urban_2002b) + factor(Workexpfactor2002) + factor(Ownership_2002b) + factor(Occup_2002b) +factor(Industry_2002b),data = DATA2002b,
                 trControl = train.control, method="lm")
model2002full

#DECOMP######
library(dineq)

#fields method

ineq_decomp2002 <- dineq_rb(TotalIncome_2002_set1~ factor(Province_urban_2002b)+ 
                              factor(CPC_2002b) + factor(Education_urban_2002b) + 
                              factor(Workexpfactor2002) + factor(Ownership_2002b) +
                              factor(Occup_2002b) +factor(Industry_2002b),
                            data = DATA2002b)
#results
ineq_decomp2002["decomposition_inequality"]

#SHAPLEY##############

library(shapley)
#function to incorporate in the shapley function
reg <- function(factors, dv, data) {
  if (length(factors) == 0) return(0)
  formula <- paste(dv, "~", paste(factor(factors), collapse = "+"))
  m <- lm(formula, data = data)
  summary(m)$r.squared
}

#create new data set with ln income
Income2002shap <- log(TotalIncome_2002_set1)
DATAREG2002shapl <- cbind(Income2002shap,Province_urban_2002b,Education_urban_2002b,CPC_2002b,Workexpfactor2002,Ownership_2002b,Occup_2002b,Industry_2002b)
DATAREG2002shapl <- na.omit(DATAREG2002shapl)
DATAREG2002shapl <- as.data.frame(DATAREG2002shapl)

#results
shapley2002 <- shapley(reg, c("factor(Province_urban_2002b)", "factor(Education_urban_2002b)", "factor(CPC_2002b)", "Workexpfactor2002", "factor(Ownership_2002b)", "factor(Occup_2002b)", "factor(Industry_2002b)"), silent = TRUE, dv="Income2002shap", data= DATAREG2002shapl)
shapley2002

weight2002s <- as.vector(shapley2002$value)
weight2002s

res2002s <- 1- sum(weight2002s)
res2002s

#export
print(xtable(as.data.frame(shapley2002),type = "latex"), file = "shapley2002.tex")

#BOOTSTRAP######

#shapley function as argument in the boot function
library(shapley)
st2 = function(y, indices){
  u <- DATAREG2002shapl[indices,]
  v = shapley(reg, c("factor(Province_urban_2002b)", "factor(Education_urban_2002b)", "factor(CPC_2002b)", "Workexpfactor2002", "factor(Ownership_2002b)", "factor(Occup_2002b)", "factor(Industry_2002b)") , silent = TRUE, dv="Income2002shap",data = u)
  return(v$value)
}
#bootstrap
#100 re-sampling
library(boot)
results2002 <- boot(data=DATAREG2002shapl,st2, R=100)
#results
print(results2002)

estimators2002 <- c( 0.12872951, 0.03851449, 0.01408269, 0.03046804, 0.04954299, 0.05521837,  0.03119975     )
bias2002 <- c(0.0011812793, 0.0010133439,  0.0009989394 ,  -0.0002088821,  0.0024702202, 0.0016525289 , 0.0044995717  )
stderror2002 <- c(0.0044995717,0.005295077,0.003078696, 0.004451449,  0.005702590, 0.006684015, 0.004607354)
unbiased2002 <- estimators2002 + bias2002
lower2002 <- unbiased2002 - ((1.96)*stderror2002) 
upper2002 <- unbiased2002 + (1.96*stderror2002)
cbind(lower2002, upper2002)

#TWOSAMPLE T TEST SHAPLEY VALUE####
#variance
var1995 <- stderror1995^2
var2002 <- stderror2002^2

#t-value
tvalueshap <- unbiased1995-unbiased2002
#stat for all the variables
tvalueshap2 <- NULL
for (i in 1:7){
  tvalueshap2[i] <- (tvalueshap[i])/(sqrt((var1995[i]) + (var2002[i])))
}
print(tvalueshap2)