I went to look into the code in detail and actually my answer was not correct. Sorry i hasted and i should not have.
The residuals that each of them calculating are different. Here is why:
The model is as follows:
y = pWy + xb + e with e ~ n(0,1)
Now if we play arround with it we get:
y = (I - pW)^-1(xb + e)
Now what Prof. LeSage does ie:
y - (I-p_hat * W)^-1 * xb_hat = (I-pW)^-1*e
So what you are getting it the residual with the auto correlation.
On the other hand, by transforming y:
y - p_hatWy = xb + e
Estimating, xb and calculation the residuals, what Bivand is doing is giving you e instead of (I-pW)^-1*e
Which one is preferred will depend on your application!
Here is some code to prove it. I am not using SPDEP directly because i am not sure how to create random maps... But that is ok the code is pretty simple anyway:
#------------------ GENERATE SAMPLE DATA
rm(list=ls()) #clean
require(igraph) #random graphs
require(AER) #get ivreg ...
n<-700 #700 locations
p=0.2
g <- erdos.renyi.game(n=n, p.or.m=p, type="gnp", directed=F, loops=F)
graph.density(g)
w <- get.adjacency(g) #get an adjacency matrix
w <- w/rowSums(w) #row standardize because of eigen vectors and eigen values
sum(rowSums(w)==0)
rho <- 0.5
intercept <- rep(1,n)
rvariable <- rnorm(n)
y <- solve(diag(n) - rho*w) %*% ( 2*intercept + 3*rvariable + rnorm(n))
After the data is generated according to a SAR LAG model we will estimated it via 2SLS (as i told you we could).
#------------------ GENERATE INSTRUMENTS
#get some instrumental variables
z0 <- w%*%rvariable
z1 <- w%*%w%*%rvariable
#check to see if there is a minimum of correlation ... it shoudl
cor(z0, w%*%y)
cor(z1, w%*%y)
The instruments work because rvariable is exogenous. So as long as w is exougenous we have a game!
#------------------ NOW ONTO ESTIMATION
#The wrong way ...
summary(out<-lm(y ~ rvariable))
confint(out)
#The not so bad, but still very wrong way
summary(out<-lm(y ~ w%*%y + rvariable))
confint(out)
#ok now this should do it ... not perfect beacuse 2sls is not efficient.
#I am doing it this way because i did not want to generate random maps...
#Plus random graphs are easily available !
summary(out<-ivreg( y ~ w%*%y + rvariable, instruments=~ z0 + z1 + rvariable ))
confint(out)
Now to what really matters, the computation of residuals:
#residuals LeSage way
y_hat0 <- solve(diag(n) - coef(out)[2]*w ) %*% ( coef(out)[1]*intercept + coef(out)[3]*rvariable )
u_hat0 <- y - y_hat0
#residuals BiVand way
y_tilda <- y - coef(out)[2]*w%*%y
summary(out_biv <- lm( y_tilda ~ rvariable ))
#ok they are not the same due to rounding error ...
coef(out)[3] == coef(out_biv)[2]; round(coef(out)[3],5) == round(coef(out_biv)[2],5)
u_hat1 <- residuals(out_biv)
u_hat1 <- solve(diag(n) - coef(out)[2]*w)%*%u_hat1
#If we give Bivand some taste of autocorrelation it is the same as LeSage ...
round( u_hat0 - u_hat1, 5)
In the end you should see the residuals difference == 0 !
A cautionary note here is that depending on the structure of w the effect might not be identifiable so the strategy of using the random graph generator might be bogus some times !
Anyway I hope this really solved your question