library(foreign)

wide<-read.dta("C:/trait ambivalence/trimmed/bush/bush_wide.dta")
attach(wide)
hist(age,xlab='Age in Years')
box()
detach(wide)
remove(list=ls())

thermometers<-read.dta("C:/trait ambivalence/trimmed/bush/bush_long.dta")
attach(thermometers)
boxplot(bush~study,ylab='Thermometer Scores',xlab='Panel Wave')
detach(thermometers)
remove(list=ls())


correlation<-seq(0.0,0.9,by=0.1)

bias.logit.b1<-c(20.63,19.3,15.64,17,16.19,15.41,12.63,11.81,5.92,4.33)
bias.qre.b1<-c(4.13,4.05,4.18,4.40,3.74,4.51,4.01,3.9,3.36,3.19)

bias.b1<-matrix(c(bias.logit.b1,bias.qre.b1),10,2)
colnames(bias.b1)<-c("Logit", "QRE")

barplot(bias.b1, beside=TRUE, col=gray(1:10/10), legend=correlation, ylab="Percent Absolute Bias", xpd=FALSE, axis.lty=1)
box()
remove(list=ls())

states<-c(25,25,20,15,10,5)
national<-c(20,20,5,20,25,10)
interests<-matrix(c(states,national),6,2)
colnames(interests)<-c("State Groups", "National Groups")
rownames(interests)<-c("Health","Manuf.","Bank","Law","Guns","Sport")
barplot(interests,col=gray(1:6/6))
par("usr")
par(xpd=TRUE)
legend("bottom",inset=-.2,horiz=TRUE,c(rownames(interests)),fill=gray(1:6/6))
box()
remove(list=ls())


##################################################################

library(car)
rm(list=ls())

#Run a Linear Model
attach(Duncan)
duncan.model<-lm(prestige~income+education)

plot(education,duncan.model$residuals)

plot(education,duncan.model$fitted.values)



#MLE probit model
attach(Mroz)

lik.probit<-function(B,X,y){
k<-ncol(X)
beta<-B[1:k]
llik<-log(pnorm(X%*%beta))*y+log(1-pnorm(X%*%beta))*(1-y)
sum(llik)
}

lfp<-recode(lfp," 'yes'=1;'no'=0 ", as.factor=F)

indvar=cbind(1,k5,k618,age,wc,hc,lwg,inc)

probit.mroz<-optim(rep(0,8),lik.probit,hessian=TRUE,method="BFGS",control=list(maxit=2000,fnscale=-1,trace=1),X=indvar,y=lfp)

probit.mroz

OI<-solve(-1*probit.mroz$hessian)
se<-sqrt(diag(OI))
se

#GLM probit model
glm.mroz<-glm(lfp~k5+k618+age+wc+hc+lwg+inc,family=binomial(link='probit'))
summary(glm.mroz)

#Plot Predicted probabilities

ruler<-seq(30,60,.5)
age.util<-function(x){
    utility<- 1.918418-0.874712*0.2377-0.038595*1.353-0.037824*x+0.365635*1.0971-0.020525*20.129
    return(utility)
}

plot(ruler, pnorm(age.util(ruler)), ylim=c(0,1), type="l")