R_Script_CCTP_UBI(1).R

rm(list = ls())

install.packages("PNADcIBGE")
install.packages("dineq")
library(PNADcIBGE)
library(survey)
library(dplyr) # summary
library(dineq) # ntiles.wtd: Weighted tiles

dadosPNADc <- get_pnadc(year=2023, interview=5, defyear=2021, deflator=TRUE) 


####################
### Income type ###
####################

# Total income
summary(dadosPNADc$variables$VD4047)
TodosPSOReal <- dadosPNADc$variables$VD4047*dadosPNADc$variables$CO1e
svymean(x=~TodosPSOReal, design=dadosPNADc, na.rm=TRUE)
###
# BPC
summary(dadosPNADc$variables$V5001A)
BPCReal <- dadosPNADc$variables$V5001A2*dadosPNADc$variables$CO1e
svymean(x=~BPCReal, design=dadosPNADc, na.rm=TRUE)
# Bolsa Familia
summary(dadosPNADc$variables$V5002A)
PBFReal <- dadosPNADc$variables$V5002A2*dadosPNADc$variables$CO1e
svymean(x=~PBFReal, design=dadosPNADc, na.rm=TRUE)
# Other government social programs
summary(dadosPNADc$variables$V5003A)
OuPSReal <- dadosPNADc$variables$V5003A2*dadosPNADc$variables$CO1e
svymean(x=~OuPSReal, design=dadosPNADc, na.rm=TRUE)
# Unemployment benefit and Seguro-Defeso
summary(dadosPNADc$variables$V5005A)
SDSDReal <- dadosPNADc$variables$V5005A2*dadosPNADc$variables$CO1e
svymean(x=~SDSDReal, design=dadosPNADc, na.rm=TRUE)
# Other income
summary(dadosPNADc$variables$V5008A)
OutRReal <- dadosPNADc$variables$V5008A2*dadosPNADc$variables$CO1e
svymean(x=~OutRReal, design=dadosPNADc, na.rm=TRUE)
# Retirement Pension
summary(dadosPNADc$variables$V5004A)
AposentReal <- dadosPNADc$variables$V5004A2*dadosPNADc$variables$CO1e
svymean(x=~AposentReal, design=dadosPNADc, na.rm=TRUE)
# Alimony
summary(dadosPNADc$variables$V5006A)
PensAlReal <- dadosPNADc$variables$V5006A2*dadosPNADc$variables$CO1e
svymean(x=~PensAlReal, design=dadosPNADc, na.rm=TRUE)
# Rental or lease income
summary(dadosPNADc$variables$V5007A)
AlugReal <- dadosPNADc$variables$V5007A2*dadosPNADc$variables$CO1e
svymean(x=~AlugReal, design=dadosPNADc, na.rm=TRUE)

### Labor income
# Usual income from all sources
rendaHabtReal <- dadosPNADc$variables$VD4019*dadosPNADc$variables$CO1
svymean(x=~rendaHabtReal, design=dadosPNADc, na.rm=TRUE)
# Effective labor income
TrabEfetReal <- dadosPNADc$variables$VD4020*dadosPNADc$variables$CO3 
svymean(x=~TrabEfetReal, design=dadosPNADc, na.rm=TRUE)

### Total effective income
# Total effective income - all sources
rendaEfetReal <- dadosPNADc$variables$VD4022*dadosPNADc$variables$CO2e
svymean(x=~rendaEfetReal, design=dadosPNADc, na.rm=TRUE)
# Total effective income - all sources
rendaEfetReal <- dadosPNADc$variables$VD4022*dadosPNADc$variables$CO2e
svymean(x=~rendaEfetReal, design=dadosPNADc, na.rm=TRUE)
# Total effective income - all sources, with usual for labor (VD4046)
rendaEfetReal46 <- dadosPNADc$variables$VD4046*dadosPNADc$variables$CO3
svymean(x=~rendaEfetReal46, design=dadosPNADc, na.rm=TRUE)


####################
# Total income #
####################

### Variable for total income (real values)

## rendaHabtReal +  TodosPSOReal
Hab_PSO <- ifelse(is.na(rendaEfetReal46), NA,
                                         ifelse(is.na(rendaHabtReal),TodosPSOReal,
                                                ifelse(is.na(TodosPSOReal), rendaHabtReal, 
                                                rendaHabtReal +TodosPSOReal)))
svymean(x=~Hab_PSO, design=dadosPNADc, na.rm=TRUE)

## ... + AposentReal
Hab_PSO_Ap <- ifelse(is.na(rendaEfetReal46), NA,
                  ifelse(is.na(Hab_PSO),AposentReal,
                         ifelse(is.na(AposentReal), Hab_PSO, 
                                Hab_PSO +AposentReal)))
svymean(x=~Hab_PSO_Ap, design=dadosPNADc, na.rm=TRUE)

