Simulating a bimodal biased IV estimator How can I simulate a bimodal biased IV estimator?
The common unimodal heavy-tailed biased estimator would be something like this:
IV_data <- function(beta=10, lambda=-20, theta=0.5, n=50, mu=100, sd=2){
  z = rnorm(n,mu, sd)
  u = rnorm(n, 0, sd)
  x = theta*z + u 
  y = beta*x + lambda*u + rnorm(n, mu, sd)
  data.frame(y, x, z, u)
}

library(sem)
set.seed(20)
beta_IV <- replicate(1000, coef(tsls(y~x, ~z, data=IV_data()))[2])
d<-density(beta_IV)
plot(d)

I guess if I tweak the parameters a bit I could get a bimodal version. Or do I have to have a complete different model structure to get it?
 A: You can get this behaviour of the 2SLS estimator if you have a high degree of endogeneity and weak instruments.
Specification of the structural system
Here is a complete specification of the simple bivariate system that I will be simulating.
$$
\begin{align}
Y_i &= X_i\beta + U_i \\
X_i &= Z_i \pi + V_i
\end{align}
$$
where the errors are jointly normally distributed, 
$$
\begin{bmatrix}U_i \\ V_i\end{bmatrix} \sim \text{N}\left(\begin{bmatrix}0 \\ 0\end{bmatrix}, \begin{bmatrix}\sigma^2_U & \rho \sigma_V\sigma_U \\ \rho \sigma_V\sigma_U & \sigma^2_V\end{bmatrix}\right)
$$
and the first-stage regression coefficient is related to the concentration parameter $\mu$, 
$$
\mu = \frac{\pi^2 \mathbf{Z}'\mathbf{Z}}{\mathbb{V}(V_i)}
$$
and $\mathbf{Z} = [Z_1, \ldots, Z_N]'$. 
Typically, the strength of instruments is specified via the concetration parameter (see, for example, Bowden and Turkington, 1984), and that is the main control parameter I will use in my simulations.
For the particular case that you want to see, I set the concentration parameter to a very low value (0.25), and the degree of endogeneity, given by the covariance of the error to a very high number (0.99).
Note that this is a simplified version for the simple bivariate case. General formulae can be found in Flores-Lagunes (2007), for example.
R implementation
Here is some code that simulates the estimates from the system above for the case that the true value of the structural parameter $\beta$ is zero. The sample size
is set to 1000, and I replicate the system 1000 times.
library(MASS)
library(sem)
library(ggplot2)

# function to return the structural parameter 
fn2SLS = function(iSampleSize, dConcParam, mCovar = NULL) {

  # generate the instrument
  vZ = rnorm(iSampleSize)
  dZZ = vZ %*% vZ

  # the value of the first stage parameter implied
  #   by the concentration parameter
  dPi = sqrt(dConcParam/dZZ)

  # covariance matrix of the errors
  if(is.null(mCovar)){
    mCovar = 
      matrix(c(1, .99, .99, 1), nrow = 2, byrow = TRUE)
  } 

  # generate the structural errors 
  mErr = mvrnorm(n = iSampleSize, mu = c(0, 0), Sigma = mCovar)

  # generate the endogenous regressor
  vX = dPi * vZ + mErr[, 1]

  # generate the outcome variable under the null of zero beta
  vY = mErr[, 2]

  # compute the tsls estimator
  dBeta = coef(tsls(vY ~ vX, ~ vZ, data = data.frame(vY, vX, vZ)))[2]
  return(dBeta)
}

# test the function
# debugonce(fn2SLS)
fn2SLS(iSampleSize = 1000, dConcParam = 0.25)

# simulate the 2SLS estimator for the given value of the concentration
#   parameter (sample size = 1000)
vBetaSim = sapply(1:1000, 
              function(x) fn2SLS(iSampleSize = 1000, dConcParam = 0.25))

# plot the kernel density of the estimated coefficients
ggplot(data = data.frame(x = vBetaSim[abs(vBetaSim) <= 5]), aes(x = x)) + geom_area(stat = 'density', fill = 'lightblue') +
  geom_line(stat = 'density', color = 'blue') + theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  geom_vline(xintercept = c(0, vBetaSim[which.max(density(vBetaSim)$y)]),
             color = 'purple') +
  ggtitle('Biased bimodal 2SLS estimator') +
  xlab('2SLS estimates of structural parameter')

If you look at the kernel density plot of the estimates, you will see that the estimates are biased, centred around the purple line (true value being
zero denoted by the red line). You will also note that there are atleast two other nodes to the right of the highest mode.

Note: I am aware that the purple line does not quite go through the mode, but it is supposed to.
