# Recall packages and the data 
library(ivreg)
library(wooldridge)
data("SchoolingReturns", package = "ivreg")
data("fish")
data("wage2")

fit_iq   <- lm(lwage ~ educ + exper + tenure + married + south + urban + black + IQ,
               data = wage2)

ols_ed   <- lm(log(wage) ~ education + experience + I(experience^2) +
                 ethnicity + smsa + south, data = SchoolingReturns)

ols_fish <- lm(ltotqty ~ lavgprc + mon + tues + wed + thurs, data = fish)

# Examples ----
## Wages and ability: IV regression for the measurement error ----
iv_iq <- ivreg(lwage ~ educ + exper + tenure + married + south +
                 urban + black + IQ |
                 educ + exper + tenure + married + south +
                 urban + black + KWW, data = wage2)

cf <- function(m) summary(m)$coefficients[c("educ", "IQ"), 1:2]
round(cbind(OLS = cf(fit_iq), IV = cf(iv_iq)), 4)
summary(iv_iq)

# Weak instrument
first <- lm(IQ ~ educ + exper + tenure + married + south + urban + black + KWW,
            data = wage2)
summary(first)$coefficients["KWW", ]     # t value^2 = F
summary(first)$coefficients["KWW", "t value"]^2

# Alternatively
first0 <- update(first, . ~ . - KWW)
anova(first0, first)
             
## Returns to schooling: the Wald estimator
m <-aggregate(cbind(logwage =log(wage), education)~nearcollege,data = SchoolingReturns, FUN = mean)
m

(m$logwage[2]-m$logwage[1])/(m$education[2]-m$education[1])
fs_ed <-lm(education~nearcollege+experience+ I(experience^2)+
             ethnicity+smsa+south, data = SchoolingReturns)
summary(fs_ed)

# IV regression
iv_ed <- ivreg(log(wage) ~ education + experience + I(experience^2) +
                 ethnicity + smsa + south |
                 nearcollege + experience + I(experience^2) +
                 ethnicity + smsa + south, data = SchoolingReturns)
rbind(OLS = summary(ols_ed)$coefficients["education", 1:2],
      IV  = summary(iv_ed)$coefficients["education", 1:2])
confint(iv_ed)["education", ]
round(summary(fs_ed)$coefficients["nearcollegeyes", ], 4)
summary(iv_ed)


## 2sls by hand 
X <- model.matrix(~ education + experience + I(experience^2) +
                    ethnicity + smsa + south, data = SchoolingReturns)

Z <- model.matrix(~ nearcollege + experience + I(experience^2) +
                    ethnicity + smsa + south, data = SchoolingReturns)

PZ <- Z %*% solve(t(Z)  %*% Z) %*% t(Z)

beta <- solve(t(X)%*% PZ %*% X)%*% t(X)%*% PZ %*% log(SchoolingReturns$wage)
round(beta, 3)

res <-  log(SchoolingReturns$wage) - X%*%beta

sum(res^2)/(nrow(SchoolingReturns) - 7)

# Fish 
fs_fish  <- lm(lavgprc ~ mon + tues + wed + thurs + wave2 + wave3, data = fish)
fs_fish0 <- lm(lavgprc ~ mon + tues + wed + thurs, data = fish)
anova(fs_fish0, fs_fish)[2, c("F", "Pr(>F)")]

iv_fish <- ivreg(ltotqty ~ lavgprc + mon + tues + wed + thurs |
                   mon + tues + wed + thurs + wave2 + wave3, data = fish)
rbind(OLS = summary(ols_fish)$coefficients["lavgprc", 1:2],
      IV  = summary(iv_fish)$coefficients["lavgprc", 1:2])
confint(iv_fish)["lavgprc", ]


iv_sup <- ivreg(ltotqty ~ lavgprc | mon + tues + wed + thurs, data = fish)
summary(iv_sup, diagnostics = TRUE)$diagnostics["Weak instruments", ]
round(summary(iv_sup)$coefficients["lavgprc", 1:2], 3)