### R code from vignette source '/home/zeileis/svn/projects/AER/Course/Ch-Microeconometrics.Rnw' ################################################### ### code chunk number 1: setup ################################################### options(prompt = "R> ", continue = "+ ", width = 64, digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) library("AER") set.seed(1071) ################################################### ### code chunk number 2: Ch-Microeconometrics.Rnw:273-276 ################################################### data("SwissLabor", package = "AER") swiss_probit_ex <- glm(participation ~ age, data = SwissLabor, family = binomial(link = "probit")) ################################################### ### code chunk number 3: Ch-Microeconometrics.Rnw:288-289 ################################################### summary(swiss_probit_ex) ################################################### ### code chunk number 4: Ch-Microeconometrics.Rnw:303-306 ################################################### swiss_probit <- glm(participation ~ . + I(age^2), data = SwissLabor, family = binomial(link = "probit")) coeftest(swiss_probit) ################################################### ### code chunk number 5: binary-spine ################################################### plot(participation ~ age, data = SwissLabor, ylevels = 2:1) ################################################### ### code chunk number 6: Ch-Microeconometrics.Rnw:343-344 ################################################### plot(participation ~ age, data = SwissLabor, ylevels = 2:1) ################################################### ### code chunk number 7: Ch-Microeconometrics.Rnw:389-391 ################################################### fav <- mean(dnorm(predict(swiss_probit, type = "link"))) fav * coef(swiss_probit) ################################################### ### code chunk number 8: Ch-Microeconometrics.Rnw:407-412 ################################################### av <- colMeans(SwissLabor[, -c(1, 7)]) av <- data.frame(rbind(swiss = av, foreign = av), foreign = factor(c("no", "yes"))) av <- predict(swiss_probit, newdata = av, type = "link") av <- dnorm(av) ################################################### ### code chunk number 9: Ch-Microeconometrics.Rnw:417-419 ################################################### av["swiss"] * coef(swiss_probit)[-7] av["foreign"] * coef(swiss_probit)[-7] ################################################### ### code chunk number 10: Ch-Microeconometrics.Rnw:450-452 ################################################### swiss_probit0 <- update(swiss_probit, formula = . ~ 1) 1 - as.vector(logLik(swiss_probit)/logLik(swiss_probit0)) ################################################### ### code chunk number 11: Ch-Microeconometrics.Rnw:481-483 ################################################### table(true = SwissLabor$participation, pred = round(fitted(swiss_probit))) ################################################### ### code chunk number 12: Ch-Microeconometrics.Rnw:486-489 ################################################### tab <- table(true = SwissLabor$participation, pred = round(fitted(swiss_probit))) tabp <- round(100 * c(tab[1,1] + tab[2,2], tab[2,1] + tab[1,2])/sum(tab), digits = 2) ################################################### ### code chunk number 13: acc (eval = FALSE) ################################################### ## library("ROCR") ## pred <- prediction(fitted(swiss_probit), ## SwissLabor$participation) ## plot(performance(pred, "acc")) ################################################### ### code chunk number 14: acc1 ################################################### library("ROCR") pred <- prediction(fitted(swiss_probit), SwissLabor$participation) plot(performance(pred, "acc")) ################################################### ### code chunk number 15: roc (eval = FALSE) ################################################### ## plot(performance(pred, "tpr", "fpr")) ## abline(0, 1, lty = 2) ################################################### ### code chunk number 16: roc1 ################################################### plot(performance(pred, "tpr", "fpr")) abline(0, 1, lty = 2) ################################################### ### code chunk number 17: Ch-Microeconometrics.Rnw:616-619 ################################################### deviance(swiss_probit) sum(residuals(swiss_probit, type = "deviance")^2) sum(residuals(swiss_probit, type = "pearson")^2) ################################################### ### code chunk number 18: Ch-Microeconometrics.Rnw:690-694 ################################################### data("MurderRates") murder_logit <- glm(I(executions > 0) ~ time + income + noncauc + lfp + southern, data = MurderRates, family = binomial) ################################################### ### code chunk number 19: Ch-Microeconometrics.Rnw:696-697 ################################################### coeftest(murder_logit) ################################################### ### code chunk number 20: Ch-Microeconometrics.Rnw:722-726 ################################################### murder_logit2 <- glm(I(executions > 0) ~ time + income + noncauc + lfp + southern, data = MurderRates, family = binomial, control = list(epsilon = 1e-15, maxit = 50, trace = FALSE)) ################################################### ### code chunk number 21: Ch-Microeconometrics.Rnw:728-729 ################################################### coeftest(murder_logit2) ################################################### ### code chunk number 22: Ch-Microeconometrics.Rnw:755-756 ################################################### table(I(MurderRates$executions > 0), MurderRates$southern) ################################################### ### code chunk number 23: Ch-Microeconometrics.Rnw:808-811 ################################################### data("RecreationDemand") rd_pois <- glm(trips ~ ., data = RecreationDemand, family = poisson) ################################################### ### code chunk number 24: Ch-Microeconometrics.Rnw:814-815 ################################################### coeftest(rd_pois) ################################################### ### code chunk number 25: odtest1 ################################################### dispersiontest(rd_pois) ################################################### ### code chunk number 26: odtest2 ################################################### dispersiontest(rd_pois, trafo = 2) ################################################### ### code chunk number 27: Ch-Microeconometrics.Rnw:895-897 (eval = FALSE) ################################################### ## rd_qpois <- glm(trips ~ ., data = RecreationDemand, ## family = quasipoisson) ################################################### ### code chunk number 28: Ch-Microeconometrics.Rnw:943-947 ################################################### library("MASS") rd_nb <- glm.nb(trips ~ ., data = RecreationDemand) coeftest(rd_nb) logLik(rd_nb) ################################################### ### code chunk number 29: Ch-Microeconometrics.Rnw:973-975 ################################################### round(sqrt(rbind(diag(vcov(rd_pois)), diag(sandwich(rd_pois)))), digits = 3) ################################################### ### code chunk number 30: Ch-Microeconometrics.Rnw:990-991 ################################################### coeftest(rd_pois, vcov = sandwich) ################################################### ### code chunk number 31: Ch-Microeconometrics.Rnw:1017-1019 ################################################### rbind(obs = table(RecreationDemand$trips)[1:10], exp = round( sapply(0:9, function(x) sum(dpois(x, fitted(rd_pois)))))) ################################################### ### code chunk number 32: countreg-plot (eval = FALSE) ################################################### ## plot(table(RecreationDemand$trips), ylab = "") ################################################### ### code chunk number 33: countreg-plot1 ################################################### plot(table(RecreationDemand$trips), ylab = "") ################################################### ### code chunk number 34: Ch-Microeconometrics.Rnw:1096-1100 ################################################### library("pscl") rd_zinb <- zeroinfl(trips ~ . | quality + income, data = RecreationDemand, dist = "negbin") summary(rd_zinb) ################################################### ### code chunk number 35: Ch-Microeconometrics.Rnw:1111-1113 ################################################### rd_zinb_out <- capture.output(summary(rd_zinb)) writeLines(rd_zinb_out[1:18]) ################################################### ### code chunk number 36: Ch-Microeconometrics.Rnw:1123-1124 ################################################### writeLines(rd_zinb_out[-(1:18)]) ################################################### ### code chunk number 37: Ch-Microeconometrics.Rnw:1129-1130 ################################################### round(colSums(predict(rd_zinb, type = "prob")[,1:10])) ################################################### ### code chunk number 38: Ch-Microeconometrics.Rnw:1201-1204 ################################################### rd_hurdle <- hurdle(trips ~ . | quality + income, data = RecreationDemand, dist = "negbin") summary(rd_hurdle) ################################################### ### code chunk number 39: Ch-Microeconometrics.Rnw:1215-1217 ################################################### rd_hurdle_out <- capture.output(summary(rd_hurdle)) writeLines(rd_hurdle_out[1:17]) ################################################### ### code chunk number 40: Ch-Microeconometrics.Rnw:1227-1228 ################################################### writeLines(rd_hurdle_out[-(1:17)]) ################################################### ### code chunk number 41: Ch-Microeconometrics.Rnw:1234-1235 ################################################### round(colSums(predict(rd_hurdle, type = "prob")[,1:10])) ################################################### ### code chunk number 42: fairmodel ################################################### data("Affairs") aff_tob_ex <- tobit(affairs ~ yearsmarried, data = Affairs) ################################################### ### code chunk number 43: fairmodel ################################################### aff_tob <- tobit(affairs ~ age + yearsmarried + religiousness + occupation + rating, data = Affairs) ################################################### ### code chunk number 44: Ch-Microeconometrics.Rnw:1351-1352 ################################################### summary(aff_tob_ex) ################################################### ### code chunk number 45: Ch-Microeconometrics.Rnw:1367-1368 ################################################### coeftest(aff_tob) ################################################### ### code chunk number 46: Ch-Microeconometrics.Rnw:1383-1385 ################################################### aff_tob2 <- update(aff_tob, right = 4) coeftest(aff_tob2) ################################################### ### code chunk number 47: Ch-Microeconometrics.Rnw:1405-1407 ################################################### linearHypothesis(aff_tob, c("age = 0", "occupation = 0"), vcov = sandwich) ################################################### ### code chunk number 48: Ch-Microeconometrics.Rnw:1465-1466 ################################################### SwissLabor$partnum <- as.numeric(SwissLabor$participation) - 1 ################################################### ### code chunk number 49: kleinspady (eval = FALSE) ################################################### ## library("np") ## swiss_bw <- npindexbw(partnum ~ income + age + education + ## youngkids + oldkids + foreign + I(age^2), data = SwissLabor, ## method = "kleinspady", nmulti = 5) ################################################### ### code chunk number 50: kleinspady1 ################################################### if(file.exists("Ch-Microeconometrics-kleinspadybw.rda")) { library("np") load("Ch-Microeconometrics-kleinspadybw.rda") } else { library("np") swiss_bw <- npindexbw(partnum ~ income + age + education + youngkids + oldkids + foreign + I(age^2), data = SwissLabor, method = "kleinspady", nmulti = 5) save(swiss_bw, file = "Ch-Microeconometrics-kleinspadybw.rda") } ################################################### ### code chunk number 51: Ch-Microeconometrics.Rnw:1500-1501 ################################################### summary(swiss_bw) ################################################### ### code chunk number 52: Ch-Microeconometrics.Rnw:1517-1519 ################################################### swiss_ks <- npindex(bws = swiss_bw, gradients = TRUE) summary(swiss_ks) ################################################### ### code chunk number 53: Ch-Microeconometrics.Rnw:1522-1525 ################################################### swiss_ks_out <- capture.output(summary(swiss_ks)) writeLines(swiss_ks_out[1:18]) writeLines("...") ################################################### ### code chunk number 54: Ch-Microeconometrics.Rnw:1539-1541 ################################################### table(Actual = SwissLabor$participation, Predicted = round(predict(swiss_probit, type = "response"))) ################################################### ### code chunk number 55: bw-tab ################################################### data("BankWages") edcat <- factor(BankWages$education) levels(edcat)[3:10] <- rep(c("14-15", "16-18", "19-21"), c(2, 3, 3)) tab <- xtabs(~ edcat + job, data = BankWages) prop.table(tab, 1) ################################################### ### code chunk number 56: bw-plot (eval = FALSE) ################################################### ## plot(job ~ edcat, data = BankWages, off = 0) ################################################### ### code chunk number 57: bw-plot1 ################################################### plot(job ~ edcat, data = BankWages, off = 0) ################################################### ### code chunk number 58: Ch-Microeconometrics.Rnw:1657-1660 ################################################### library("nnet") bank_mnl <- multinom(job ~ education + minority, data = BankWages, subset = gender == "male", trace = FALSE) ################################################### ### code chunk number 59: Ch-Microeconometrics.Rnw:1665-1666 ################################################### coeftest(bank_mnl) ################################################### ### code chunk number 60: Ch-Microeconometrics.Rnw:1699-1703 ################################################### library("MASS") bank_polr <- polr(job ~ education + minority, data = BankWages, subset = gender == "male", Hess = TRUE) coeftest(bank_polr) ################################################### ### code chunk number 61: Ch-Microeconometrics.Rnw:1712-1714 ################################################### AIC(bank_mnl) AIC(bank_polr)