# Biased estimates when intercept is included in a linear regression

I am simulating 10000 data-sets, each of length 20, that follow an autoregressive model with lag 1, using the following code:

set.seed(1)
N = 20
n.reps <- 10000
burn.in = 50
total <- N + burn.in
x <- matrix(NA, n.reps, total)
x[ , 1] <- 0.1
noise <- rnorm(n.reps*total, 0, 0.1)

for(j in 1:n.reps)
for(i in 2:total)
x[j, i] <- 0.5*x[j, i-1] + noise[j + i]

x <- x[ , -(1:burn.in) ]


Then I estimated the autoregressive coefficients (0.5 in this case) for each observed path, using two linear regressions:

res1 <- res2 <- numeric(n.reps)

for(i in 1:n.reps)
{
res1[i] <- lm(x[i, 2:N] ~ -1 + I(x[i, 1:(N-1)]))$coef res2[i] <- lm(x[i, 2:N] ~ I(x[i, 1:(N-1)]))$coef[2]
}

mean(res1) # 0.4687619
mean(res2) # 0.3845817


In the first line I'm fitting the right linear model, while in the second I'm including an intercept that is not present in the data generating process. I expected both methods to give me unbiased estimates of the autoregressive coefficient, but it looks like including the intercept makes the estimates be biased downward. The bias disappears as I increase the sample size N.

The second linear model nests the real one, so I expected higher variance in the estimates when the larger model is fitted, but not bias. So my question is: am I doing something wrong? Thanks!

• Note that you don't need to burn in the simulations, as you're not doing MCMC. Also how much bias are you seeing? And how many simulations are you doing? Your simulation based estimate of bias could have a large simulation variance if the number of simulations is not large enough Commented Dec 29, 2012 at 14:28
• I'm including a burn in because I don't want my answer to depend on the initialization of the AR(1) process. I've used fairly big sets of simulations (> 10000), and the bias can be around 0.1 if the sample size is small (N = 20). The bias progressively disappears as N increases. Commented Dec 29, 2012 at 14:39
• @prob: The point of the "burn in" is that this process does not start at the stationary distribution and so, just as in MCMC, using some burn in period will get you somewhat closer. In contrast to the general MCMC situation, it's easy to characterize the rate of convergence and also almost as easy to just start right at the stationary distribution (which the OP hasn't done). This also relates to the concept of transient response in signal processing, which the OP might want to investigate. Commented Dec 29, 2012 at 14:41
• How are you assessing the bias in this simulation? I've replicated your code and not found any indication that they give significantly different results. (10000 reps) Commented Dec 29, 2012 at 19:34
• @AdamO I have included the whole Monte Carlo experiment I have used. But more importantly I found a mistake in my assumptions (see the answer). Commented Dec 31, 2012 at 15:12

There is a mistake in my question: OLS doesn't give an unbiased estimator for the autoregressive coefficients of an $\mathrm{AR}(p)$ process! In fact the OLS estimator is biased toward zero as explained in this article:
A. Maeshiro (2000), An Illustration of the Bias of OLS for $Y_t = \lambda Y_{t-1} + U_t$, J. Econ. Education, vol. 31, no. 1, 76–80.
Then it appears that including an intercept makes things worse, since the estimates are even more biased. Also: the OLS is consistent, that's why the bias is disappearing as $N$ increases.