Added distance matrix in temp.r

temp.r

@@ -1,17 +1,88 @@
 library(tidyverse)
 library(fixest)
+library("geosphere")
+
 DF <- readRDS("Data/Output_Data/Full_Data_Set.rds") %>% mutate(year=as.numeric(year)) %>% group_by(wdid) %>% mutate(SBD1=max(SBD1)) %>% ungroup
+EMPTY_START <- DF %>% filter(year<=2010) %>% group_by(wdid) %>% filter(sum(AF)==0) %>% select(wdid) %>% unique
+#DF <- DF %>% anti_join(EMPTY_START)
+DF[DF$CREP_any==1,'SBD1'] <- 1
+DF[DF$CREP_any==1 &DF$SBD1_year==Inf,'SBD1_year'] <- 2011
+
 DF$SBD1 <- ifelse(DF$SBD1_year!=Inf,1,0)
+#DF <- DF %>% select(-ditch_2200627,-ditch_2200541,-ditch_3500570)
 DF$SBD1_year <- ifelse(DF$SBD1_year==2009,2011,DF$SBD1_year)
 DF$POST <- ifelse(DF$year>=DF$SBD1_year,1,0)
-feols(AF~SBD1*POST+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD2+SBD3+SBD4+SBD5+SBD6+wdid+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF)
+DF <- DF %>% left_join(DF %>% group_by(wdid,CREP_any) %>% summarize(CREP_year=min(year)) %>% filter(CREP_any==1) %>% select(-CREP_any)) %>% mutate(CREP_year=ifelse(is.na(CREP_year),Inf,CREP_year))
+AUG <- read_csv("Data/Augmentation/Augmentation.csv") %>% mutate(wdid=as.character(wdid)) %>% select(wdid)
+DF <- DF %>% anti_join(AUG)
+#feols(AF~SBD1*POST+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD2+SBD3+SBD4+SBD5+SBD6+wdid+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF)
 ###########
 #DF$SBD1_year <- ifelse(DF$SBD1_year==Inf,10000,DF$SBD1_year)
-SUN_MOD <- feols(AF~sunab(SBD1_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_alfalfa+per_potatoes|SBD1+SBD2+SBD3+SBD4+SBD5+SBD6+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF)
 
-SUN_MOD <- feols(AF~sunab(SBD1_year,year)|wdid+SBD1+SBD2+SBD3+SBD4+SBD5+SBD6+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF %>% filter(SBD1_year==2011|SBD1_year==Inf))
+
-etable(SUN_MOD,agg='cohort')
+
-coefplot(SUN_MOD)
+
-#######CHECK COHORTS IN HYDROBASE LOOKS LIKE MISSING DATA IN SOME LATER COHORTS. DATA MAY NEED WORK
+SUN_MOD <- feols(AF~SBD1+sunab(SBD1_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_alfalfa+per_potatoes|wdid+year+SBD2^year+SBD3^year+SBD4^year+SBD5^year+SBD6^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year,data=DF %>% filter(SBD1_year==Inf|SBD1_year==2011) )
-GRAPH_DATA <- DF %>% group_by(SBD1_year,year) %>% summarize(AF=mean(AF)/max(AF))
+#coefplot(SUN_MOD)
-ggplot(GRAPH_DATA,aes(x=year,y=AF,color=as.factor(SBD1_year),group=SBD1_year))+geom_point()+geom_line()
+#coefplot(feols(AF~sunab(CREP_year,year)|wdid+year,data=DF))
+coefplot(feols(AF~sunab(CREP_year,year)|wdid+year+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF))
+
+
+
+
+
+
+GET_PROGRAM <-function(PROGRAM_NUMBER,DATA,CUTOFF=0.25,YEAR_RANGE=2009:2025){
+#DATA <- DF
+#PROGRAM_NUMBER <- 4
+#YEAR_RANGE=2009:2025
+#CUTOFF <- 0.25
+	DATA_ORIG <- DATA
+	MASTER_WELL_LIST <- DATA %>% select(wdid,longitude,latitude) %>% unique
+	WELLS_IN_PROGRAM <-	DATA[DATA[,PROGRAM_NUMBER]==1,]
+	PROGRAM_START <- min(WELLS_IN_PROGRAM$year)
+	PROGRAM_YEAR_RANGE <- min(WELLS_IN_PROGRAM$year):max(YEAR_RANGE)
+	for(YEAR in PROGRAM_YEAR_RANGE){
+#		YEAR <-2014
+	C_PROGRAM_WELLS <-	WELLS_IN_PROGRAM[WELLS_IN_PROGRAM[,"year"]==YEAR,c("longitude","latitude")] %>% unique %>% as.matrix
+RES_DISTANCE <- c()
+for(i in 1:nrow(MASTER_WELL_LIST)){
+	RES_DISTANCE[i] <-ifelse((2.38084*min(distHaversine(MASTER_WELL_LIST[i,c("longitude","latitude")],C_PROGRAM_WELLS))/5280)<=CUTOFF,1,0)
+}
+	C_RES <- cbind(MASTER_WELL_LIST[,1],RES_DISTANCE) %>% as_tibble
+	colnames(C_RES) <-c("wdid",paste0("close_",colnames(DATA_ORIG)[PROGRAM_NUMBER]))
+	C_RES$year <- YEAR
+	if(exists("RES")){RES <- rbind(RES,C_RES)} else{RES <- C_RES}
+}
+	if(PROGRAM_START>min(YEAR_RANGE)){
+FILL_IN_YEARS <- min(YEAR_RANGE):(PROGRAM_START-1)
+length(FILL_IN_YEARS )
+5*(nrow(MASTER_WELL_LIST))
+
+MASTER_WELL_LIST[,1]
+FILL_IN_WELLS <- 	rep(t(MASTER_WELL_LIST[,1]),length(FILL_IN_YEARS))
+ZEROS <- rep(0,length(FILL_IN_WELLS))
+FILL_IN_YEARS <- rep( FILL_IN_YEARS,nrow(MASTER_WELL_LIST)) %>% sort 
+FILL_IN <- cbind(FILL_IN_WELLS,ZEROS,FILL_IN_YEARS) %>% as_tibble
+FILL_IN <- FILL_IN %>% mutate(ZEROS=as.numeric(ZEROS),FILL_IN_YEARS=as.numeric(FILL_IN_YEARS))
+
+colnames(FILL_IN) <- c("wdid",paste0("close_",colnames(DATA_ORIG)[PROGRAM_NUMBER]),"year")
+
+RES <- rbind(FILL_IN,RES) %>% unique
+	}
+	RES <- RES[,c(1,3,2)]
+	return(RES)
+}
+
+DF <- DF %>% left_join(GET_PROGRAM(4,DF)) %>%left_join(GET_PROGRAM(5,DF)) %>% left_join(GET_PROGRAM(6,DF) )
+CLOSE_YEAR <- 
+CLOSE_YEAR <-DF %>% filter(close_CREP_any==1) %>% group_by(wdid) %>% summarize(close_CREP_year=min(year)) %>% unique
+DF <- DF %>% left_join(CLOSE_YEAR) %>% mutate(close_CREP_year=ifelse(is.na(close_CREP_year),Inf,close_CREP_year))
+
+
+summary(DF$close_CREP_year)
+coefplot(feols(AF~sunab(close_CREP_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD1^year+SBD2^year+SBD3^year+SBD4^year+SBD5^year+SBD6^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year,data=DF ))
+
+etable(feols(AF~sunab(close_CREP_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD1^year+SBD2^year+SBD3^year+SBD4^year+SBD5^year+SBD6^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year,data=DF ),agg="cohort")
+
+