# Lecture 7: Endogeneity

# True model: $$y^*_i = 1 + 2x^*_i + \varepsilon^*_i$$ 
# $x^*_i \sim \mathcal{N}(0,1)$ 
# $\varepsilon^{*}_i \sim \mathcal{N}(0,1)$ 
# measurement error $\nu_i \sim \mathcal{N}(0, 0.5)$ (the variance is $0.5$)

# A small simulation
set.seed(1)
beta_1 <- 1 
beta_2 <- 2 
s2x <- 1 
s2nu <- 0.5 
nobs <- 200
B <- 1000
est <- matrix(NA, B, 3, 
              dimnames = list(NULL, c("no err", "err y", "err x")))
r2  <- matrix(NA, B, 2, 
              dimnames = list(NULL, c("no err", "err y")))


for (b in 1:B) {
  xs    <- rnorm(nobs, 0, sqrt(s2x))
  ys    <- beta_1 + beta_2 * xs + rnorm(nobs)
  y_obs <- ys + rnorm(nobs, 0, sqrt(s2nu))
  x_obs <- xs + rnorm(nobs, 0, sqrt(s2nu))
  est[b, ] <- c(coef(lm(ys ~ xs))[2], 
                coef(lm(y_obs ~ xs))[2],
                coef(lm(ys ~ x_obs))[2])
  r2[b, ]  <- c(summary(lm(ys ~ xs))$r.squared,
                summary(lm(y_obs ~ xs))$r.squared)
}

colMeans(est)                         
beta_2 * (1 - s2nu / (s2x + s2nu))     
colMeans(r2)         


boxplot(est, ylab = expression(hat(beta)[2]))
abline(h = beta_2, col = "red")


est <- rep(NA, B)
size <- c(50, 500, 5000)
for (nn in 1:length(size)) {
  for(b in 1:B){
    xs <- rnorm(size[nn])
    x_obs <- xs + rnorm(size[nn], 0, sqrt(s2nu))
    ys <- beta_1 + beta_2 * xs + rnorm(size[nn])
    est[b] <- coef(lm(ys ~ x_obs))[2]
  } 
  print(paste0("n = ", size[nn], " mean of hat_b2 = ", round(mean(est), 3))) 
} 


# Let's do togheter

# True model: $$y_i = 1 + 2x_i + 3z_i + \varepsilon_i$$
# $z_i, \varepsilon_i, u_i \sim \mathcal{N}(0,1)$, independent of each other
# $x_i = \rho z_i + \sqrt{1-\rho^2}\,u_i$, with $\rho = 0.5$ or $\rho = 0$. 
# Then $\text{var}(x) = \rho^2 + (1-\rho^2) = 1$ and $\text{cov}(x,z) = \rho\,\text{var}(z) = \rho$
# Three fits: 
#  1.  full model with $x$ and $z$; 
#  2. $z$ omitted and $\text{cov}(x,z)=0.5$; 
#  3. $z$ omitted and $\text{cov}(x,z)=0$


# Econometric example
data <- read.table("kpt.dat")
colnames(data) <- c("wage", "educ", "AFQT", "educSibl", "exper",
                    "tenure", "educM", "educF", "urban", "home")

data$exper <- data$exper/52

fit1 <- lm(wage ~ educ + poly(exper, 2, raw = T), data = data)
summary(fit1)$coefficients


with(data, cor(educ, AFQT))


fit2 <- lm(wage ~ educ + poly(exper, 2, raw = T) + AFQT, data = data)
summary(fit2)$coefficients


# Simultaneity bias 
set.seed(1)
n <- 500
bd <- c(10, -1.0, 1.0); bs <- c(2, 1.5, 1.0)     # (intercept, price, shifter)
d <- rnorm(n); s <- rnorm(n); ed <- rnorm(n); es <- rnorm(n)
p <- ((bd[1] - bs[1]) + bd[3]*d - bs[3]*s + ed - es) / (bs[2] - bd[2])
q <- bs[1] + bs[2]*p + bs[3]*s + es              # equilibrium quantity

plot(p, q, col = "grey60", xlab = "price", ylab = "quantity")
abline(bd[1], bd[2], col = "blue", lwd = 2)      # demand (d = 0, no shock)
abline(bs[1], bs[2], col = "red",  lwd = 2)      # supply (s = 0, no shock)
abline(lm(q ~ p), lwd = 2, lty = 2)              # OLS
legend("topleft", c("demand", "supply", "OLS"), lwd = 2, lty = c(1, 1, 2),
       col = c("blue", "red", "black"), bty = "n", cex = 0.7)

coef(lm(q ~ p))["p"]       # true demand slope: -1.0; true supply slope: +1.5
cor(p, ed); cor(p, es)


ols_slope <- function(sd_d, sd_s, n = 1e5) {
  d <- rnorm(n)
  s <- rnorm(n)
  ed <- rnorm(n, 0, sd_d)
  es <- rnorm(n, 0, sd_s)
  p <- ((bd[1] - bs[1]) + bd[3]*d - bs[3]*s + ed - es) / (bs[2] - bd[2])
  q <- bs[1] + bs[2]*p + bs[3]*s + es
  coef(lm(q ~ p))["p"]
}
sapply(c(0.3, 1, 3), function(sd_d) ols_slope(sd_d, 1))   # varying sd of eps^d