Huge bias in IV during Monte Carlo Simulations

Question

I am trying to see how IV performs with Monte Carlo Simulations. My model is: $y = X \beta + p \alpha + \xi + \epsilon $. In this model $ \xi $ is not observed and $p$ is correlated with $\xi$. So I use an instrumental variable $z$ which is correlated with $p$ but is not correlated with $\xi + \epsilon$. However, for $\alpha$, I get erratic results. For example, I set $\alpha = -5$, however, I get values as high as $\pm 400$ for $\alpha$ in some of the simulations. This occurs when I use around 100 observations which to me is surprising. I couldn't figure out what I am missing. Any help will be appreciated.

If you also have a reach a well-performing IV Monte Carlo Simulation, it will help me to spot my problem. My Code is:

#Cleaning
rm(list = ls())

library(MASS)
library(AER)

#Number of observations
obs <- 100
#Number of characteristics
k <- 4

#Assign the mean values of X
mu.x <- rep(0, k)

#Assign the variance covariance matrix of X
sigma.x <- matrix(diag(rep(2, k)), ncol = k)

#Sigma for correlated price and instrumental variable
sigma.ins <- matrix(c(1, 0.5, 0.5, 1), ncol = 2)


nsim = 500

#slope values for x characteristics
coeff.x <- as.vector(c(4, 4, 2, 2, rep(0, k-4)))
beta.p <- -5

#Assign names for the columns
xnam <- paste("X", 1:k, sep = "")

#OLS regression formula
fm.olsreg = paste0("y ~ ", 
                   paste(xnam, collapse = " + "), " + ", 
                   paste("p", collapse = " + "))

#IV regression formula
fm.ivreg = paste0("y ~ ", 
                  paste(xnam, collapse = " + "), " + ", 
                  paste("p", collapse = " + "),  " | ",
                  paste(xnam, collapse = " + "), " + ",
                  paste("z", collapse = " + "))


#Storing p coefficient for IV and OLS
p.values <- matrix(NaN, nrow = nsim, ncol = 2)

for (sim in 1:nsim) {
  #Creating x characteristics for each of j products in t markets
  X.jt <- mvrnorm(obs, mu = mu.x, Sigma = sigma.x)
  
  #The base values for price and instrumental variable
  p.raw.z <- mvrnorm(obs, mu = c(0,0), Sigma = sigma.ins)
  
  #Creating the endogenous part
  xi.jt <- as.matrix(rnorm(obs, 0, 3/2))
  
  #real p data corrlated with xi
  p.jt <- p.raw.z[,1] + xi.jt
  
  #instrumental variable correlated with p
  z.jt <- p.raw.z[,2] + rnorm(obs, 0, 0.3)
  
  #creating errors
  e.jt <- rnorm(obs, 0, 0.8)
  
  #Producing y variables
  y <- X.jt %*% coeff.x + p.jt*beta.p + xi.jt +  e.jt
  
  #putting everything together
  data <- as.matrix(cbind(X.jt, p.jt, xi.jt, z.jt, y))
 
  #Assigning names to related columns
  colnames(data)[1:as.integer(k+4)] <- c(xnam, 'p', 'xi', 'z', 'y')
  
  #OLS regression
  ols = lm(data = as.data.frame(data),
           formula = fm.olsreg)
  
  #IV regression
  mo.ivreg <- ivreg(data = as.data.frame(data),
                    formula = fm.ivreg)
  
  #Store p values for both method
  p.values[sim,1] <- mo.ivreg$coefficients['p']
  p.values[sim, 2] <- ols$coefficients['p']
  
}

plot(density(p.values[,1]))
plot(density(p.values[,2]))

summary(p.values)

You may get for $\alpha$ something like

$$\begin{array}{c|c|c|} & \text{IV} & \text{OLS} \\ \hline \text{Min} & -15.843 & -4.530 \\ \hline \text{Mean} & -4.575 & -4.309 \\ \hline \text{Max} & 285.922 & -4.083 \\ \hline \end{array}$$

Thanks for any help.