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

sacsarlm(formula, data = list(), listw, listw2 = NULL, na.action, type="sac", method = "eigen", quiet = NULL, zero.policy = NULL, tol.solve = 1e-10, llprof=NULL, interval1=NULL, interval2=NULL, trs1=NULL, trs2=NULL, control = list())
a symbolic description of the model to be fit. The details of model specification are given for lm()
an optional data frame containing the variables in the model. By default the variables are taken from the environment which the function is called
a listw object created for example by nb2listw
a listw object created for example by nb2listw, if not given, set to the same spatial weights as the listw argument
a function (default options("na.action")), can also be na.omit or na.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 no-neighbour observations. Note that only weights lists created without using the glist argument to nb2listw may be subsetted.
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 row-standardised weights, but otherwise included
"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 log-determinant methods.
default NULL, use !verbose global option value; if FALSE, reports function values during optimization.
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
the tolerance for detecting linear dependencies in the columns of matrices to be inverted - passed to solve() (default=1.0e-10). This may be used if necessary to extract coefficient standard errors (for instance lowering to 1e-12), but errors in solve() 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
default NULL, can either be an integer, to divide the feasible ranges into a grid of points, or a two-column 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
list of extra control arguments - see section below

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.


A list object of class sarlm

Control arguments


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), 1-36.

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), 150-179.

See Also

lm, lagsarlm, errorsarlm, summary.sarlm, eigenw, impacts.sarlm

  • sacsarlm
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)
Documentation reproduced from package spdep, version 0.6-9, License: GPL (>= 2)

Community examples

Looks like there are no examples yet.