###################################################
### chunk number 1: resdens
###################################################
#line 126 "Lecture10_2010.Rnw"
library(sm)
Weakliem <- read.table("Weakliem.txt")
mod1 <- lm(secpay~gini + democrat, data=Weakliem)
sm.density(rstudent(mod1), model="normal")


###################################################
### chunk number 2: qq1
###################################################
#line 151 "Lecture10_2010.Rnw"
library(car)
qq.plot(mod1)


###################################################
### chunk number 3: qq2
###################################################
#line 165 "Lecture10_2010.Rnw"
mod2 <- lm(secpay ~ gini + democrat, data=Weakliem,
    subset=-c(25,49))
qq.plot(mod2, simulate=T, labels=FALSE)


###################################################
### chunk number 4: resdens2
###################################################
#line 170 "Lecture10_2010.Rnw"
sm.density(rstudent(mod2), model="normal")


###################################################
### chunk number 5: assesshet1
###################################################
#line 240 "Lecture10_2010.Rnw"
plot(fitted.values(mod1), rstudent(mod1),
    main="Studentized Residuals versus Fitted Values")


###################################################
### chunk number 6: remout
###################################################
#line 249 "Lecture10_2010.Rnw"
infl.outliers <- which(rownames(Weakliem) %in% 
    c("Slovakia", "CzechRepublic"))
mod2 <- lm(secpay~gini*democrat, data=Weakliem,
    subset=-infl.outliers)

plot(fitted.values(mod2), rstudent(mod2),
    main="Studentized Residuals versus Fitted Values")
abline(h=0, lty=2)


###################################################
### chunk number 7: dhs
###################################################
#line 278 "Lecture10_2010.Rnw"
library(apsrtable)
library(foreign)
dat <- read.dta("dhs_sl.dta")
mod <- lm(ceb ~ age + christian + educ + notv, data=dat)
apsrtable(mod, model.names="", Sweave=T, digits=3)


###################################################
### chunk number 8: rvf1
###################################################
#line 301 "Lecture10_2010.Rnw"
plot(fitted.values(mod), rstudent(mod),
	main="Studentized Residuals vs Fitted Values")
abline(h=0, lty=2)


###################################################
### chunk number 9: spreadlevel1
###################################################
#line 306 "Lecture10_2010.Rnw"
spread.level.plot(mod)


###################################################
### chunk number 10: sl2
###################################################
#line 309 "Lecture10_2010.Rnw"
sl2 <- spread.level.plot(mod)


###################################################
### chunk number 11: ncv
###################################################
#line 365 "Lecture10_2010.Rnw"
ncv.test(mod, ~age + christian + educ + notv, data=dat)


###################################################
### chunk number 12: wls
###################################################
#line 418 "Lecture10_2010.Rnw"
mod.wls <- lm(ceb ~ age + christian + educ + notv, data=dat, 
     weight=1/age)
with(summary(mod.wls), coefficients)


###################################################
### chunk number 13: fgls
###################################################
#line 477 "Lecture10_2010.Rnw"
mod1.ols <- lm(ceb ~ age + christian + educ + notv, data=dat)
aux.mod1 <- lm(log(resid(mod1.ols)^2) ~ age + christian + educ + notv, data=dat)
h <- exp(predict(aux.mod1))
mod.fgls <- lm(ceb ~ age + christian + educ + notv, data=dat, weights=1/h)
with(summary(mod.fgls), coefficients)


###################################################
### chunk number 14: hetreg
###################################################
#line 517 "Lecture10_2010.Rnw"
source("~/R/MLhetreg.R")
X <- model.matrix(mod1.ols)[,-1]
Z <- X
het.mod <- MLhetreg(dat$ceb, X, Z, )
summary(het.mod)


###################################################
### chunk number 15: compmod
###################################################
#line 531 "Lecture10_2010.Rnw"
b1 <- coef(mod1.ols)
b2 <- coef(mod.wls)
b3 <- coef(mod.fgls)
b4 <- het.mod$par[1:length(b1)]
v1 <- sqrt(diag(vcov(mod1.ols)))
v2 <- sqrt(diag(vcov(mod.wls)))
v3 <- sqrt(diag(vcov(mod.fgls)))
v4 <- sqrt(diag(solve(het.mod$hessian)))[1:length(b1)]
coefs <- cbind(b1,b2,b3,b4)
coefs <- apply(coefs, 2, function(x)sprintf("%.3f", x))
vars <- cbind(v1, v2, v3, v4)
vars <- apply(vars, 2, function(x)paste("(", sprintf("%.3f", x), ")", sep=""))
out <- do.call(rbind, lapply(1:nrow(coefs), function(x)rbind(coefs[x, ], vars[x, ])))
rownames(out) <- sapply(1:nrow(out), function(x)paste(rep(" ", x), collapse=""))
rownames(out)[seq(1,nrow(out), by=2)] <- names(b1)
colnames(out) <- c("OLS", "WLS", "FGLS", "HetReg")


