# Quant 3: Section on spatial econometrics
# Alex Herzog
# Dec 7, 2007


# Note: R Code in this file is (partially) taken from Kristian Skrede Gleditsch's 
# webpage "Supplementary material for Spatial Regression Models"
# (http://privatewww.essex.ac.uk/~ksg/srm_book.html)

# Data on democracy and development is from the same website as above. 

# Install Roger Bivand's (http://www.nhh.no/sam/cv/bivand-roger.html)
# `spdep' package
# install.packages("spdep", dep=TRUE) 

# load library
library(spdep)

# Set working directory
#setwd("PATH-TO-YOUR-FILES")

# Democracy and GDP data from chapter 1 of the book
source("spatial_data.Rdata")

objects()
summary(sldv)


### MOTIVATING EXAMPLE ###

# Estimate OLS of democracy on gdp
dem <- sldv$democracy
log.gdp <- log(sldv$gdp.2002/sldv$population)

ols.out <- lm(dem ~ log.gdp)
summary(ols.out)  

# Plot of democracy versus gdp together with regression line
plot(log.gdp,dem,abline(glm(dem ~ log.gdp)))

# Distribution of residuals
resid <- ols.out$residuals
plot(density(resid))





### SPATIAL OBJECTS ###

# What data do we have?
# wmat: row normalized connectivty matrix based on a 200 km threshold
round(wmat,digits=2)[1:10,1:10]
rowSums(wmat)

# nblist: neighbor list; a list of connectivities based on a 200 km threshold
nblist	
attributes(nblist)
summary(nblist)

length(nblist)
nblist[[1]]
nblist[[1]][2]
nblist[sldv$tla=="UKG"]

# rs.nblist: neighbors list with row standardized spatial weights
rs.nblist <- nb2listw(nblist)



### SPATIAL LAG VARIABLE AND TEST FOR SPATIAL CORRELATION ###

# calculate spatial lag of democracy (average value of neighbors democracy)
dem.sl <- as.vector(wmat %*% dem)

wmat.uk <- wmat[rownames(wmat)=="UKG",]
wmat.uk[wmat.uk != 0]

dem[wmat.uk != 0]
mean(dem[wmat.uk != 0])

dem.sl[rownames(wmat)=="UKG"]


# Moran scatterplot: y against the spatial lag variable
# i) `by hand'
plot(dem,dem.sl,abline(glm(dem.sl ~ dem)))
abline(v=mean(dem),lty=3)
abline(h=mean(dem.sl),lty=3)

	# Moran's I test for spatial autocorrelation using a spatial weights matrix in weights list form
	summary(glm(dem.sl ~ dem))
	moran.test(dem,rs.nblist)


# ii) with `moran.plot()'
moran.plot(dem,rs.nblist)


# iii) with standardized variables
dem.std <- (dem - mean(dem)) / sqrt(var(dem))
dem.sl.std <- (dem.sl - mean(dem.sl)) / sqrt(var(dem.sl))

plot(dem.std,dem.sl.std,
	xlim=c(-2,2), ylim=c(-2,2),
	abline(glm(dem.sl.std ~ dem.std)))
abline(v=mean(dem.std),lty=3)
abline(h=mean(dem.sl.std),lty=3)



### SPATIALLY LAGGED DEPENDENT VARIABLE MODEL AND SPATIAL ERROR MODEL ###
ols.out2 <- glm(dem ~ log.gdp + dem.sl)

ols.out$coef
ols.out2$coef

# Estimation of the spatial simultaneous autoregressive lag model
sldvm.out <- lagsarlm(dem ~ log.gdp,rs.nblist,data=sldv)
summary(sldvm.out)

# Estimation of the spatial error model
sem.out <- errorsarlm(dem ~ log.gdp,rs.nblist,data=sldv)
summary(sem.out)

print(ols.out2$coef)
print(c(sldvm.out$coef,sldvm.out$rho))
print(c(sem.out$coef,sem.out$lambda))