spsurml: Maximum likelihood estimation of spatial SUR model.

Description

This function estimates spatial SUR models using maximum-likelihood methods.The number of equations, time periods and cross-sectional units is not restricted.The user can choose between different spatial specifications as described below. The estimation procedure allows for the introduction of linear restrictions on the $\beta$ parameters associated to the regressors.

Usage

spsurml(formula = NULL, data = NULL, na.action,
               listw = NULL, type = "sim", Durbin = NULL,
               method = "eigen", zero.policy = NULL, interval = NULL,
               trs = NULL, R = NULL, b = NULL, X = NULL, Y = NULL, 
               G = NULL, N = NULL, Tm = NULL,p = NULL, 
               control = list() )

Value

Object of spsur class with the output of the maximum-likelihood estimation of the specified spatial SUR model. A list with:

`call`	Matched call.
`type`	Type of model specified.
`method`	Value of `method` argument to compute the Jacobian
`Durbin`	Value of `Durbin` argument.
`coefficients`	Estimated coefficients for the regressors.
`deltas`	Estimated spatial coefficients.
`rest.se`	Estimated standard errors for the estimates of beta.
`deltas.se`	Estimated standard errors for the estimates of the spatial coefficients (`deltas`).
`resvar`	Estimated covariance matrix for the estimates of beta's and spatial coefficients (`deltas`).
`LL`	Value of the likelihood function at the maximum-likelihood estimates.
`R2`	Coefficient of determination for each equation, obtained as the squared of the correlation coefficient between the corresponding explained variable and its estimate. `spsurml` also shows a global coefficient of determination obtained, in the same manner, for the set of the G equations.
`Sigma`	Estimated covariance matrix for the residuals of the G equations.
`fdHess`	Logical value of `fdHess` argument when computing numerical covariances.
`residuals`	Residuals of the model.
`df.residuals`	Degrees of freedom for the residuals.
`fitted.values`	Estimated values for the dependent variables.
`BP`	Value of the Breusch-Pagan statistic to test the null hypothesis of diagonality among the errors of the G equations.
`LMM`	Marginal Lagrange Multipliers, LM($\rho$\|$\lambda$) and LM($\lambda$\|$\rho$), to test for omitted spatial effects in the specification.
`G`	Number of equations.
`N`	Number of cross-sections or spatial units.
`Tm`	Number of time periods.
`p`	Number of regressors by equation (including intercepts).
`Y`	Vector Y of the explained variables of the SUR model.
`X`	Matrix X of the regressors of the SUR model.
`W`	Spatial weighting matrix.
`zero.policy`	Logical value of `zero.policy` .
`interval`	Search interval for spatial parameter.
`listw_style`	Style of neighborhood matrix `W`.
`trs`	Either `NULL` or vector of powered spatial weights matrix traces output by `trW`.
`insert`	Logical value to check if `is.null(trs)`.

Control arguments

`tol`	Numerical value for the tolerance for the estimation algorithm until convergence. Default = 1e-3.
`maxit`	Maximum number of iterations until convergence; it must be an integer value. Default = 200.
`trace`	A logical value to show intermediate results during the estimation process. Default = `TRUE`.
`fdHess`	Compute variance-covariance matrix using the numerical hessian. Suited for large samples. Default = `FALSE`
`Imult`	default 2; used for preparing the Cholesky decompositions for updating in the Jacobian function
`super`	if `NULL` (default), set to `FALSE` to use a simplicial decomposition for the sparse Cholesky decomposition and method "Matrix_J", set to as.logical(NA) for method "Matrix", if `TRUE`, use a supernodal decomposition
`cheb_q`	default 5; highest power of the approximating polynomial for the Chebyshev approximation
`MC_p`	default 16; number of random variates
`MC_m`	default 30; number of products of random variates matrix and spatial weights matrix
`spamPivot`	default "MMD", alternative "RCM"
`in_coef`	default 0.1, coefficient value for initial Cholesky decomposition in "spam_update"
`type`	default "MC", used with method "moments"; alternatives "mult" and "moments", for use if trs is missing
`correct`	default `TRUE`, used with method "moments" to compute the Smirnov/Anselin correction term
`trunc`	default `TRUE`, used with method "moments" to truncate the Smirnov/Anselin correction term
`SE_method`	default "LU", may be "MC"
`nrho`	default 200, as in SE toolbox; the size of the first stage lndet grid; it may be reduced to for example 40
`interpn`	default 2000, as in SE toolbox; the size of the second stage lndet grid
`SElndet`	default `NULL`, may be used to pass a pre-computed SE toolbox style matrix of coefficients and their lndet values to the "SE_classic" and "SE_whichMin" methods
`LU_order`	default `FALSE`; used in "LU_prepermutate", note warnings given for lu method
`pre_eig`	default `NULL`; may be used to pass a pre-computed vector of eigenvalues

