The primary motivation for the sarml command is to provide a convenient way to estimate multiple spatial AR models
without having to calculate the eigenvalues of W each time. Under the assumption that the errors, u, are
independently and identically distributed normal, the log-likelihood function for the spatial AR model is
$$lnl = -\frac{n}{2}log(\pi) - \frac{n}{2}log(\sigma^2) - \frac{1}{2 \sigma^2}\sum_i u_i^2 + \sum_ilog(1 - \rho*eigvar_i)$$
where eigvar is the vector of eigenvalues of W.
Though spdep provides a convenient and fast method for calculating the eigenvalues from a shape file,
the calculation can nonetheless take a very long time for large data sets.
The sarml command allows the user to input the vector of eigenvalues directly,
which saves time when several models are estimated using the same W matrix.
Unless a vector of eigenvalues is provided using the eigvar option,
the eigenvalues are calculated from the shape file (provided using the shpfile option) using the spdep package.
Conditional on the value of $\rho$, the maximum likelihood estimate of $\beta$ is simply the vector of coefficients from a regression of $Y - \rho WY$ on X.
The estimate of the error variance also has a closed form solution: $\sigma^2 = \sum u^2/n$.
Substituting these estimates into the log-likelihood functions leads to the following concentrated log-likelihood function:
$$lc = -\frac{n}{2}(log(\pi)+1) - \frac{n}{2}log(\sum_iu_i^2) + \sum_i log(1-\rho*eigvar_i)$$
Working with the concentrated likelihood function reduces the optimization problem to a one-dimensional search for the value of $\rho$ that maximizes lc.
Unless a value is provided for startrho, the sarml procedure begins by using the optimize command to find the value of $\rho$
that maximizes lc.
This estimate of $\rho$ (or the value provided by the startrho option) is then used to calculate the implied values of $\beta$ and $\sigma^2$,
and these values are used as starting values to maximize lnl using the nlm command.
The covariance matrix for the estimates of $\beta$ and $\rho$, vmat, is the inverse of $(1/\sigma^2)V$.
V has partitions $V11 = X'X$, $V12 = X'WY$, $V21 = Y'W'X$,
and $V22 = Y'W'WY + s2*\sum (eigvar/(1-rho*eigvar))^2$.