## a.a. 2023--2024
## Esercitazione aula 
## Data 
## Data Documentation on cps2014class CPS data set (created 10/1/2015, Barry Hirsch)
## The data set cps2014class (available in SAS and Stata formats) includes 137,374 non-student wage and salary workers, ages 18 to 65, in 2014, extracted from the Census Bureau's Current Population Survey Outgoing Rotation Group files (i.e., the 1⁄4 sample asked earnings and union questions each month). Below is a brief description of the variables you will need to work with in your take-home exercise (a small number of other variables are included, but not defined below). Variable names in Stata dataset are in lower case.
#AGE: 18 to 65            
#STATE: 
#SCHOOL: Years of schooling completed (0-20) 
#FEMALE: 1=females, 0=male
#BLACK: 1=black, 0=other race  
#HISPANIC: 1=Hispanic; 0=other ethnicity
#FORBORN: 1=foreign born (includes citizens and non-citizens); 0=native born
#MARRIED: 1=married spouse present; 0=ever married, spouse not present; never married
#UNION: 1=union member; 0=not a member
#WKEARN: Usual weekly earnings ($1 to estimated mean above the topcode)
#IMPUTE: 1=earnings non-respondent, earnings imputed; 0=earnings respondent, not imputed
#WKHOURS: Usual hours worked per week (1-99)
#HRSLSTWK: Hours worked last per week (1-99)
#PARTTIME: 1=less than 35 usual hours worked per week; 0=35+ usual hours
#MSASIZE: . =not in MSA or small MSA not identified; 2=100t-250t, 3=250t-500t, 4=500t-1m, 5=1m-2.5m, 6=2.5m-5m, 7=5m+
#State codes (11-95); regional breakdown shown below
#Northeast
#Middle Atlantic East North Central West North Central South Atlantic
#East South Central West South Central Mountain
#Pacific
#10<state<20 20<state<30 30<state<40 40<state<50 50<state<60 60<state<70 70<state<80 80<state<90 90<state<100

## Importo l funzioni
library(arm)
LeUni <- function(vec) length(unique(vec))
sumNA<-function(vec){
	
sum(is.na(vec))
}


##Leggo i dati
salaries <- readstata13::read.dta13("/Users/mariagraziapittau/Library/CloudStorage/Dropbox-Universita'DegliStudiDiRomaLaSapienza/grazia.pittau@uniroma1.it/Lezioni/Statistica_Aziendale/R/Regressione/cps2014class.dta", generate.factors=T, nonint.factors=T) 

salaries[1:5,]
apply(salaries, 2, sumNA)

## Salari orari
salaries[, "hrly_earn"]<-salaries[, "wkearn"]/salaries[, "wkhours"]


##### Range e media del salario
round(range(salaries[, "hrly_earn"]), 2) ## 3.00--297 $orari
round(mean(salaries[, "hrly_earn"]), 2)
## 24 dollari l'ora 

##medie condizionate
prop.table(table(salaries[, "female"]))

round(mean(salaries[, "hrly_earn"][salaries$female==1]), 2) ## 21.55 $/h
round(mean(salaries[, "hrly_earn"][salaries$female==0]), 2) ## 27.20 $/h

## Figure corrispondenti
par(mfrow=c(2,2))
hist(salaries[, "hrly_earn"], freq = F, col = "gray", main = "", xlab = "Salario", ylab = "Densità")  
lines(density(salaries[, "hrly_earn"]), col="red") 

## in log
hist(log(salaries[, "hrly_earn"]), freq = F, col = "gray", main = "", xlab = "Salario", ylab = "Densità")  
lines(density(log(salaries[, "hrly_earn"])), col="red") 

dev.off()

par(mfrow = c(2, 2))
plot(jitter(salaries$female, 0.25), salaries[, "hrly_earn"], cex = 0.8, xaxt = "n", xlim = c(-0.5, 1.5), xlab = "Genere", ylab = "Salario") 
axis(1, c(0, 1), labels = c("U", "D")) 
text(1.2, 200, paste("Media=", mean(round(mean(salaries[, "hrly_earn"][salaries$female==1]), 2)), "$/h"))
text(0.3, 200, paste("Media=", mean(round(mean(salaries[, "hrly_earn"][salaries$female==0]), 2)), "$/h"))

salario.m <- density(salaries[, "hrly_earn"][salaries$female==0])
salario.d <- density(salaries[, "hrly_earn"][salaries$female==1])

plot(salario.m, main = "Distribuzione dei salari orari per sesso", xlab = "Salari (dollari l'ora)", ylab = "Densità", col="blue", ylim=c(0, 1)*range(range(salario.d[[2]]),range(salario.m[[2]])))
lines(salario.d, main = "", xlab = "Salari", ylab = "Densità", col="red")
text(200, 0.02, "donne", col="red")
text(200, 0.025, "uomini", col="blue")


##relazioni con le altre variabili
## years of potential experience measured by [age-schooling-6] or [age-18], whichever is smaller (if any
#experience values are negative, set them to zero)

salaries[, "experience"]<-salaries[, "age"]-salaries[, "school"]-6
range(salaries[, "experience"])
salaries[, "experience"]<-ifelse(salaries[, "experience"]<=0, 0, salaries[, "experience"])


plot(salaries[, "school"], salaries[, "hrly_earn"], xlab="Anni di istruzione", ylab = "Salario (dollari/ora)", cex.lab = 1.2, pch = 19, 
             col = ifelse(salaries[, "female"] == 0, "blue", "red"), cex=0.3)
legend(7,250, paste(c("donne", "uomini")), col=c("red", "blue"), pch = 16, cex = 0.9)


plot(salaries[, "experience"], salaries[, "hrly_earn"], xlab="Anni di esperienza", ylab = "Salario (dollari/ora)", cex.lab = 1.2, pch = 19, 
             col = ifelse(salaries[, "female"] == 0, "blue", "red"), cex=0.3)
legend(7,250, paste(c("donne", "uomini")), col=c("red", "blue"), pch = 16, cex = 0.9)

## Stimate un'equazione dei salari in cui il logaritmo dei salari orari
# (calcolato come il salario medio settimanale diviso il numero di ore lavorate mediamente nella settimana) 
#viene messo in relazione con gli anni di scuola, gli anni di potenziale esperienza lavorativa, 
# con l'esperienza misurata come differenza tra [età anagrafica - anni di scuola - 6], 
# se negativa impone che sia pari a 0; utilizzate come variabili di controllo almeno il genere, l'etnia e l'iscrizione al sindacato. 

## Regressione salari~ sesso
z.exp<-rescale(salaries[, "experience"])
z.exp2<-rescale(salaries[, "experience"])^2

## Punto 1 Intero dataset
fit.1 <- lm(log(hrly_earn) ~ rescale(school) +rescale(experience)+ I(rescale(experience)^2) +

female + black + hispanic + union, data = salaries)
display(fit.1)

## Punto 2 solo dati non imputati (parlare di subset)
##Now estimate the exact same regression, but omit all imputed earners.
fit.2 <- lm(log(hrly_earn) ~ rescale(school) + rescale(experience) + I(rescale(experience)^2) + 
female + black + hispanic + union, data = salaries, subset=impute==0)
display(fit.2)

## Punto 3 da fare in aula (ragazzi)