Skip to content

Instantly share code, notes, and snippets.

@devmag devmag/spatial.R
Last active Jan 30, 2018

Embed
What would you like to do?
R scripts to compute Conley Type Standard Errors
#####VECTORISED FUNCTION
library(data.table)
library(geosphere)
library(foreign)
library(lfe)
library(reshape)
iterateObs<-function(y1,e1,X1,fordist,coefficients,cutoff=250000) {
##recognise whether it is lat/long or single dimension (i.e. time) for distance computation
if(ncol(fordist)==2) {
distances<-lapply(1:nrow(X1), function(k) distHaversine(fordist[k,],as.matrix(fordist)))
XeeXhs<-lapply(1:nrow(X1), function(k) ( t(t(X1[k,])) %*% matrix(nrow=1,ncol=nrow(X1),data=e1[k]) * (matrix(nrow=length(coefficients),ncol=1,1)%*% (t(e1) * (distances[[k]]<=cutoff)))) %*% X1)
} else {
abstimediff <-lapply(1:nrow(fordist), function(k) abs(fordist[k]-fordist) )
window<-lapply(1:nrow(fordist), function(k) ((abstimediff[[k]] <= cutoff) * (1-abstimediff[[k]])/(cutoff+1)) * (fordist!=fordist[k]) )
XeeXhs<-lapply(1:nrow(X1), function(k) ( t(t(X1[k,])) %*% matrix(nrow=1,ncol=nrow(X1),data=e1[k]) * (matrix(nrow=length(coefficients),ncol=1,1)%*% (t(e1) * (t(window[[k]]))))) %*% X1)
}
XeeXhs <- Reduce("+", XeeXhs)
XeeXhs
}
###try to vectorise the function
conleyHAC<-function(y,yhat,X,coords,coefficients=NULL,timeid,panelid, dist_cutoff=200000,lag_cutoff=3) {
e<-y-yhat
XeeX = matrix(data=0, nrow=length(coefficients),ncol=length(coefficients))
times<-names(table(timeid))
n=length(y)
X<-as.matrix(X)
###spatial correlation correction
XeeXhs<-lapply(times, function(x) iterateObs(y[timeid==x],e[timeid==x],matrix(X[timeid==x,]),as.matrix(coords[timeid==x,]),coefficients,cutoff=dist_cutoff))
##first reduce
XeeX <- Reduce("+", XeeXhs)
invXX = solve(t(X)%*%X) * n
XeeX_spatial = XeeX / n
V = (invXX %*% XeeX_spatial %*% invXX) / n
V
#####serial correlation correction
panel<-names(table(panelid))
XeeXhs<-lapply(panel, function(x) iterateObs(y[panelid==x],e[panelid==x],matrix(X[panelid==x,]),matrix(timeid[panelid==x]),coefficients,cutoff=lag_cutoff))
XeeX <- XeeX+Reduce("+", XeeXhs)
XeeX_spatial_HAC = XeeX / n
V = (invXX %*% XeeX_spatial_HAC %*% invXX) / n
V
}
###function call
setwd("/Users/thiemo/AeroFS/Research/Fracking Growth/SPATIAL")
df<-data.table(read.dta(file="testspatial.dta"))
##download sample data from http://freigeist.devmag.net/wp-content/testspatial.dta_.zip
##create residuals
df.lm<-felm(formula = log(EmpClean) ~ post08shale | FIPS + year | (0) | STATE_FIPS, data = df)
###CONLEY ERRORS for COEFFICIENT post08shale
conleyHAC(df.lm$response,df.lm$fitted.values,df$post08anywell,df[,c("latitude","longitude"),with=F],coefficients=df.lm$coefficients,timeid=df$year,panelid=df$FIPS,dist_cutoff=10000,lag_cutoff=10)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.