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_hat*W*y = 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/rowSum(w)      #row standardize because of eigen vectors and eigen values
    sum(rowSums(w)==0) 
    
    rho <- 0.9
    intercept   <- rep(1,n)
    rvariable   <- runif(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
    z2 <- w%*%w%*%w%*%rvariable

    #check to see if there is a minimum of correlation ... it shoudl
    cor(z0, w%*%y)
    cor(z1, w%*%y)
    cor(z2, 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)
    
    #The right way if instruments permit 
    #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 + z2 + 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)