###################################################
### chunk number 16: compmod_tab
###################################################
#line 554 "Lecture10_2010.Rnw"
library(xtable)
print(xtable(out), only.contents=T)


###################################################
### chunk number 17: varres
###################################################
#line 568 "Lecture10_2010.Rnw"
gamma.num <- het.mod$par[(length(b1)+1):length(het.mod$par)]
gamma.var.num <- sqrt(diag(solve(het.mod$hessian)))[(length(b1)+1):length(het.mod$par)]

gamma <-sprintf("%.3f", gamma.num)
gamma.var <- paste("(", sprintf("%.3f", gamma.var.num), ")", sep="")
var.mat <- matrix(sapply(1:length(gamma), function(x)c(gamma[x], gamma.var[x])), ncol=1)
rownames(var.mat) <- sapply(1:nrow(out), function(x)paste(rep(" ", x), collapse=""))
rownames(var.mat)[seq(1,nrow(out), by=2)] <- names(b1)


###################################################
### chunk number 18: print_varres
###################################################
#line 582 "Lecture10_2010.Rnw"
print(xtable(var.mat), only.contents=T, include.colnames=F)


###################################################
### chunk number 19: varfig
###################################################
#line 588 "Lecture10_2010.Rnw"
newdat <- data.frame(
	age = seq(min(dat$age), max(dat$age), length=100), 
	christian = 1, 
	notv = 0, 
	educ = factor(1, levels=1:length(levels(dat$educ)), labels=levels(dat$educ)))
Z <- model.matrix(as.formula("~age + christian + educ + notv"), data=newdat)
eff <- exp(Z %*% gamma.num)
pdf("varplot.pdf", height=6, width=6)
plot(newdat$age, eff, type="l", xlab="Age", ylab= "Predicted Residual Variance")
lines(newdat$age, sqrt(eff), lty=2)
legend("topleft", c("Variance", "Standard Deviation"), lty=c(1,2), inset=.01)
invisible(dev.off())


###################################################
### chunk number 20: robse.fun
###################################################
#line 723 "Lecture10_2010.Rnw"
robust.se <- function( model, type="hc3") {
  require(car)
  s <- summary( model)
  wse <- sqrt( diag( hccm( model, type=type)))
  t <- model$coefficients/wse
  p <- 2*pt(abs( t), model$df.residual, lower.tail=F)
  results <- round(cbind( model$coefficients, wse, t, p), 3)
  dimnames(results) <- dimnames( s$coefficients)
  results
}


###################################################
### chunk number 21: coef1
###################################################
#line 741 "Lecture10_2010.Rnw"
with(summary(mod), coefficients)


###################################################
### chunk number 22: coef2
###################################################
#line 749 "Lecture10_2010.Rnw"
robust.se(mod)


###################################################
### chunk number 23: rep1
###################################################
#line 854 "Lecture10_2010.Rnw"
## Generate Data
set.seed(1234)
n <- 100
x <- matrix(rlnorm(2*n), ncol=2)
y <- 0 + 0*x[,1] + 0*x[,2] + x[,1]^2*rnorm(100)

## Run Models
mod <- lm(y ~ x[,1] + x[,2])
restricted.mod <- lm(y ~ x[,2])

## Save Residuals
u.hat <- residuals(restricted.mod)

## Obtain W and other model quantities
W <- cbind(1, x, 1/x, x[,1]/x[,2], 
	x[,2]/x[,1], 1/(x[,1]*x[,2]), 
	x[,1]*x[,2], x^2)
X <- model.matrix(mod)
a <- 1/(1-hatvalues(restricted.mod))
omega.hat <- diag(a^2*u.hat^2)

b.cragg <- solve(t(X)%*%W%*%solve(t(W)%*%
	omega.hat%*%W)%*%t(W)%*%X)%*%t(X)%*%W%*%
	solve(t(W)%*%omega.hat%*%W)%*%t(W)%*%y

v.cragg <- solve(t(X)%*%W%*%solve(t(W)%*%
	omega.hat%*%W)%*%t(W)%*%X)

t1 <- as.vector(b.cragg[2]/sqrt(diag(v.cragg))[2])


###################################################
### chunk number 24: wildboot
###################################################
#line 895 "Lecture10_2010.Rnw"
res <- NULL
for(i in 1:1000){
y.star <- restricted.mod$fitted + 
	restricted.mod$residuals * a * 
	sample(c(-1,1), n, replace=T)
b.cragg2 <- solve(t(X)%*%W%*%solve(t(W)%*%omega.hat%*%
	W)%*%t(W)%*%X)%*%t(X)%*%W%*%solve(t(W)%*%omega.hat%*%
	W)%*%t(W)%*%y.star
v.cragg2 <- solve(t(X)%*%W%*%solve(t(W)%*%omega.hat%*%W)%*%
	t(W)%*%X)
t.star <- b.cragg2[2]/sqrt(diag(v.cragg2))[2]
res <- c(res, t.star)
}


###################################################
### chunk number 25: pval
###################################################
#line 913 "Lecture10_2010.Rnw"
mean(res > t1)