## ... + AlugReal
Hab_PSO_Ap_Al <- ifelse(is.na(rendaEfetReal46), NA,
                     ifelse(is.na(Hab_PSO_Ap),AlugReal,
                            ifelse(is.na(AlugReal), Hab_PSO_Ap, 
                                   Hab_PSO_Ap +AlugReal)))
svymean(x=~Hab_PSO_Ap_Al, design=dadosPNADc, na.rm=TRUE)

## ... + PensAlReal
Hab_PSO_Ap_Al_Pe <- ifelse(is.na(rendaEfetReal46), NA,
                        ifelse(is.na(Hab_PSO_Ap_Al),PensAlReal,
                               ifelse(is.na(PensAlReal), Hab_PSO_Ap_Al, 
                                      Hab_PSO_Ap_Al +PensAlReal)))
svymean(x=~Hab_PSO_Ap_Al_Pe, design=dadosPNADc, na.rm=TRUE)





####################
####  QUANTILES ####
####################

dadosPNADc$variables$Hab_PSO_Ap_Al_Pe<-Hab_PSO_Ap_Al_Pe
dadosPNADc$variables$AposentReal<-AposentReal
dadosPNADc$variables$rendaHabtReal<-rendaHabtReal
dadosPNADc$variables$BPCReal<-BPCReal
dadosPNADc$variables$PBFReal<-PBFReal
dadosPNADc$variables$OuPSReal<-OuPSReal
dadosPNADc$variables$SDSDReal<-SDSDReal
dadosPNADc$variables$OutRReal<-OutRReal
dadosPNADc$variables$PensAlReal<-PensAlReal
dadosPNADc$variables$AlugReal<-AlugReal

# Attributing zero value to NAN observations
dadosPNADc$variables$AposentReal[is.na(dadosPNADc$variables$AposentReal)]<-0
dadosPNADc$variables$rendaHabtReal[is.na(dadosPNADc$variables$rendaHabtReal)]<-0
dadosPNADc$variables$BPCReal[is.na(dadosPNADc$variables$BPCReal)]<-0
dadosPNADc$variables$PBFReal[is.na(dadosPNADc$variables$PBFReal)]<-0
dadosPNADc$variables$OuPSReal[is.na(dadosPNADc$variables$OuPSReal)]<-0
dadosPNADc$variables$SDSDReal[is.na(dadosPNADc$variables$SDSDReal)]<-0
dadosPNADc$variables$OutRReal[is.na(dadosPNADc$variables$OutRReal)]<-0
dadosPNADc$variables$PensAlReal[is.na(dadosPNADc$variables$PensAlReal)]<-0
dadosPNADc$variables$AlugReal[is.na(dadosPNADc$variables$AlugReal)]<-0

### Definition of Quantiles
quantis1 <-c(0.5,0.9,0.99,1)


### Average
# Matrix DISTR
DISTR <- matrix(0,length(quantis1),13)
colnames(DISTR) <- c("Quantil", "Lim_Quantil", "Rend_Tot", "Trab_Hab", "BPC", "Aposent", 
                     "PBF", "OuPS", "SDSD", "OutR", "PensAl", "Alug", "Trab_Hab_Pub")
QUANTI <- svyquantile(x=~Hab_PSO_Ap_Al_Pe, design=dadosPNADc, quantiles=quantis1, na.rm=TRUE) # endpoints of the interval