Details

The list of (spatial) models that can be estimated with the spsurml function are:

"sim": SUR model with no spatial effects $$ y_{tg} = X_{tg} \beta_{g} + \epsilon_{tg} $$
"slx": SUR model with spatial lags of the regressors $$ y_{tg} = X_{tg} \beta_{g} + WX_{tg} \theta_{g} + \epsilon_{tg} $$
"slm": SUR model with spatial lags of the explained variables $$y_{tg} = \rho_{g} Wy_{tg} + X_{tg} \beta_{g} + \epsilon_{tg} $$
"sem": SUR model with spatial errors $$ y_{tg} = X_{tg} \beta_{g} + u_{tg} $$ $$ u_{tg} = \lambda_{g} Wu_{tg} + \epsilon_{tg} $$
"sdm": SUR model of the Spatial Durbin type $$ y_{tg} = \rho_{g} Wy_{tg} + X_{tt} \beta_{g} + WX_{tg} \theta_{g} + \epsilon_{tg} $$
"sdem": SUR model with spatial lags of the regressors and spatial errors $$ y_{tg} = X_{tg} \beta_{g} + WX_{tg} \theta_{g} + u_{tg} $$ $$ u_{tg} = \lambda_{g} W u_{tg} + \epsilon_{tg} $$
"sarar": SUR model with spatial lags of the explained variables and spatial errors $$ y_{tg} = \rho_{g} Wy_{tg} + X_{tg} \beta_{g} + u_{tg} $$ $$ u_{tg} = \lambda_{g} W u_{tg} + \epsilon_{tg} $$
"gnm": SUR model with spatial lags of the explained variables, regressors and spatial errors $$ y_{tg} = \rho_{g} Wy_{tg} + X_{tg} \beta_{g} + WX_{tg} \theta_{g} + u_{tg} $$ $$ u_{tg} = \lambda_{g} W u_{tg} + \epsilon_{tg} $$

References

Anselin, L. (1988). Spatial econometrics: methods and models. Dordrecht: Kluwer
Bivand, R.S. and Piras G. (2015). Comparing Implementations of Estimation Methods for Spatial Econometrics. Journal of Statistical Software, 63(18), 1-36. https://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.
Breusch T, Pagan A (1980). The Lagrange multiplier test and its applications to model specification in econometrics. Rev Econ Stud 47: 239-254
Cliff, A. D. and Ord, J. K. (1981). Spatial processes: Models and applications, Pion.
LeSage J and Pace, R.K. (2009). Introduction to Spatial Econometrics. CRC Press, Boca Raton.
L<U+00F3>pez, F.A., Mur, J., and Angulo, A. (2014). Spatial model selection strategies in a SUR framework. The case of regional productivity in EU. Annals of Regional Science, 53(1), 197-220.
Mur, J., L<U+00F3>pez, F., and Herrera, M. (2010). Testing for spatial effects in seemingly unrelated regressions. Spatial Economic Analysis, 5(4), 399-440.
Ord, J. K. (1975). Estimation methods for models of spatial interaction, Journal of the American Statistical Association, 70, 120-126;

Examples

Run this code

# NOT RUN {
#################################################
######## CROSS SECTION DATA (G>1; Tm=1) ########
#################################################

#### Example 1: Spatial Phillips-Curve. Anselin (1988, p. 203)
rm(list = ls()) # Clean memory
data(spc)
Tformula <- WAGE83 | WAGE81 ~ UN83 + NMR83 + SMSA | UN80 + NMR80 + SMSA
spcsur.sim <-spsurml(formula = Tformula, data = spc, type = "sim")
summary(spcsur.sim)
# All the coefficients in a single table.
print(spcsur.sim)
# Plot of the coefficients of each equation in different graphs
plot(spcsur.sim) 

## A SUR-SLX model 
## (listw argument can be either a matrix or a listw object )
spcsur.slx <-spsurml(formula = Tformula, data = spc, type = "slx", 
  listw = Wspc)
summary(spcsur.slx)
# All the coefficients in a single table.
print(spcsur.slx)
# Plot of the coefficients in a single graph
if (require(gridExtra)) {
  pl <- plot(spcsur.slx, viewplot = FALSE)
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
                        nrow = 2)
} 

