# Simultaneity bias: a small simulation ----

# Generate price and then quantity 
set.seed(1)
n <- 500
bd <- c(10, -1.0, 1.0) 
bs <- c(2, 1.5, 1.0)   
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

# OLS line is neither the demand nor the supply curve
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)

# the weight of the shocks
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"]
}

# varying sd of eps^d
sapply(c(0.3, 1, 3, 10), 
       function(sd_d) ols_slope(sd_d, sd_s = 1))   

# varying sd of eps^s
sapply(c(0.3, 1, 3, 10, 100), 
       function(sd_s) ols_slope(sd_d = 1, sd_s))   


# Examples ----

## Errors in the variables: ability measured by IQ ----
library(wooldridge)
data("wage2")
str(wage2)
fit_iq <- lm(lwage ~ educ + exper + tenure + 
               married + south + urban + black + IQ,
             data = wage2)
round(summary(fit_iq)$coefficients[c("educ", "IQ"), ], 4)

with(wage2, cor(IQ, KWW))

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

## Returns from education (having a proxy of ability) ----
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


## Returns from education (without having a proxy of ability) ----

data("SchoolingReturns", package = "ivreg")

ols_ed <- lm(log(wage) ~ education + experience + I(experience^2) +
               ethnicity + smsa + south, data = SchoolingReturns)
round(summary(ols_ed)$coefficients["education", ], 4)

## Fulton fish market
data("fish")
ols_fish <- lm(ltotqty ~ lavgprc + mon + tues + wed + thurs, data = fish)
round(summary(ols_fish)$coefficients["lavgprc", ], 4)