for(i in 1:length(quantis1)){ # For Hab_PSO_Ap_Al_Pe
  DISTR[i,1] <- quantis1[i] #Quantile
  DISTR[i,2] <- QUANTI$Hab_PSO_Ap_Al_Pe[i,1] #Total income - upper limit of the quantile
  ifelse(i== 1, 
   DISTR[i,3] <- svymean(x=~Hab_PSO_Ap_Al_Pe, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
   DISTR[i,3] <- svymean(x=~Hab_PSO_Ap_Al_Pe, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
)
}
for(i in 1:length(quantis1)){ # Usual labor income
  ifelse(i== 1, 
         DISTR[i,4] <- svymean(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,4] <- svymean(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # BPC
  ifelse(i== 1, 
         DISTR[i,5] <- svymean(x=~BPCReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,5] <- svymean(x=~BPCReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Retirement Pension
  ifelse(i== 1, 
         DISTR[i,6] <- svymean(x=~AposentReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,6] <- svymean(x=~AposentReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # BFP
  ifelse(i== 1, 
         DISTR[i,7] <- svymean(x=~PBFReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,7] <- svymean(x=~PBFReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Other government social programs
  ifelse(i== 1, 
         DISTR[i,8] <- svymean(x=~OuPSReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,8] <- svymean(x=~OuPSReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Unemployment benefit and Seguro-Defeso
  ifelse(i== 1, 
         DISTR[i,9] <- svymean(x=~SDSDReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,9] <- svymean(x=~SDSDReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Other income
  ifelse(i== 1, 
         DISTR[i,10] <- svymean(x=~OutRReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,10] <- svymean(x=~OutRReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Alimony
  ifelse(i== 1, 
         DISTR[i,11] <- svymean(x=~PensAlReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,11] <- svymean(x=~PensAlReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Rent
  ifelse(i== 1, 
         DISTR[i,12] <- svymean(x=~AlugReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTR[i,12] <- svymean(x=~AlugReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Public Sector
  ifelse(i== 1, 
         DISTR[i,13] <- svymean(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]& VD4008=="Empregado no setor público (inclusive servidor estatutário e militar)"), na.rm=TRUE),
         DISTR[i,13] <- svymean(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2] & VD4008=="Empregado no setor público (inclusive servidor estatutário e militar)"), na.rm=TRUE)
  )
}


### SUM
# Matrix DISTR
DISTRS <- matrix(0,length(quantis1),13)
colnames(DISTRS) <- c("Quantil", "Lim_Quantil", "Rend_Tot", "Trab_Hab", "BPC", "Aposent", 
    "PBF", "OuPS", "SDSD", "OutR", "PensAl", "Alug", "Trab_Hab_Pub")

for(i in 1:length(quantis1)){ # Para Hab_PSO_Ap_Al_Pe
  DISTRS[i,1] <- quantis1[i] #Quantil
  DISTRS[i,2] <- QUANTI$Hab_PSO_Ap_Al_Pe[i,1] #Total income - upper limit of the quantile
  ifelse(i== 1, 
         DISTRS[i,3] <- svytotal(x=~Hab_PSO_Ap_Al_Pe, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,3] <- svytotal(x=~Hab_PSO_Ap_Al_Pe, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Usual labor income
  ifelse(i== 1, 
         DISTRS[i,4] <- svytotal(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,4] <- svytotal(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # BPC
  ifelse(i== 1, 
         DISTRS[i,5] <- svytotal(x=~BPCReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,5] <- svytotal(x=~BPCReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Retirement pension
  ifelse(i== 1, 
         DISTRS[i,6] <- svytotal(x=~AposentReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,6] <- svytotal(x=~AposentReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # BFP
  ifelse(i== 1, 
         DISTRS[i,7] <- svytotal(x=~PBFReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,7] <- svytotal(x=~PBFReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Other government social programs
  ifelse(i== 1, 
         DISTRS[i,8] <- svytotal(x=~OuPSReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,8] <- svytotal(x=~OuPSReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Unemployment benefit and Seguro-Defeso
  ifelse(i== 1, 
         DISTRS[i,9] <- svytotal(x=~SDSDReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,9] <- svytotal(x=~SDSDReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Other income
  ifelse(i== 1, 
         DISTRS[i,10] <- svytotal(x=~OutRReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,10] <- svytotal(x=~OutRReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Alimony
  ifelse(i== 1, 
         DISTRS[i,11] <- svytotal(x=~PensAlReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,11] <- svytotal(x=~PensAlReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Rent
  ifelse(i== 1, 
         DISTRS[i,12] <- svytotal(x=~AlugReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2]), na.rm=TRUE),
         DISTRS[i,12] <- svytotal(x=~AlugReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2]), na.rm=TRUE)
  )
}
for(i in 1:length(quantis1)){ # Public Sector
  ifelse(i== 1, 
         DISTRS[i,13] <- svytotal(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & VD4008=="Empregado no setor público (inclusive servidor estatutário e militar)"), na.rm=TRUE),
         DISTRS[i,13] <- svytotal(x=~rendaHabtReal, design=subset(dadosPNADc, Hab_PSO_Ap_Al_Pe<= DISTR[i,2] & Hab_PSO_Ap_Al_Pe> DISTR[i-1,2] & VD4008=="Empregado no setor público (inclusive servidor estatutário e militar)"), na.rm=TRUE)
  )
}


### SUM
DISTRS_TOT <- matrix(0,1,13)
colnames(DISTRS_TOT) <- c("Quantil", "Lim_Quantil", "Rend_Tot", "Trab_Hab", "BPC", "Aposent", 
                          "PBF", "OuPS", "SDSD", "OutR", "PensAl", "Alug", "Trab_Hab_Pub")
  DISTRS_TOT[1] <- 1 #Quantile
  DISTRS_TOT[2] <- max(Hab_PSO_Ap_Al_Pe, na.rm = TRUE)
  DISTRS_TOT[3] <- svytotal(x=~Hab_PSO_Ap_Al_Pe, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[4] <- svytotal(x=~rendaHabtReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[5] <- svytotal(x=~BPCReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[6] <- svytotal(x=~AposentReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[7] <- svytotal(x=~PBFReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[8] <- svytotal(x=~OuPSReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[9] <- svytotal(x=~SDSDReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[10] <- svytotal(x=~OutRReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[11] <- svytotal(x=~PensAlReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[12] <- svytotal(x=~AlugReal, design=dadosPNADc, na.rm=TRUE)
  DISTRS_TOT[13] <- svytotal(x=~rendaHabtReal, design=subset(dadosPNADc, VD4008=="Empregado no setor público (inclusive servidor estatutário e militar)"), na.rm=TRUE)