Los residuos que calcula cada uno de ellos son diferentes. He aquí por qué:
El modelo es el siguiente:
$y = \rho Wy + xb + e$ con $e \sim n(0,1)$
Ahora bien, si jugamos con ello obtenemos:
$y = (I - \rho W)^{-1}(xb + e)$
Ahora lo que el Prof. LeSage hace es:
$y - (I-\hat{\rho} \cdot W)^{-1} \cdot x\hat{b} = (I-\rho W)^{-1}\cdot e$
Así que lo que se obtiene es el residuo con la correlación automática.
Por otra parte, transformando y:
$y - \hat{\rho}\cdot W \cdot y = xb + e$
Estimación, $xb$ y calcular los residuos, lo que Bivand está haciendo es darte $e$ en lugar de $(I-\rho W)^{-1}\cdot e$
La elección dependerá de su aplicación.
He aquí un código que lo demuestra. No estoy usando SPDEP directamente porque no estoy seguro de cómo crear mapas aleatorios... Pero eso está bien el código es bastante simple de todos modos:
#------------------ 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))
Una vez generados los datos según un modelo SAR LAG, los estimaremos mediante 2SLS (como ya te dije que podíamos hacer).
#------------------ 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)
Los instrumentos funcionan porque rvariable es exógena. Por tanto, mientras w es exógena ¡tenemos un juego!
#------------------ 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)
Ahora pasemos a lo que realmente importa, el cálculo de los residuos:
#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)
¡Al final debería ver la diferencia de residuos == 0 !
Hay que tener en cuenta que, dependiendo de la estructura del $W$ el efecto podría no ser identificable, por lo que la estrategia de utilizar el generador de gráficos aleatorios podría ser falsa en algunas ocasiones.
De todos modos espero que esto realmente resuelto su pregunta