sacsarlm
Spatial simultaneous autoregressive SAC model estimation
Maximum likelihood estimation of spatial simultaneous autoregressive “SAC/SARAR” models of the form:
$$y = \rho W1 y + X \beta + u, u = \lambda W2 u + \varepsilon$$
where $rho$ and $lambda$ are found by nlminb
or optim()
first, and $beta$ and other parameters by generalized least squares subsequently
 Keywords
 spatial
Usage
sacsarlm(formula, data = list(), listw, listw2 = NULL, na.action, type="sac", method = "eigen", quiet = NULL, zero.policy = NULL, tol.solve = 1e10, llprof=NULL, interval1=NULL, interval2=NULL, trs1=NULL, trs2=NULL, control = list())
Arguments
 formula
 a symbolic description of the model to be fit. The details
of model specification are given for
lm()
 data
 an optional data frame containing the variables in the model. By default the variables are taken from the environment which the function is called
 listw
 a
listw
object created for example bynb2listw
 listw2
 a
listw
object created for example bynb2listw
, if not given, set to the same spatial weights as thelistw
argument  na.action
 a function (default
options("na.action")
), can also bena.omit
orna.exclude
with consequences for residuals and fitted values  in these cases the weights list will be subsetted to remove NAs in the data. It may be necessary to set zero.policy to TRUE because this subsetting may create noneighbour observations. Note that only weights lists created without using the glist argument tonb2listw
may be subsetted.  type
 default "sac", may be set to "sacmixed" for the Manski model to include the spatially lagged independent variables added to X using
listw
; when "sacmixed", the lagged intercept is dropped for spatial weights style "W", that is rowstandardised weights, but otherwise included  method
 "eigen" (default)  the Jacobian is computed as the product
of (1  rho*eigenvalue) using
eigenw
, and "spam" or "Matrix" for strictly symmetric weights lists of styles "B" and "C", or made symmetric by similarity (Ord, 1975, Appendix C) if possible for styles "W" and "S", using code from the spam or Matrix packages to calculate the determinant; "LU" provides an alternative sparse matrix decomposition approach. In addition, there are "Chebyshev" and Monte Carlo "MC" approximate logdeterminant methods.  quiet
 default NULL, use !verbose global option value; if FALSE, reports function values during optimization.
 zero.policy
 default NULL, use global option value; if TRUE assign zero to the lagged value of zones without
neighbours, if FALSE (default) assign NA  causing
sacsarlm()
to terminate with an error  tol.solve
 the tolerance for detecting linear dependencies in the columns of matrices to be inverted  passed to
solve()
(default=1.0e10). This may be used if necessary to extract coefficient standard errors (for instance lowering to 1e12), but errors insolve()
may constitute indications of poorly scaled variables: if the variables have scales differing much from the autoregressive coefficient, the values in this matrix may be very different in scale, and inverting such a matrix is analytically possible by definition, but numerically unstable; rescaling the RHS variables alleviates this better than setting tol.solve to a very small value  llprof
 default NULL, can either be an integer, to divide the feasible ranges into a grid of points, or a twocolumn matrix of spatial coefficient values, at which to evaluate the likelihood function
 trs1, trs2
 default NULL, if given, vectors for each weights object of powered spatial weights matrix traces output by
trW
; when given, used in some Jacobian methods  interval1, interval2
 default is NULL, search intervals for each weights object for autoregressive parameters
 control
 list of extra control arguments  see section below
Details
Because numerical optimisation is used to find the values of lambda and rho, care needs to be shown. It has been found that the surface of the 2D likelihood function often forms a “banana trench” from (low rho, high lambda) through (high rho, high lambda) to (high rho, low lambda) values. In addition, sometimes the banana has optima towards both ends, one local, the other global, and conseqently the choice of the starting point for the final optimization becomes crucial. The default approach is not to use just (0, 0) as a starting point, nor the (rho, lambda) values from gstsls
, which lie in a central part of the “trench”, but either four values at (low rho, high lambda), (0, 0), (high rho, high lambda), and (high rho, low lambda), and to use the best of these start points for the final optimization. Optionally, nine points can be used spanning the whole (lower, upper) space.
Value

A list object of class
sarlm
Control arguments
References
Anselin, L. 1988 Spatial econometrics: methods and models. (Dordrecht: Kluwer); LeSage J and RK Pace (2009) Introduction to Spatial Econometrics. CRC Press, Boca Raton.
Roger Bivand, Gianfranco Piras (2015). Comparing Implementations of Estimation Methods for Spatial Econometrics. Journal of Statistical Software, 63(18), 136. http://www.jstatsoft.org/v63/i18/.
Bivand, R. S., Hauke, J., and Kossowski, T. (2013). Computing the Jacobian in Gaussian spatial autoregressive models: An illustrated comparison of available methods. Geographical Analysis, 45(2), 150179.
See Also
lm
, lagsarlm
, errorsarlm
,
summary.sarlm
, eigenw
, impacts.sarlm
Examples
data(oldcol)
COL.sacW.eig < sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,
nb2listw(COL.nb, style="W"))
summary(COL.sacW.eig, correlation=TRUE)
W < as(nb2listw(COL.nb, style="W"), "CsparseMatrix")
trMatc < trW(W, type="mult")
summary(impacts(COL.sacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE)
COL.msacW.eig < sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,
nb2listw(COL.nb, style="W"), type="sacmixed")
summary(COL.msacW.eig, correlation=TRUE)
summary(impacts(COL.msacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE)