# MLE estimation of spatial effects radius

I am trying to identify the maximum likelihood estimates in an SDM model (a hedonic home price model, with observations being 5,000 individual homes), $$y=\rho W y + X\beta + WX\lambda + \epsilon$$

The likelihood function of this model is $$\ln(\mathcal{L}) = \dfrac{n\ln(\pi\sigma^2)}{2} + \ln|I - \rho W| + \dfrac{e'e}{2\sigma^2}$$ with $e= y - \rho W y - X\beta - WX\lambda$

I have been able to successfully estimate parameters $\rho, \beta,\lambda$ using the lagsarlm function in the spdep package. In my case, $W$ is an inverse distance matrix limited to a radius $r$, or $$W_{ij} = \begin{cases} 1/d_{ij} & \text{for d_{ij} \leq r miles}\\ 0 & \text{for d_{ij} > r miles} \end{cases}$$ I would like to identify $r$ simultaneously with $\rho, \beta,\lambda$. I have even written a likelihood function to do this in R (making use of the neighbors and lag functions in spdep)

sdm.lik <- function(params, y, X, GeoPoints){
n <- length(y)
rho <- params[1]
sigma2 <- params[3]
delta <- params[-c(1,2,3)]
# make W subject to radius
dists      <- nbdists(points.dnn, GeoPoints)
dists.inv  <- lapply(dists, function(x) 1/x)
Wlist <- nb2listw(points.dnn, glist=dists.inv, style="W", zero.policy=TRUE)
W <- listw2mat(Wlist)
# Lag X variables
Z <- cbind(X[,1], lag.listw(Wlist, X[,1], zero.policy=TRUE))
for(i in 2:ncol(X)){
Z <- cbind(Z, X[,i], lag.listw(Wlist, X[,1], zero.policy=TRUE))
}
# Calculate log likelihood
e <- y - ((rho * W) %*% y) - (Z %*% delta)
lagmatrix <- diag(1,nrow=length(y), ncol=length(y)) - rho*W
logL <- ( -0.5*n*log(pi * sigma2) + log( det(lagmatrix) )
-(t(e) %*% e) / (2*sigma2) )
return(-logL)
}


Optimizing this function is far too slow to be practical, however. Unsurprisingly, the highest cost operation is the log-determinant.

1. Do you have tips for computing this determinant more quickly?
2. Does lagsarlm already have functionality to do this?
3. Can I hack it to make it do this?
• I have managed to cut some time off the determinant: because W is symmetric and real, I can use the Cholesky decomposition to cut the time of det() from 28 seconds to 11. Would still like to be faster. – gregmacfarlane Dec 13 '12 at 21:32
• What is 'SDM' in SDM Model? Also, the topic of this work highly interests me. Can you give a reference to these models that use inverse distance matrices? – hearse Dec 13 '12 at 21:35
• SDM-Spatial Durbin Model. The text I am working from is "Introduction to Spatial Econometrics" by LeSage and Pace, who do an excellent job of organizing and presenting theory developed largely by Luc Anselin. – gregmacfarlane Dec 13 '12 at 22:44

This turns out to be a two-part answer

### 1. Log-Determinant

The spdep package does contain tools to efficiently calculate the log-determinant of a sparse weights matrix. The following lines of code drop the calculation time from 11 seconds to about 2. The code is not very clean because the methods are not intended for use by typical statistical consumers.

# Calculate log-determinant of (I- rho*W) - using methods in the spdep package
env <- new.env(parent=globalenv())
assign("listw", Wlist, envir=env)
assign("can.sim", can.sim, envir=env)
assign("similar", FALSE, envir=env)
assign("family", "SAR", envir=env)
assign("n", n, envir=env)
Matrix_setup(env, Imult=2, super=FALSE)
get("similar", envir=env)
logdet <- do_ldet(rho, env)