Q:

Why is the intercept of linear regression biased?

Out of curiosity, I conducted the following simulation (code below). Why is it that when the variance of the error term is large coefficient associated with the intercept is biased? Can you recommend some reference that discusses this? Or better yet is there a formal proof of that?
rm(list=ls())
set.seed(12345)
m  <- 10000
x1 <- runif(m,0,100)       # random numbers from uniform distribution
x2 <- 1:10000
u  <- rnorm(m,0,100)       # random numbers from standard normal distribution
y  <- 5*x1 + 2*x2 + 10 + u # generating y series
data  <- cbind(x1, x2, y)
beta1 <- c()
beta2 <- c()
beta3 <- c()
R2    <- c()
n       <- 1000 # number of loops
ksubset <- 100  # length of subset
for (i in 1:n){
  datam <- data.frame(data[sample(nrow(data),ksubset), ])
  ols   <- summary(lm(y~x1+x2, data=datam))
  beta1 <- append(beta1, ols$coefficients[1])
  beta2 <- append(beta2, ols$coefficients[2])
  beta3 <- append(beta3, ols$coefficients[3])
  R2    <- append(R2,    ols$r.squared)
}

results <- c(mean(beta1), mean(beta2), mean(beta3))
results
## [1] 7.290909 5.027431 2.000345

The code takes thousand random samples of size 100 from the population and calculates a regression in each. Then I take the average of each of the estimated coefficients, which in theory should be equal to the original model. It works great for the slopes, but not for the intercept.
Update
I have noticed that with the population set at 10,000, the bias persists. If I increase the size of the population to 1,000,000, the bias disappears. The size of the target population was insufficient. Anyway, @Maarten's answer is a step in the right direction.
Update2
Entire population results:
> summary(lm(y~x1+x2))       
Call:
lm(formula = y ~ x1 + x2)

Residuals:
    Min      1Q  Median      3Q     Max 
-367.24  -66.37   -0.33   66.64  385.26 

Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept) 7.0388165  2.6626801    2.644  0.00822 ** 
x1          5.0296554  0.0348480  144.331  < 2e-16 ***
x2          2.0002959  0.0003465 5772.764  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 100 on 9997 degrees of freedom
Multiple R-squared:  0.9997,  Adjusted R-squared:  0.9997 
F-statistic: 1.667e+07 on 2 and 9997 DF,  p-value: < 2.2e-16

I updated the question, including simulation results (result). I also present the model results with the entire population, which in effect is biased due to the huge variance. The simulation, as expected, reproduce those results correctly. I repeat, my mistake is in the size of the target population.

A:

The role of the constant is to force the mean of the errors to be 0. In your program the mean of the error term is 0 in your population, but the variance is very high, so it will deviate a lot from 0 in many of your samples. This means that in each sample the constant can deviate a lot from the population value. So much so that a 1000 replications is not enough to find that on average you will get the right constant. Try and increase the number of replications to say 10,000 or 100,000 and see what happens.

Update
As promised I tried to replicate Héctor's results, but I gues I did something wrong the first time I tried, as I can't reproduce it anymore. Maybe I set the residual standard deviation at 1000 instead of 100, maybe I tried 100 replications instead of a 1000, or some silly typo like that. Anyhow here is my simulation with a 1000 replications in Stata, and a constant that is close enough to 10 given the spread. 
. set seed 123456

. clear all

. program define sim
  1.         drop _all
  2.         set obs 100
  3.         gen x1 = runiform() < .5
  4.         gen x2 = rnormal()
  5.         gen y  = 5*x1 + 2*x2 + 10 + rnormal(0,100)
  6.         reg y x1 x2
  7. end

. simulate b=_b[_cons], reps(1000) nodots : sim

      command:  sim
            b:  _b[_cons]

. sum

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
           b |      1,000     9.94305    14.99877  -46.69164   60.36443

