Consider the model in spatial econometrics denoted SAR by James P LeSage:

$y = \rho W y + X \beta+ \epsilon$

I use the R package spdep and the econometric toolbox by LeSage from www.spatial-econometrics.com. However, when comparing the way these two compute the residuals there are differences. The R package by Roger Bivand computes them by

$\dot{y}=y - \rho W y$

and then linearly regressing $\dot{y}$ on the $X$ matrix. Obtaining residuals from this fit. (lm(y-rWy~X-1). Yet LeSage computes the residuals the way (I originally thought would be the right way):

$y-(\rho W y)^{-1}Xb$

So here're my questions:

  1. Why do the two authors compute the residuals differently?
  2. Is there an objection against one or the other way of computing them?

Thank you for kind help!


1 Answer 1


The residuals that each of them calculating are different. Here is why:

The model is as follows:

$y = \rho Wy + xb + e$ with $e \sim n(0,1)$

Now if we play arround with it we get:

$y = (I - \rho W)^{-1}(xb + e)$

Now what Prof. LeSage does is:

$y - (I-\hat{\rho} \cdot W)^{-1} \cdot x\hat{b} = (I-\rho W)^{-1}\cdot e$

So what you are getting it the residual with the auto correlation.

On the other hand, by transforming y:

$y - \hat{\rho}\cdot W \cdot y = xb + e$

Estimating, $xb$ and calculation the residuals, what Bivand is doing is giving you $e$ instead of $(I-\rho W)^{-1}\cdot 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
g <- erdos.renyi.game(n=n, p.or.m=p, type="gnp", directed=F, loops=F)

w <- get.adjacency(g) #get an adjacency matrix
w <- w/rowSums(w)     #row standardize because of eigen vectors and eigen values

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
cor(z0, w%*%y)
cor(z1, w%*%y)

The instruments work because rvariable is exogenous. So as long as w is exogenous we have a game!

#------------------ NOW ONTO ESTIMATION

#The wrong way ...
summary(out<-lm(y ~ rvariable)) 

#The not so bad, but still very wrong way
summary(out<-lm(y ~ w%*%y + rvariable)) 

#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 )) 

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

  • $\begingroup$ i understand you reasoning but on the other hand the matrix $(I-\rho W)$ needs to bee inverted for estimation anyway - right? at least it gets inverted within the lagsarlm() function implemented. $\endgroup$
    – Seb
    Jan 14, 2012 at 18:53
  • $\begingroup$ No that is not right (I-pW) does not need to be inverted. It might be the case that is is in that code, I don't know (never looked), but this model can be estimated by "2SLS" as in Prucha with no need to inver (I-pW) since you instrument Wy with Wx ... by maximum likelihood in which case you maximize the concentrated likelihood function lnL(ρ) = κ + ln|In − ρW| − (n/2) ln(S(ρ)) over "p" or via MCMC as in LeSage where they avoid direct inversion solvin (In−ρW)μ = X via cholsky decomposition ... $\endgroup$
    – mmgm
    Jan 14, 2012 at 19:18
  • $\begingroup$ thanks for your extensive answer and help. if i may summarise the bottomline is: the bivand way is computationally simpler and therefore preferable. well, that's a good hint - since the results in my application differ quite substantially - which i find rather strange but is probably due to misspecification (?). anyway, i would have hoped for a statistical reasoning ;) $\endgroup$
    – Seb
    Jan 15, 2012 at 20:09
  • 1
    $\begingroup$ I am revising my answer since when i went to get you the statistical reasoning i realized i hasted into a conclusion that was not right $\endgroup$
    – mmgm
    Jan 16, 2012 at 10:41
  • 1
    $\begingroup$ @Seb i have posted some code to prove my new answer. Sorry if I mislead you, Miguel $\endgroup$
    – mmgm
    Jan 16, 2012 at 10:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.