Spatial econometrics -- computing residuals

Question

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!

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