## VIP: The output of the whole set of the examples can be examined 
## by executing demo(demo_spsurml, package="spsur")
  
# }
# NOT RUN {
### A SUR-SLM model
spcsur.slm <-spsurml(formula = Tformula, data = spc, type = "slm", 
                     listw = Wspc)
summary(spcsur.slm)
if (require(gridExtra)) {
  pl <- plot(spcsur.slm, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

#' ### A SUR-SEM model
spcsur.sem <-spsurml(formula = Tformula, data = spc, type = "sem", 
                     listw = Wspc)
summary(spcsur.sem)                      
print(spcsur.sem)
if (require(gridExtra)) {
  pl <- plot(spcsur.sem, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

### A SUR-SDM model
spcsur.sdm <-spsurml(formula = Tformula, data = spc, type = "sdm", 
                     listw = Wspc)
summary(spcsur.sdm)
print(spcsur.sdm)
if (require(gridExtra)) {
  pl <- plot(spcsur.sdm, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

## A SUR-SDM model with different spatial lags in each equation
TformulaD <- ~ UN83 + NMR83 + SMSA | UN80
spcsur.sdm2 <-spsurml(formula = Tformula, data = spc, type = "sdm", 
                      listw = Wspc, Durbin = TformulaD)
summary(spcsur.sdm2)
if (require(gridExtra)) {
  pl <- plot(spcsur.sdm2, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}
### A SUR-SDEM model
spcsur.sdem <-spsurml(formula = Tformula, data = spc, type = "sdem", 
                      listw = Wspc)
print(spcsur.sdem)
if (require(gridExtra)) {
  pl <- plot(spcsur.sdem, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

### A SUR-SARAR model
spcsur.sarar <-spsurml(formula = Tformula, data = spc, type = "sarar", 
                       listw = Wspc, control = list(tol = 0.1))
print(spcsur.sarar)
if (require(gridExtra)) {
  pl <- plot(spcsur.sarar, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

### A SUR-GNM model
spcsur.gnm <-spsurml(formula = Tformula, data = spc, type = "gnm", 
                       listw = Wspc, control = list(tol = 0.1))
print(spcsur.gnm)
if (require(gridExtra)) {
  pl <- plot(spcsur.gnm, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

## A A SUR-GNM model model with different spatial lags in each equation
TformulaD <- ~ UN83 + NMR83 + SMSA | UN80
spcsur.gnm2 <-spsurml(formula = Tformula, data = spc, type = "gnm", 
                       listw = Wspc, Durbin = TformulaD,
                       control = list(tol = 0.1))
print(spcsur.gnm2)
if (require(gridExtra)) {
  pl <- plot(spcsur.gnm2, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

# }
# NOT RUN {
##################################################
#########  G=1; Tm>1         ########
##################################################
#
##### Example 2: Homicides + Socio-Economics (1960-90)
## Homicides and selected socio-economic characteristics for continental
## U.S. counties.
## Data for four decennial census years: 1960, 1970, 1980 and 1990.
## \url{https://geodacenter.github.io/data-and-lab/ncovr/}

# }
# NOT RUN {
### It usually requires 1-2 minutes maximum...
rm(list = ls()) # Clean memory
### Read NCOVR.sf object
data(NCOVR, package = "spsur")
nbncovr <- spdep::poly2nb(NCOVR.sf, queen = TRUE)
### Some regions with no links...
lwncovr <- spdep::nb2listw(nbncovr, style = "W", zero.policy = TRUE)
Tformula <- HR80  | HR90 ~ PS80 + UE80 | PS90 + UE90
### A SUR-SIM model
NCOVRSUR.sim <-spsurml(formula = Tformula, data = NCOVR.sf, type = "sim")
summary(NCOVRSUR.sim)
if (require(gridExtra)) {
  pl <- plot(NCOVRSUR.sim, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], nrow = 3)
}
### A SUR-SLX model
NCOVRSUR.slx <-spsurml(formula = Tformula, data = NCOVR.sf, type = "slx", 
                       listw = lwncovr, zero.policy = TRUE)
print(NCOVRSUR.slx)
if (require(gridExtra)) {
  pl <- plot(NCOVRSUR.slx, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], nrow = 2)
}

### A SUR-SLM model
### method = "Matrix" (Cholesky) instead of "eigen"
### (fdHess = TRUE to compute numerical covariances )
NCOVRSUR.slm <-spsurml(formula = Tformula, data = NCOVR.sf, 
                       type = "slm", listw = lwncovr, method = "Matrix", 
                       zero.policy = TRUE, control = list(fdHess = TRUE))
summary(NCOVRSUR.slm)

if (require(gridExtra)) {
  pl <- plot(NCOVRSUR.slm, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

# LR test for nested models
#lrtest(NCOVRSUR.sim, NCOVRSUR.slm)

### A SUR-SDM model with different spatial lags in each equation
### Analytical covariances (default)
TformulaD <- ~ PS80 + UE80 | PS90 
NCOVRSUR.sdm <-spsurml(formula = Tformula, data = NCOVR.sf, 
                       type = "sdm", listw = lwncovr, method = "Matrix",
                       Durbin = TformulaD, zero.policy = TRUE)
print(NCOVRSUR.sdm)
if (require(gridExtra)) {
  pl <- plot(NCOVRSUR.sdm, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}
### A SUR-SEM model
NCOVRSUR.sem <- spsurml(formula = Tformula, data = NCOVR.sf, 
                        type = "sem", listw = lwncovr, method = "Matrix",
                        zero.policy = TRUE, control = list(fdHess = TRUE))
print(NCOVRSUR.sem)
if (require(gridExtra)) {
  pl <- plot(NCOVRSUR.sem, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}

### A SUR-SDEM model
NCOVRSUR.sdem <-spsurml(formula = Tformula, data = NCOVR.sf, 
                        type = "sdem", listw = lwncovr, method = "Matrix",
                        zero.policy = TRUE, control = list(fdHess = TRUE))
print(NCOVRSUR.sdem)
if (require(gridExtra)) {
  pl <- plot(NCOVRSUR.sdem, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$pldeltas, nrow = 3)
}
# }
# NOT RUN {
###############################################
######### SUR with G>1; Tm>1  ###
###############################################
##### Reshape NCOVR in panel format
# }
# NOT RUN {
N <- nrow(NCOVR.sf)
Tm <- 4
index_time <- rep(1:Tm, each = N)
index_indiv <- rep(1:N, Tm)
pHR <- c(NCOVR.sf$HR60, NCOVR.sf$HR70, NCOVR.sf$HR80, NCOVR.sf$HR90)
pPS <- c(NCOVR.sf$PS60, NCOVR.sf$PS70, NCOVR.sf$PS80, NCOVR.sf$PS90)
pUE <- c(NCOVR.sf$UE60, NCOVR.sf$UE70, NCOVR.sf$UE80, NCOVR.sf$UE90)
pDV <- c(NCOVR.sf$DV60, NCOVR.sf$DV70, NCOVR.sf$DV80, NCOVR.sf$DV90)
pFP <- c(NCOVR.sf$FP59, NCOVR.sf$FP70, NCOVR.sf$FP80, NCOVR.sf$FP90)
pSOUTH <- rep(NCOVR.sf$SOUTH, Tm)
pNCOVR <- data.frame(indiv = index_indiv, time = index_time,
                     HR = pHR, PS = pPS, UE = pUE, DV = pDV,
                     FP = pFP, SOUTH = pSOUTH)
pform <- HR | DV | FP ~ PS + UE | PS + UE + SOUTH | PS
### SIM 
### Remark: It is necessary to provide Tm value as argument 
### when G>1 && Tm>1
pNCOVRSUR.sim <- spsurml(formula = pform, data = pNCOVR, 
                         type = "sim", Tm = Tm)
print(pNCOVRSUR.sim)
if (require(gridExtra)) {
  pl <- plot(pNCOVRSUR.sim, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$lplbetas[[3]], nrow = 3)
}
# SLM
pNCOVRSUR.slm <- spsurml(formula = pform, data = pNCOVR, 
                         listw = lwncovr, type = "slm", method = "Matrix", Tm = Tm,
                         zero.policy = TRUE, control= list(fdHess = TRUE))
print(pNCOVRSUR.slm)
if (require(gridExtra)) {
  pl <- plot(pNCOVRSUR.slm, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$lplbetas[[3]], pl$pldeltas, nrow = 4)
}

pNCOVRSUR.sem <- spsurml(formula = pform, data = pNCOVR, 
                         listw = lwncovr, type = "sem", method = "Matrix", Tm = Tm,
                         zero.policy = TRUE, control= list(fdHess = TRUE))
print(pNCOVRSUR.sem) 
if (require(gridExtra)) {
  pl <- plot(pNCOVRSUR.sem, viewplot = FALSE) 
  grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]], 
               pl$lplbetas[[3]], pl$pldeltas, nrow = 4)
}
# }

Run the code above in your browser using DataLab