McSpatial (version 1.1.1)

sarml: Spatial AR Maximum-Likelihood Estimation

Description

Estimates the model $Y = \rho WY + X \beta + u$ by maximizing the log-likelihood function.

Usage

sarml(form,wmat=NULL,shpfile=NULL,wy=NULL,eigvar=NULL,startrho=NULL,
  print=TRUE,data=NULL)

Arguments

form

Model formula

wmat

The W matrix. If not specified, W will be calculated from the shape file. Default: W = NULL.

shpfile

Shape file. Needed unless (1) wy and eigvar are both provided, or (2) wmat and eigvar are provided

wy

The WY variable. Default: not specified; program attempts to calculate WY using wmat or shpfile.

eigvar

The vector of eigenvalues for W. Default: not provided. shpfile must be specified to calculate the eigenvalues within the sarml command.

startrho

A starting value for $\rho$. Default: startrho=0. Estimation will generally be faster if startrho = 0.

print

If print=F, no results are printed. Default: print=T.

data

A data frame containing the data. Default: use data in the current working directory

Value

• beta
  The estimated vector of coefficients, $\beta$.
• rho
  The estimated value of $\rho$.
• sig2
  The estimated error variance, $\sigma^2$.
• vmat
  The covariance matrix for $(\beta, \sigma^2)$.
• eigvar
  The vector of eigenvalues

Details

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: $\hat{\sigma^2} = \sum_iu_i^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 $V_{11} = X'X$, $V_{12} = X'WY$, $V_{21} = Y'W'X$, and $V_{22} = Y'W'WY + \sigma^2 \sum_i \frac{eigvar_i}{(1 - \rho*eigvar_i)^2 }$.

See Also

makew qregspiv qregbmat qregsim1 qregsim2 qregcpar qreglwr

Examples

library(spdep)

cmap <- readShapePoly(system.file("maps/CookCensusTracts.shp",
  package="McSpatial"))
cmap <- cmap[cmap$CHICAGO==1&cmap$CAREA!="O'Hare",]
samppop <- cmap$POPULATION>0&cmap$AREA>0
cmap <- cmap[samppop,]
cmap$lndens <- log(cmap$POPULATION/cmap$AREA)
lmat <- coordinates(cmap)
cmap$LONGITUDE <- lmat[,1]
cmap$LATITUDE  <- lmat[,2]
cmap$dcbd <- geodistance(longvar=cmap$LONGITUDE,latvar=cmap$LATITUDE,
  lotarget=-87.627800,latarget=41.881998)$dist

fit <- makew(shpfile=cmap,eigenvalues=TRUE)
wmat <- fit$wmat
eigvar <- fit$eigvar

# input w, calculate eigvar within sarml
fit <- sarml(lndens~dcbd,wmat=wmat,eigvar=eigvar,data=cmap)