ML estimate of neighbors 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$$

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]
 radius <- params[2]
 sigma2 <- params[3]
 delta <- params[-c(1,2,3)]
 # make W subject to radius
 points.dnn <- dnearneigh(GeoPoints, 0, radius*5280)
 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?
  2.a). Can I hack it to make it do this?

Greg Macfarlane
Georgia Tech | Civil Engineering


	[[alternative HTML version deleted]]