#############################################################################
#
# PROBIT MODEL
#
#############################################################################

# Married Women’s Annual Labor Supply
#
# Wooldridge (5th edition, Chapter 7, Example 17.2)
#
# inlf  (“in the labor force”) be a binary variable indicating
# labor force participation by a married woman during 1975: 
# inlf 1 if the woman reports working for a wage outside the 
# home at some point during the year, and zero otherwise.
#
# It is assumed that labor force participation depends on other 
# sources of income, including husband’s earnings (nwifeinc , 
# measured in thousands of dollars), years of education (educ ),
# past years of labor market experience (exper ), age , 
# number of children less than six years old (kidslt6 ), 
# and number of kids between 6 and 18 years of age (kidsge6 ). 
# Using the data in MROZ.RAW from Mroz (1987) 428 of the 753 
# women in the sample report being in the labor force at some 
# point during 1975.
#
# Obs:   753
#
#  1. inlf          =1 if in labor force, 1975
#  2. hours         hours worked, 1975
#  3. kidslt6       # kids < 6 years
#  4. kidsge6       # kids 6-18
#  5. age           woman's age in yrs
#  6. educ          years of schooling
#  7. wage          estimated wage from earns., hours
#  8. repwage       reported wage at interview in 1976
#  9. hushrs        hours worked by husband, 1975
#  10. husage        husband's age
#  11. huseduc      husband's years of schooling
#  12. huswage      husband's hourly wage, 1975
#  13. faminc       family income, 1975
#  14. mtr          fed. marginal tax rate facing woman
#  15. motheduc     mother's years of schooling
#  16. fatheduc     father's years of schooling
#  17. unem         unem. rate in county of resid.
#  18. city         =1 if live in SMSA
#  19. exper        actual labor mkt exper
#  20. nwifeinc     (faminc - wage*hours)/1000
#  21. lwage        log(wage)
#  22. expersq      exper^2
#
# Mroz (1987) The Sensitivity of an Empirical Model of Married 
# Women’s Hours of Work to Economic and Statistical Assumptions, 
# Econometrica, 55, 765–799.
#
#############################################################################
rm(list=ls())

pdf(file="Probit-model.pdf",width=12,height=8)

data = read.table("MROZ.txt",header=TRUE)
attach(data)
n = nrow(data)

####################################################
# Bayesian inference when y:hours and x:exper
####################################################
x     = exper
y     = inlf
xs = sort(unique(x))
nx = length(xs)
count = matrix(0,nx,2)
for (i in 1:nx){
  count[i,1] = sum(x[y==0]==xs[i])
  count[i,2] = sum(x[y==1]==xs[i])
}

# MCMC setup
M0    = 10000
M     = 2000
L     = 10
niter = M0+L*M

# Bayesian probit model
X     = cbind(1,x)
XtX   = t(X)%*%X
iXtX  = solve(XtX)
Xty   = t(X)%*%y
tciXtX= t(chol(iXtX))
A     = iXtX%*%t(X)

# Initial values
b   = lm(y~x)$coef
bs  = matrix(0,niter,2)
y11 = rep(0,n)
set.seed(2352)
for (iter in 1:niter){
  m = b[1]+b[2]*x
  u = runif(n)
  y11[y==0] = qnorm(u[y==0]*pnorm(0,m[y==0],1),m[y==0],1)
  y11[y==1] = qnorm(1-u[y==1]*(1-pnorm(0,m[y==1],1)),m[y==1],1)
  b = A%*%y11 + tciXtX%*%rnorm(2)
  bs[iter,] = b
}
bs = bs[seq(M0+1,niter,by=L),]

par(mfrow=c(2,3))
ts.plot(bs[,1],xlab="iteration",ylab="",main=expression(beta[0]))
acf(bs[,1],main="")
hist(bs[,1],xlab="",main="",prob=TRUE)
ts.plot(bs[,2],xlab="iteration",ylab="",main=expression(beta[1]))
acf(bs[,2],main="")
hist(bs[,2],xlab="",main="",prob=TRUE)

N   = 100
xxx = seq(min(x),max(x),length=N)
prob = matrix(0,N,M)
for (i in 1:N)
  prob[i,] = pnorm(-bs[,1]+bs[,2]*xxx[i])
qprob = apply(prob,1,quantile,c(0.025,0.5,0.975))

par(mfrow=c(1,1))
plot(xxx,qprob[2,],ylab="#Prob of being in labor force, 1975",
     xlab="actual labor mkt exper",ylim=c(0.4,1.1),col=2,
     type="l",lwd=2,axes=FALSE)
lines(xxx,qprob[1,],col=2,lwd=2,lty=2)
lines(xxx,qprob[3,],col=2,lwd=2,lty=2)
axis(2);box();axis(1,at=xs)
for (i in 1:nx){
  text(xs[i],0.4,count[i,1],cex=0.95)
  text(xs[i],1.1,count[i,2],cex=0.95)
}

par(mfrow=c(1,1))
plot(density(prob[6,]),xlim=c(0.4,1),ylim=c(0,15),main="",
     xlab="Prob of being in labor force, 1975",lwd=2)
lines(density(prob[11,]),col=2,lwd=2)
lines(density(prob[16,]),col=3,lwd=2)
lines(density(prob[21,]),col=4,lwd=2)
lines(density(prob[26,]),col=5,lwd=2)
lines(density(prob[31,]),col=6,lwd=2)
legend(0.4,12,legend=c("5y","10y","15y","20y","25y","30y"),
       col=1:6,lwd=2,bty="n")
text(0.45,12,"actual labor mkt exper") 

dev.off()