I am running some simulations, trying to use lagsarlm (from the spdep package) to recover the parameters used to generate the data. In a basic simulation I am running, I am finding that I am able to recover rho almost perfectly, and all but one of the betas perfectly. However, the beta attached to the constant is substantially biased, for reasons I cannot understand.

The code I am using is below. The spatial weights matrix is from the 48 contiguous U.S. states (rows sum to 1). Can anyone see where I am going wrong? Or is the biased B0 coefficient somehow a consequence of the particular neighbourhood structure I'm using? Any help or tips would be much appreciated.

- Malcolm

n <- dim(W)[1] # sample size

Bs <- c(2, 5, 3, -2) # vector of Beta coefficients

rho <- 0.2 # set autocorrelation coefficient

bres <- matrix(NA, nrow00, ncol=5)

for (i in 1:1000) {

+ e <- rnorm(n, mean=0, sd=4)Bs <- c(2, 5, 3, -2) # vector of Beta coefficients

rho <- 0.2 # set autocorrelation coefficient

bres <- matrix(NA, nrow00, ncol=5)

for (i in 1:1000) {

+ X1 <- rnorm(n, 4, 2) # create some independent variables

+ X2 <- rnorm(n, 2, 2)

+ X3 <- rnorm(n, -3, 1)

+ X <- cbind(rep(1, n), X1, X2, X3)

+ y <- (solve(diag(n)-rho*W)) %*% ((X%*%Bs)+e) # generate lagged Ys

+ data <- as.data.frame(cbind(y, X))

+ lagmod <- lagsarlm(y ~ X1 + X2 + X3, dataÚta, oxw.listw1)

+ bres[i,] <- coefficients(lagmod)

+ }

apply(bres, 2, mean)

[1] 0.1876292 2.6275783 4.9897827 3.0060831 -1.9883803

apply(bres, 2, median)

[1] 0.1907057 2.4000895 4.9887496 3.0043258 -2.0076960[1] 0.1876292 2.6275783 4.9897827 3.0060831 -1.9883803

apply(bres, 2, median)