# 0. Data Generation n <- 500 x1 <- rnorm(n,sd=3) x2 <- runif(n,min=-2,max=2) def<- 6 ys <- 1 + x1 - x2 + rt(n,df=def) c <- (ys > 0) y <- pmax(ys,0) plot(1+x1-x2, c) plot(1+x1-x2, y) # 1. Binary Response Model (logit vs. probit vs. cauchit) logit <- glm(c~x1+x2,family=binomial(link="logit")) probit <- glm(c~x1+x2,family=binomial(link="probit")) cauchit <- glm(c~x1+x2,family=binomial(link="cauchit")) summary(logit) logit$coef vcov(logit) fitted(logit) plot(1 + x1 - x2,fitted(logit)) Prob <- array(0,c(n,4)) Prob[,1] <- 1-pt(-cbind(1,x1,x2)%*%c(1,1,-1),df=def) Prob[,2] <- fitted(logit) Prob[,3] <- fitted(probit) Prob[,4] <- fitted(cauchit) par(mfrow=c(2,2)) for(i in 2:4){ plot(Prob[,1],Prob[,i],xlab="True Probability",ylab="Estimated Probability") abline(0,1) } # Poisson model glm(formula,family=poisson) # 2. Sample Selection library(quantreg) # 2.1. Without any correction lm(y ~ x1 + x2) rq(y ~ x1 + x2) # 2.2. Heckman's 2 step procedure lambda <- dnorm(cbind(1,x1,x2)%*%(probit$coef))/pnorm(cbind(1,x1,x2)%*%(probit$coef)) lm(y ~ x1 + x2 + lambda) # See also "hectic" function in R package, micEcon # 2.3. Powell Estimator aa <- crq.fit.pow(cbind(1,x1,x2),y,1-c,left=TRUE) aa$coef