#Binary Choice Example library("AER") data("SwissLabor") head(SwissLabor) summary(SwissLabor) #Run Regressions for ols(LPM), probit, and logit models #Convert part to a binary numeric variable part <- as.numeric(SwissLabor$participation)-1 #------------------------------------------------------------------------------------------------------------- #OLS regression swiss_ols <- lm(part ~ income + age + education + youngkids + oldkids + foreign + I(age^2), data=SwissLabor) summary(swiss_ols) #------------------------------------------------------------------------------------------------------------- #Probit regression using glm() swiss_probit <- glm(participation ~ . + I(age^2), data=SwissLabor, family=binomial(link = "probit")) summary(swiss_probit) #Why relationships might be nonlinear plot(participation ~ age, data=SwissLabor) plot(participation ~ education, data=SwissLabor) #Computing the marginal impact at the mean imp <- mean(dnorm(predict(swiss_probit, type="link"))) swiss_probit_int <- imp*coef(swiss_probit) swiss_probit_int #Psuedo R squared measure (McFadden's LRI) swiss_probit0 <- update(swiss_probit, formula = . ~ 1) LRI = 1 - as.vector(logLik(swiss_probit)/logLik(swiss_probit0)) LRI #------------------------------------------------------------------------------------------------------------- #Probit regression using own function foreign <- as.numeric(SwissLabor$foreign) -1 age_2 <- (SwissLabor$age)^2 xmat <- cbind(1, SwissLabor$income, SwissLabor$age, SwissLabor$education, SwissLabor$youngkids, SwissLabor$oldkids, foreign, age_2) #Probit function probit <- function(b,x,y) { P <- pnorm(x%*%b) LL <- -y%*%log(P) - (1 - y)%*%log(1-P) return(LL) } #Starting values for optimization are from OLS results b0 <- swiss_ols$coef #Initializes the optimizer my_probit <- optim(b0, probit, x=xmat, y=part, method=c("BFGS"), hessian=TRUE, control=list(maxit=10000)) #Compute standard errors, t-values, and LL value my_probit_se <- sqrt(diag(solve(my_probit$hessian))) my_logit_t <- my_probit$par/my_probit_se LL <- my_probit$value #------------------------------------------------------------------------------------------------------------- #Logit regression using glm() swiss_logit <- glm(participation ~ . + I(age^2), data=SwissLabor, family=binomial(link = "logit")) summary(swiss_logit) #------------------------------------------------------------------------------------------------------------- #Compare OLS, Probit, and Logit results b_results <- cbind(swiss_ols$coefficients, swiss_probit$coefficients, swiss_logit$coefficients) b_results #-------------------------------------------------------------------------------------------------------------