Simulate spatially autocorrelated data

Given a spatial weights matrix and degree of autocorrelation, returns autocorrelated data.

sim_sar(
  m = 1,
  mu = rep(0, nrow(w)),
  rho,
  sigma = 1,
  w,
  type = c("SEM", "SLM"),
  approx = FALSE,
  K = 20,
  ...
)

Source

LeSage, J. and Pace, R. K. (2009). An Introduction to Spatial Econometrics. CRC Press.

Arguments

m: The number of samples required. Defaults to m=1 to return an n-length vector; if m>1, an m x n matrix is returned (i.e. each row will contain a sample of auto-correlated values).
mu: An n-length vector of mean values. Defaults to a vector of zeros with length equal to nrow(w).
rho: Spatial autocorrelation parameter in the range (-1, 1). A single numeric value.
sigma: Scale parameter (standard deviation). Defaults to sigma = 1. A single numeric value.
w: n x n spatial weights matrix; typically row-standardized.
type: Type of SAR model: spatial error model ("SEM") or spatial lag model ("SLM").
approx: Use power of matrix W to approximate the inverse term?
K: Number of matrix powers to use if approx = TRUE.
...: further arguments passed to MASS::mvrnorm.

Value

If m = 1 then sim_sar returns a vector of the same length as mu, otherwise an m x length(mu) matrix with one sample in each row.

Details

This function takes n = nrow(w) draws from the normal distribution using rnorm to obtain vector x; if type = 'SEM', it then pre-multiplies xby the inverse of the matrix(I - rho * W)to obtain spatially autocorrelated values. Fortype = 'SLM', the multiplier matrix is applied to x + mu to produce the desired values.

The approx method approximates the matrix inversion using the method described by LeSage and Pace (2009, p. 40). For high values of rho, larger values of K are required for the approximation to suffice; you want rho^K to be near zero.

Examples

# spatially autocorrelated data on a regular grid
library(sf)
row = 10
col = 10
sar_parts <- prep_sar_data2(row = row, col = col)
w <- sar_parts$W
x <- sim_sar(rho = 0.65, w = w)
dat <- data.frame(x = x)

# create grid 
sfc = st_sfc(st_polygon(list(rbind(c(0,0), c(col,0), c(col,row), c(0,0)))))
grid <- st_make_grid(sfc, cellsize = 1, square = TRUE)
st_geometry(dat) <- grid
plot(dat)

# draw form SAR (SEM) model
z <- sim_sar(rho = 0.9, w = w)
moran_plot(z, w)
grid$z <- z

# multiple sets of observations
# each row is one N-length draw from the SAR model
x <- sim_sar(rho = 0.7, w = w, m = 4)
nrow(w)
dim(x)
apply(x, 1, aple, w = w)
apply(x, 1, mc, w = w)

# Spatial lag model (SLM): y = rho*Wy + beta*x + epsilon
x <- sim_sar(rho = 0.5, w = w)
y <- sim_sar(mu = x, rho = 0.7, w = w, type = "SLM")

# Spatial Durbin lag model (SLM with spatial lag of x)
# SDLM: y = rho*Wy + beta*x + gamma*Wx + epsilon
x = sim_sar(w = w, rho = 0.5)
mu <- -0.5*x + 0.5*(w %*% x)[,1]
y <- sim_sar(mu = mu, w = w, rho = 0.6, type = "SLM")