library(Matrix)set.seed(2)# number of observationsn <-100# true parametersbeta <- c(0,1,-1)rho <-0.75# design matrix with two standard normal variates as "covariates"X <- cbind(intercept=1, x=rnorm(n), y=rnorm(n))# sparse identity matrixI_n <- sparseMatrix(i=1:n, j=1:n, x=1)# number of nearest neighbors in spatial weight matrix Wm <-6# spatial weight matrix with m=6 nearest neighbors# W must not have non-zeros in the main diagonal!lat <- rnorm(n)long <- rnorm(n)W <- kNearestNeighbors(lat, long, k=6)# innovationseps <- rnorm(n=n, mean=0, sd=1)# generate data from model S <- I_n - rho * W
z <- solve(qr(S), X %*% beta + eps)y <- as.vector(z >=0)# 0 or 1, FALSE or TRUE# estimate SAR probit modelfit1 <- sar_probit_mcmc(y, X, W, ndraw=100, thinning=1, prior=NULL)plot(fit1, which=c(1,3), trueparam = c(beta, rho))