rm(list = ls())       # delete everything from the memory

library(AER)        # package with IV regression 
library(stargazer)

n <- 1000

x1. <- rnorm(n,0,2)
e <- 2*x1. + rnorm(n, 0, 3) 

z <- rnorm(n,0,5)
x1 <- x1. + z

x2 <- rchisq(n,2)

y <- 5 + 2*x1 + 3*x2 + e 


mean(x1)
var(x1)
sd(x1)


cov(x1,e)
cov(x1,e) / (sd(x1) * sd(e))

cor(x1,e)
cor(x2,e)
cor(z,e)


lm(y ~ x1 + x2)
x1hat <- lm(x1 ~ z + x2)$fitted.values
fit <- lm(y ~ x1hat + x2)

stargazer(fit, type="text")
summary(fit)




# Using AER package --------------------

ivreg(y ~ x1 + x2, ~ x2 + z)
iv.model <- ivreg(y ~ x1 + x2 | x2 + z)

summary(iv.model)
summary(iv.model, vcov = sandwich, diagnostics = TRUE)


# Minimal Working Example of Monte Carlo experiment in R -----------

M <- 10000
stats <- matrix(NA, nrow=M, ncol=2)
for (i in 1:M){
  n <- 200
  x <- rnorm(n)
  stats[i,1] <- mean(x)
  stats[i,2] <- sd(x)
}

mean(stats[,1])
mean(stats[,2])
apply(stats, 2, mean)
apply(stats, 2, sd)



# Monte Carlo experiment compare OLS and 2SLS  --------------
rm(list = ls()) 

M <- 1000
stats <- matrix(NA, nrow=M, ncol=2)

n <- 1000

for (i in 1:M){  
  
  x1. <- rnorm(n,0,2)
  e <- 2*x1. + rnorm(n, 0, 3) 
  z <- rnorm(n,0,5)
  x1 <- x1. + z
  x2 <- rchisq(n,2)
  y <- 5 + 2*x1 + 3*x2 + e 
  
  ols <- lm(y ~ x1 + x2)
  stats[i,1] <- ols$coefficients[2]
  x1hat <- lm(x1 ~ z + x2)$fitted.values
  two.sls <- lm(y ~ x1hat + x2)
  stats[i,2] <- two.sls$coefficients[2]
}


apply(stats, 2, mean)
apply(stats, 2, sd)

par(mfrow=c(1,2))
hist(stats[,1],breaks=20, main="OLS")
hist(stats[,2],breaks=20, main="2SLS")



# 10. Simulate z and e in a way that they have some dependence

# 11. Create x1 in a way that the dependece in z is negligible. 
#     The greater their correlation the stronger the instrument. 
#     The less their correlation the weaker the instrument. 
#     Identify the nature of problems with this.


# 12. Although there is no assumption about the instrument and 
#     the exogenous variables in the model, having a very high 
#     correlation will cause problems like collinearity. 
#     Note this is not strict, it's not needed to be equal 0,
#     just not that strong.