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
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
- an optional data frame containing the variables in the model. By default the variables are taken from the environment which the function is called
listwobject created for example by
listwobject created for example by
nb2listw, if not given, set to the same spatial weights as the
- a function (default
options("na.action")), can also be
na.excludewith 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
nb2listwmay 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
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. 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), 150-179.
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)