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

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?
  • $\begingroup$ 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. $\endgroup$ – gregmacfarlane Dec 13 '12 at 21:32
  • $\begingroup$ 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? $\endgroup$ – hearse Dec 13 '12 at 21:35
  • 2
    $\begingroup$ 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. $\endgroup$ – 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)

2. Maximum Likelihood Radius

Roger Bivand, the package author, explained to me that such a method is not really valid, and I discovered this myself when I couldn't find any traction on the likelihood function. What he instead recommended was to simply estimate the same model specification over a grid of radii and select the radius that provided the maximum point. The figure below shows how the likelihood value changes for higher sample sizes and different radii on my dataset.

Likelihood versus radius as sample size increases


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.