## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.6-0 (2020-12-14) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code
##
## Attaching package: 'rdist'
## The following object is masked from 'package:fields':
##
## rdist
## ERROR. SPATIAL DEPENDENCE IN ERROR.
library(spdep)
n <- 100 #number of observation
r <- 1000 #number of models for each level of rho
rho <- c(0, 0.2, 0.7, 0.9) #Vector of Rhos to test
X <- rnorm(n, 0, 3) #simulate Xs
w <- cell2nb(10, 10) #nb object for simulating autocorrelated data
bres <- data.frame() #holds the results of the simulatation
ctr <- 0 #counter
for (j in 1:length(rho)) {
iw <- invIrM(w, rho[j])
for (i in 1:r) {
ctr <- ctr + 1
e <- iw %*% rnorm(n)
y <- X + e
lm1 <- lm(y ~ X)
bres[ctr, 1] <- lm1$coef[2]
bres[ctr, 2] <- paste("Spatial AC Errors (rho=", rho[j], ")", sep = "")
}
}
## PLOTTING
library(ggplot2)
bres$V2 <- as.factor(bres$V2)
bres$V2 <- relevel(bres$V2, 4)
p1 <- ggplot(data = bres, aes(x = V1)) + geom_histogram() + facet_grid(. ~ V2)
p1 + geom_vline(xintercept=1, col='red')-When rho > 0 OLS estimates are both biased and inefficient.
### LAG. SPATIAL DEPENDENCE IN Y.
n <- 100
r <- 1000
rho <- c(0, 0.2, 0.7, 0.9)
X <- rnorm(n, 0, 3)
w <- cell2nb(10, 10)
bres <- data.frame()
ctr <- 0
for (j in 1:length(rho)) {
iw <- invIrM(w, rho[j])
for (i in 1:r) {
ctr <- ctr + 1
e <- rnorm(n)
y <- X + e
y <- iw %*% y
lm1 <- lm(y ~ X)
bres[ctr, 1] <- lm1$coef[2]
bres[ctr, 2] <- paste("Spatial AC Y (rho =", rho[j], ")", sep = "")
}
}
## PLOTTING
library(ggplot2)
bres$V2 <- as.factor(bres$V2)
bres$V2 <- relevel(bres$V2, 4)
p1 = ggplot(data = bres, aes(x = V1)) + geom_histogram() + facet_grid(. ~ V2)
p1 + geom_vline(xintercept=1, col='red')