#  monte carlo experiment: multicollinearity

# attach library that contain multivariate normal proc

library(MASS)

# set sample size

n <- 200

b <- as.matrix(c(1,2,3))

# define 2x2 variance-covariance matrix for X variables

SigmaX <- matrix(c(1,0,0,1),2,2,byrow=T)

# define number of replications

sim <- 1000

# define vectors to hold OLS estimates

 b1 <- rep(0,sim)
 b2 <- rep(0,sim)
 corrvec <- rep(0,sim)

# start monte carlo loop

for(i in 1:sim){

 # create a matrix from multivariate normal draws, adding a column
 # of ones

 X <- cbind( rep(1,n),mvrnorm(n=200, rep(0, 2), SigmaX))

 e <- rnorm(n)

 # the DGP

 y <- X %*% b + e

 # do OLS

 bols <- (solve(t(X) %*% X)) %*% (t(X) %*% y)

 # save the estimates

 b1[i]<-bols[2]
 b2[i]<-bols[3]
 corrvec[i] <- cor(X[,2],X[,3])
}

# check for bias by taking means across replications
 
cat("mean of b1:", mean(b1), "\n")
cat("mean of b2:", mean(b2), "\n")

# demonstrate bias by substracting true value from means across
# replications 

cat("Bias is b1:", mean(b1)-b[2], "\n")
cat("Bias is b2:", mean(b2)-b[3], "\n")

cat("Mean correlation b/t X1 and X2:", mean(corrvec), "\n")