#ECNS 562 - Eric J. Belasco

#Simulation Exercise

#In-class demo



#Some initial assumptions

n = 5000

m <- 100

b0 <- 2; b1 <- 0.5

#Empty matrices to be populated

b <- matrix(0,m,2)

biv <- matrix(0,m,2)

pvals <- matrix(0,m,2)

r2 <- matrix(0,m,1)



#Start loop

i <- 1

while( i <= m){

	e <- rnorm(n,0,2)

	x <- runif(n,0,10)

	y <- b0 + b1*x + e

	ols <- summary(lm(y~x))





	b[i,] <- ols$coef[,1]

	pvals[i,] <- ols$coef[,4]

	r2[i] <- ols$r.squared



	i <- i + 1

}





#Plot results

par(mfrow=c(3,2))

hist(b[,1], main = "b0 hat");

hist(pvals[,1], main = "b0 pval");

hist(b[,2], main = "b1 hat");

hist(pvals[,2], main = "b1 pval");

hist(r2, main = "r-squared");





###########################################################################

#Instrumental Variables Excercies

############################################################################

i <- 1

while( i <= m){

	e <- rnorm(n,0,9)

	v <- runif(n,0,10)

	x <- .2*e + v

	y <- b0 + b1*x + e

	ols <- summary(lm(y~x))



	z <- 1 + .2*v + rnorm(1)



	b[i,] <- ols$coef[,1]



	xhat <- lm(x~z)$fitted.values

	iv <- summary(lm(y~xhat))

	biv[i,] <- iv$coef[,1]



	i <- i + 1

}





#Plot results

par(mfrow=c(2,1))

hist(b[,2], main = "b1 hat - OLS");

hist(biv[,2], main = "b1 hat - IV");