Rare event logistic regression bias: how to simulate the underestimated p's with a minimal example? CrossValidated has several questions on when and how to apply the rare event bias correction by King and Zeng (2001).  I am looking for something different:  a minimal simulation-based demonstration that the bias exists.
In particular, King and Zeng state

"... in rare events data the biases in probabilities can be
  substantively meaningful with sample sizes in the thousands and are in
  a predictable direction: estimated event probabilities are too small."

Here is my attempt to simulate such a bias in R:
# FUNCTIONS
do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)  
    # Extract the fitted probability.
    #    If p is constant, glm does y ~ 1, the intercept-only model.
    #    If p is not constant, assume its smallest value is p[1]:
    glm(y ~ p, family = 'binomial')$fitted[1]
}
mean.of.K.estimates = function(p, K){
    mean(replicate(K, do.one.sim(p) ))
}

# MONTE CARLO
N = 100
p = rep(0.01, N)
reps = 100
# The following line may take about 30 seconds
sim = replicate(reps, mean.of.K.estimates(p, K=100))
# Z-score:
abs(p[1]-mean(sim))/(sd(sim)/sqrt(reps))
# Distribution of average probability estimates:
hist(sim)

When I run this, I tend to get very small z-scores, and the histogram of estimates is very close to centered over the truth p = 0.01.  
What am I missing?  Is it that my simulation is not large enough show the true (and evidently very small) bias?  Does the bias require some kind of covariate (more than the intercept) to be included?
Update 1:  King and Zeng include a rough approximation for the bias of $\beta_0$ in equation 12 of their paper.  Noting the N in the denominator, I drastically reduced N to be 5 and re-ran the simulation, but still no bias in the estimated event probabilities is evident.  (I used this only as inspiration. Note that my question above is about estimated event probabilities, not $\hat \beta_0$.)
Update 2: Following a suggestion in the comments, I included an independent variable in the regression, leading to equivalent results:
p.small = 0.01
p.large = 0.2
p = c(rep(p.small, round(N/2) ), rep(p.large, N- round(N/2) ) )
sim = replicate(reps, mean.of.K.estimates(p, K=100))

Explanation: I used p itself as the independent variable, where p is a vector with a repetitions of a small value (0.01) and a larger value (0.2).  In the end, sim stores only the estimated probabilities corresponding to $p = 0.01$ and there is no sign of bias.
Update 3 (May 5, 2016):
This does not noticeably change the results, but my new inner simulation function is
do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(0)
    }else{
        # Extract the fitted probability.
        #    If p is constant, glm does y ~ 1, the intercept only model.
        #    If p is not constant, assume its smallest value is p[1]:
        return(glm(y ~ p, family = 'binomial')$fitted[1])
    }
}

Explanation:  The MLE when y is identically zero does not exist (thanks to comments here for the reminder).  R fails to throw a warning because its "positive convergence tolerance" actually gets satisfied.  More liberally speaking, the MLE exists and is minus infinity, which corresponds to $p=0$; hence my function update.  The only other coherent thing I can think of doing is discarding those runs of the simulation where y is identically zero, but that would clearly lead to results even more counter to the initial claim that "estimated event probabilities are too small".
 A: Figure 7 in the paper seems to most-directly address the question of bias in the predictions. I don't fully understand the figure (specifically, the interpretation "estimated event probabilities are too small" seems like an oversimplification) but I managed to reproduce something similar to it based on their terse description of their simulation in Section 6.1:
n_grid = 40
x_grid = seq(0, 7, length.out = n_grid)
beta0 = -6
beta1 = 1

inverse_logit = function(x) 1/(1 + exp(-x))

do.one.sim = function(){
    N = 5000
    x = rnorm(N)
    p = inverse_logit(beta0 + beta1*x)
    # Draw fake data based on probabilities p
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(rep(0, n_grid))
    }else{
        # Extract the error
        mod = glm(y ~ x, family = 'binomial')
        truth = inverse_logit(beta0 + beta1*x_grid)
        pred = predict(mod, newdata = data.frame(x = x_grid),
            type = 'response')
        return(pred - truth)
    }
}
mean.of.K.estimates = function(K){
    rowMeans(replicate(K, do.one.sim()))
}

set.seed(1)
bias = replicate(10, mean.of.K.estimates(100))
maxes = as.numeric(apply(bias, 1, max))
mins = as.numeric(apply(bias, 1, min))

par(mfrow = c(3, 1), mar = c(4,4,2,2))
plot(x_grid, rowMeans(bias), type = 'l',
    ylim = c(min(bias), max(bias)),
    xlab = 'x', ylab = 'bias')
lines(x_grid, maxes, lty = 2)
lines(x_grid, mins, lty = 2)
plot(x_grid, dnorm(x_grid), type = 'l',
    xlab = 'x', ylab = 'standard normal density')
plot(x_grid, inverse_logit(beta0 + beta1*x_grid),
    xlab = 'x', ylab = 'true simulation P(Y = 1)',
    type = 'l')

The first plot is my replication of their figure 7, with the addition of dashed curves representing the full range of results over 10 trials. 
As per the paper, x here is a predictor variable in the regression drawn from a standard normal. Thus, as illustrated in the second plot, the the relative frequency of observations for x > 3 (where the most visible bias occurs in the first plot) is diminishingly small.
The third plot shows the "true" simulation probabilities in the generating process as a function of x. It appears that the greatest bias occurs where x is rare or non-existent.
Taken together, these suggest that the OP entirely misinterpreted the central claim of the paper by confusing "rare event" (i.e. x > 3) with "event for which P(Y = 1) is very small". Presumably the paper concerns the former instead of the latter.

A: Rare events bias only occurs when there are regressors. It won't occur in an intercept-only model like the one simulated here. See this post for details:
 http://statisticalhorizons.com/linear-vs-logistic#comment-276108
