rm(list=ls()) library(foreign) require(sem) seed(3980123) # Create an endogenous fake data set W<-runif(2000, min=-5, max=5) # this is the instrumental variable Z<-runif(2000, min=-5, max=5) y<-c() X<-c() # y = 3 + X + 2Z + u1 # X = 4 + 2Y + W - 0.5Z + u2 # [ y -x] = [ 3 + 2Z + u1] # [ -2y x] = [ 4 - 0.5Z + W + u2] for(i in 1:2000){ mat<-matrix(data=c(1, -1, -2, 1), nrow=2, ncol=2, byrow=T) YY<-qr.solve(mat, b=c(3 + 2*Z[i] + rnorm(1, mean=0, sd=2), 4 - 0.5*Z[i] + W[i] + rnorm(1, mean=0, sd=2))) y[i]<-YY[1] X[i]<-YY[2] } dat<-data.frame(y=y, X=X, W=W, Z=Z) write.dta(dat, file="model-checking.dta") ################################################## rm(list=ls()) library(foreign) dat <- read.dta("model-checking.dta") attach(dat) # two standard OLS models mod.1<-lm(y~X) mod.2<-lm(y~X+Z) AIC(mod.1) AIC(mod.2) BIC(mod.1) BIC(mod.2) # dffits and dfbetas mod.2.dfb <- dfbetas(mod.2) head(mod.2.dfb) plot(mod.2.dfb[,1] ~ X) plot(mod.2.dfb[,1] ~ Z) mod.2.dff <- dffits(mod.2) head(mod.2.dff) plot(mod.2.dff ~ X) plot(mod.2.dff ~ Z) # 2SLS model library(sem) mod.iv.2 <- tsls(y~X+Z, ~Z+W) # Hausman test from http://www.klein.co.uk/downloads/MPhil/Lent/LabSessionLent3_R.R cf_diff <- coef(mod.iv.2) - coef(mod.2) vc_diff <- vcov(mod.iv.2) - vcov(mod.2) x2_diff <- as.vector(t(cf_diff) %*% solve(vc_diff) %*% cf_diff) # is Chi^2 with degrees of freedom (df) = rank of matrix vc_diff pchisq(q = x2_diff, df = 2, lower.tail = FALSE) detach(dat) # model fit quality? rm(list=ls()) library(foreign) dat2 <-read.dta("505-homework3-data.dta") x <- seq(from=8, to=10, by=0.1) cauch <- dcauchy(x) prob <- dnorm(x) logist <- dlogis(x) plot(cauch ~ x, type="l", lwd=3, col="darkgray", ylim=c(0.00, 0.006)) lines(prob ~ x, lty = 2) lines(logist ~ x, lty=3) legend("topleft", legend=c("cauchy", "probit/normal","logistic"), lty=c(1, 2, 3), lwd=c(3,1,1), col=c("darkgray","black","black")) mod.logit <- glm(y~x+z, data=dat2, family=binomial(link="logit")) mod.probit <- glm(y~x+z, data=dat2, family=binomial(link="probit")) mod.cauchit <- glm(y~x+z, data=dat2, family=binomial(link="cauchit")) library(pscl) vuong(mod.logit,mod.probit) vuong(mod.logit,mod.cauchit) BIC(mod.logit) BIC(mod.probit) BIC(mod.cauchit) AIC(mod.logit) AIC(mod.probit) AIC(mod.cauchit) # dffits and dfbetas mod.2.dfb <- dfbetas(mod.cauchit) head(mod.2.dfb) plot(mod.2.dfb[,1] ~ dat2$x) plot(mod.2.dfb[,1] ~ dat2$z) mod.2.dff <- dffits(mod.cauchit) head(mod.2.dff) plot(mod.2.dff ~ dat2$x) plot(mod.2.dff ~ dat2$z)