rm(list=ls())

pdf(file="simulating-IVscenarios.pdf",width=12,height=9)
par(mfrow=c(2,3))

#simulation
set.seed(121053)
n     = 200
beta  = 2
N     = 1000
b.ols = rep(0,N)
b.2sls  = rep(0,N)

rho=0.1;delta=0.1
rho=0.1;delta=0.5
rho=0.5;delta=0.1
rho=0.5;delta=0.5

for (rho in c(0.1,0.5))
for (delta in c(0.1,0.5,1.0)){
tau   = sqrt(1-rho^2)
for (i in 1:N){
  w = rnorm(n)
  e = rnorm(n,rho*w,tau)
  z = rnorm(n)
  x = delta*z + w
  y = beta*x  + e
  b.ols[i] = lm(y~x-1)$coef
  xhat = lm(x~z-1)$fit
  b.2sls[i] = lm(y~xhat-1)$coef
}

plot(density(b.ols),xlim=c(0,4),main="",xlab="beta hat",lwd=2)
lines(density(b.2sls),lty=2,lwd=2)
points(beta,0,col=2,pch=16,cex=2)
title(paste("rho=",rho," - delta=",delta,"\n OLS (solid), 2SLS (dashed)",sep=""))
}


dev.off()