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!