How to simulate a SAR model? I'm trying to simulate a SAR data generating process but I'm not getting the excepted result. To be clear, by SAR model I mean:
$y = \rho Wy + X\beta + e$
which can be rewritten as
$y = (I - \rho W)^{-1} (X\beta + e)$
I want to show that using OLS will result in biased coefficients (see this answer). But when I run this code I get unbiased estimates.
x1s_ols <- NULL
x2s_ols <- NULL

for (i in 1:1000) {
  N <- 100
  
  x1 <- rnorm(N)
  x2 <- rnorm(N, mean = 1 + 2*x1)
  e <- rnorm(N, 0, 2)
  
  W <- matrix(abs(rnorm((N)^2)), nrow=N, ncol=N)
  diag(W) <- 0 #remove diagonal
  W[lower.tri(W)] <- t(W)[lower.tri(W)] #force symmetric
  #row standardise
  W <- W / rowSums(W)
  rho <- 0.75 #runif(1)
  
  #calculate spatial process - solve = inverse
  y <- solve(diag(N) - rho * W) %*%  (2 + 3*x1 + 4*x2 + e)
  
  df <- data.frame(x1, x2, y) 
  
  m <- lm(y ~ x1 + x2, data=df)
  x1s_ols[i] <- coef(summary(m))[2]
  x2s_ols[i] <- coef(summary(m))[3]
}

#both unbiased
print(mean(x1s_ols))
print(mean(x2s_ols))

Could someone explain what I'm doing wrong, or what I need to change in my code to get obviously biased coefficients?

I think maybe what is happening is with your random weights matrix, it is very diffuse, so the bias overall is not very large? (Although I am getting betas a tinge lower than expected, which was contra to what I thought would happen.)
If I do a fixed grid with Queen's contiguity (row-normalized) weight matrix I get more expected results with beta's too large:
set.seed(10)
sim <- 1000
x1s_ols <- rep(-999,1000)
x2s_ols <- rep(-999,1000)

# Lets do a fixed W on a grid
N <- 100
x <- 1:10
y <- 1:10
grid <- expand.grid(x,y)
W <- as.matrix(dist(grid,diag=TRUE,upper=TRUE))
W <- 1*(W < 2) # queen
diag(W) <- 0
W <- W/rowSums(W)
rho <- 0.75
left <- solve(diag(N) - rho * W)

for (i in 1:sim) {

  x1 <- rnorm(N)
  x2 <- rnorm(N, mean = 1 + 2*x1)
  e <- rnorm(N, 0, 2)
  
  #calculate spatial process - solve = inverse
  y <- left %*%  (2 + 3*x1 + 4*x2 + e)
  
  df <- data.frame(x1, x2, y)
  
  m <- lm(y ~ x1 + x2, data=df)
  x1s_ols[i] <- coef(summary(m))[2]
  x2s_ols[i] <- coef(summary(m))[3]
}

#These are now too large
print(mean(x1s_ols))
print(mean(x2s_ols))

Which I am getting:
# > print(mean(x1s_ols))
# [1] 3.424781
# > print(mean(x2s_ols))
# [1] 4.561689

So maybe if you really want to simulate a random weights matrix, simulate some real spatial data on the X/Y plane and redo the W matrix? Continuous distances should be more like rgamma or rnorm(N^2)^2 (right skewed), than abs(rnorm(N^2)) I think.
