if (FALSE) {
require(spatialprobit)
# number of observations
n <- 10
# true parameters
beta <- 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 matrix
I_n <- sparseMatrix(i=1:n, j=1:n, x=1)
# number of nearest neighbors in spatial weight matrix W
m <- 6
# spatial weight matrix with m=6 nearest neighbors
W <- sparseMatrix(i=rep(1:n, each=m),
j=replicate(n, sample(x=1:n, size=m, replace=FALSE)), x=1/m, dims=c(n, n))
# innovations
eps <- 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 model
set.seed(12345)
sarprobit.fit1 <- sar_probit_mcmc(y, X, W, ndraw=500, burn.in=100,
thinning=1, prior=NULL, computeMarginalEffects=TRUE)
summary(sarprobit.fit1)
# print impacts
impacts(sarprobit.fit1)
################################################################################
#
# Example from LeSage/Pace (2009), section 10.3.1, p. 302-304
#
################################################################################
# Value of "a" is not stated in book!
# Assuming a=-1 which gives approx. 50
library(spatialprobit)
a <- -1 # control degree of censored observation
n <- 1000
rho <- 0.7
beta <- c(0, 2)
sige <- 0.5
I_n <- sparseMatrix(i=1:n, j=1:n, x=1)
x <- runif(n, a, 1)
X <- cbind(1, x)
eps <- rnorm(n, sd=sqrt(sige))
param <- c(beta, sige, rho)
# random locational coordinates and 6 nearest neighbors
lat <- rnorm(n)
long <- rnorm(n)
W <- kNearestNeighbors(lat, long, k=6)
y <- as.double(solve(I_n - rho * W) %*% (X %*% beta + eps))
table(y > 0)
# set negative values to zero to reflect sample truncation
ind <- which(y <=0)
y[ind] <- 0
# Fit SAR Tobit (with approx. 50% censored observations)
fit_sartobit <- sartobit(y ~ x,W,ndraw=1000,burn.in=200, showProgress=TRUE)
# print impacts
impacts(fit_sartobit)
}
Run the code above in your browser using DataLab