// Using Sp commands in Stata 15

// written by Robert Grant at BayesCamp for Timberlake Consultants, August 2017.

clear all
version 15

grmap, activate
cd "~/Dropbox/BayesCamp/clients/Timberlake/marketing/spatial"

// We are going to attach a local authority shapefile to a file of data on CO2 emissions

unzipfile "Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain.zip", replace
spshape2dta Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain, replace
// check it
use "Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain.dta", clear
describe
assert lad16cd!=""
bysort lad16cd: assert _N==1 // The variable lad16cd is the identifier. 
encode lad16cd, gen(lacode)
spset lacode, modify replace
save, replace

// now we open the carbon dioxide data
use "LA-CO2-2014.dta", clear
merge 1:1 lad16cd using "Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain.dta"

keep if _merge==3 // Sorry Northern Ireland! Just keeping it simple.
drop _merge
save "CO2-merged.dta", replace

import excel "File_10_ID2015_Local_Authority_District_Summaries.xlsx", sheet("IMD") firstrow clear
keep LocalAuthorityDistrictcode2 IMDAveragescore
rename LocalAuthorityDistrictcode2 lad16cd
merge 1:1 lad16cd using "CO2-merged.dta"
keep if _merge==3
drop _merge
grmap domestictotal, clnumber(9) fcolor(PuRd) name(dommap, replace) title("Total CO2")

// generate adjacency matrix
spmatrix create contiguity W, rook
spmatrix create idistance W2

// take logarithms of domestictotal and population
generate logdom = log(domestictotal)
generate logpop = log(population)

// make per capita CO2
generate pc = logdom-logpop
grmap pc, clnumber(9) fcolor(PuRd) name(pcmap, replace) title("Per capita CO2")

// exploratory graphs
twoway (scatter pc IMDA, msym(Oh)) (lowess pc IMDA) (lfit pc IMDA), name(pc_imd, replace)
twoway (scatter logdom IMDA, msym(Oh)) (lowess logdom IMDA) (lfit logdom IMDA), name(logdom_imd, replace)
twoway (scatter logdom logpop, msym(Oh)) (lowess logdom logpop) (lfit logdom logpop), name(logdom_logpop, replace)

// Linear regression model
regress pc IMDAveragescore
predict regres, residuals
scatter regres domest, msymbol(Oh) name(regres, replace)
grmap regres, clnumber(9) fcolor(BuYlRd) name(regresmap, replace)

// SAR model with adjacency weighting
spregress pc IMDAveragescore, gs2sls dvarlag(W)
predict spres, residuals
scatter spres domest, msymbol(Oh) name(spres, replace)
grmap spres, clnumber(9) fcolor(BuYlRd) name(spresmap, replace) title("SAR (adjacency) residuals")
twoway (scatter spres regres, ms(Oh) ///
			ytitle("SAR (adjacency) residuals") xtitle("Linear regression residuals")) ///
	   (function y=x, range(regres)), legend(off) name(resvres, replace)
	   
// SAR model with inverse distance weighting
spregress pc IMDAveragescore, gs2sls dvarlag(W2)
predict spres2, residuals
scatter spres2 domest, msymbol(Oh) name(spres2, replace)
grmap spres2, clnumber(9) fcolor(BuYlRd) name(spresmap2, replace) title("SAR (inverse-distance) residuals")
twoway (scatter spres2 regres, ms(Oh) ///
			ytitle("SAR (inverse-distance) residuals") xtitle("Linear regression residuals")) ///
	   (function y=x, range(regres)), legend(off) name(resvres2, replace)