#Download female wage data and make some adjustments

data("PSID1976", package="AER")

PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0,

levels = c(FALSE, TRUE), labels = c("no", "yes")))

PSID1976$nwincome <- with(PSID1976, (fincome - hours * earnings)/1000)

PSID1976$partnum <- as.numeric(PSID1976$participation) - 1

PSID1976$lhours <- log(PSID1976$hours+1)

#---------------------------------------------------------------------------------------------------------------

Plot histogram of hours and log(1+hours) next to eachother in single figure

par(mfrow=c(1,2))

hist(PSID1976$hours, main="Hours")

hist(PSID1976$lhours, main="log(1+Hours)")

#---------------------------------------------------------------------------------------------------------------

#Run OLS and Tobit model and compare

library(AER)

form <- lhours ~ nwincome + education + experience + I(experience^2) + age + youngkids + oldkids 

wl_ols <- lm(form, data = PSID1976)

wl_tobit <- tobit(form, data = PSID1976)

summary(wl_ols)

summary(wl_tobit)



mult <- mean(pnorm(predict(wl_tobit)/wl_tobit$scale))

wl_tobit_slope <- mult*wl_tobit$coefficients



cbind(wl_ols$coef,wl_tobit_slope)

#---------------------------------------------------------------------------------------------------------------

#Run Heckman sample selection model

library("sampleSelection")



wl_heckman <- heckit(partnum ~ nwincome + education + experience + I(experience^2) + age + youngkids + oldkids,

	lhours ~ nwincome + education + experience + I(experience^2) + age + youngkids + oldkids , 

	PSID1976, method="2step")



k <- (length(wl_heckman$coef)-3)/2

cbind(round(wl_ols$coef,digits=5),round(wl_tobit_slope,digits=5),round(wl_heckman$coef[(k+1):(2*k)],digits=5)) 

#---------------------------------------------------------------------------------------------------------------