######################
# heteroskedasticity #
######################
# part 1
rm(list = ls())
ad = "https://www.fsb.miamioh.edu/lij14/411_smoke.txt"
data = read.table(url(ad), header=T)  
head(data)
attach(data)
summary(lm(cigs~lincome+lcigpric+educ+age+agesq+restaurn))$coef

# part 2
plot(cigs~restaurn)
sd(cigs[restaurn==0])
sd(cigs[restaurn==1])

m1 = lm(cigs~lincome+lcigpric+educ+age+agesq+restaurn)
uhat = resid(m1)
uhatsq = uhat^2
summary(lm(uhatsq~lincome+lcigpric+educ+age+agesq+restaurn))
install.packages("lmtest")
library(lmtest)
bptest(m1)

# part 3
install.packages("sandwich")
library("sandwich")
coeftest(m1, vcov = vcovHC(m1, type = "HC0"))

# part 4
luhatsq = log(uhatsq)
m3 = lm(luhatsq~lincome+lcigpric+educ+age+agesq+restaurn)
ghat = fitted(m3)
hhat = exp(ghat)
we = 1/hhat
wls_m = lm(cigs~lincome+lcigpric+educ+age+agesq+restaurn, weights=we)
summary(wls_m)