Major cleanup of scripts

1B_Run_Full_Simulation.r

@@ -4,7 +4,7 @@ library(fixest) #Estimating a model of birth rates, to provide variance in the b
 library(forecast) #Fore ARIMA migration simulations
 library(parallel)
 library(uuid) #To add a index to each batch
-####If the prelimnary data needs to be reloaded run the supplied bash script to download, process, and generate all needed data sets for the Monte Carlo population Simulation. Otherwise skip this step to save time
+####If the preliminary data needs to be reloaded run the supplied bash script to download, process, and generate all needed data sets for the Monte Carlo population Simulation. Otherwise skip this step to save time
 RELOAD_DATA <- FALSE 
 if(RELOAD_DATA){system("bash Prelim_Process.sh")}
 
@@ -115,16 +115,25 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL
 ###############NOTE NEED TO USE THIS AT END TO ADJUST THE RESULTS WHILE LEAVING THE DEMOGRAPHIC MATRIX
 	TERRA_POWER_EFFECT[1:7]	<- 	TERRA_POWER_EFFECT[1:7]+CONSTRUCTION_MIGRATION 
 	#Dry Pine Migration
-	DRY_PINE <- round(runif(1,0.15,0.45)*c(126,176,-302))
+	DRY_PINE <- round(runif(1,0.15,0.45)*c(126,176,-302)*1.15)
 	DRY_PINE[3] <- DRY_PINE[3]-sum(DRY_PINE)
 	OTHER_PROJECT_EFFECTS_TEMP[2:4] <- DRY_PINE
 	OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[2:4]  <-  cumsum(DRY_PINE)
 	#Data Center Simulation
 	if(runif(1)>0.67){
+	  MW_CAPACITY <- runif(1,20,120)
+	  PERM_WORKERS <- MW_CAPACITY/8+30 #Seems to be a good ration 32.5 workers at 20 capacity and 45 at max this is a heurstic I came up with
 	  OTHER_PROJECT_EFFECTS_PERM
 	  DATA_CENTER_START_YEAR <- which.max(runif(YEARS_AHEAD-6))+3
-	  PERM_DATA_CENTER_OP <- 40+2.515334+19.52286*ADDED_FROM_BASELINE
+	  PERM_DATA_CENTER_OP <- PERM_WORKERS*(1+(2.515334+19.52286*ADDED_FROM_BASELINE)/40) #Add induced
-	  TEMP_DATA_MIGRATION <- round(runif(1,0.5,1)*1.15*557*runif(1,0.41,1))#1.15 people per worker when including a family of 3.05 for 5% of construction workers. 41%-100% local migration. 557 from IMPLAN. Size of facility between 0.5 and 100% of estimated large capacity facility
+	 	 #Estimate construction workers on data center
+	  AVG_CONST_COST_PER_MW <- (10.4+13.2)/2
+	 RELATIVE_COST_DRAW <- runif(1,10.4*0.75,13.2*1.25)/AVG_CONST_COST_PER_MW #Draw from the range of reported construction costs per MW of capacity in Cheyene. See "DATA CENTER DEVELOPMENT COST GUIDE"
+	 PERCENT_IN_KEMMERER <- runif(1,0.41,1) #What percent of the in migration will live in Kemmerer vs other areas. Taken from TerraPower assumptions
+	 POP_PER_YEAR <- (1.15*46.02*MW_CAPACITY)/3 #MW capacity ranges from 18 to 120 based on historic Sabre projects. IMPLAN estimates 46.02 total jobs, per MW of capacity (using average cost per MWh). 1.15 accounts for 5% of workers bringing families
+	 TEMP_DATA_MIGRATION <- round(PERCENT_IN_KEMMERER*RELATIVE_COST_DRAW*(POP_PER_YEAR))
+
+
 	  OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+3] <- OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+3]+TEMP_DATA_MIGRATION #Enter
 	  OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+6] <- OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+6]-TEMP_DATA_MIGRATION #Leave
 	  OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]<-  OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]+TEMP_DATA_MIGRATION

1C_Run_Simulation_Upper_Bound.r

@@ -4,7 +4,7 @@ library(fixest) #Estimating a model of birth rates, to provide variance in the b
 library(forecast) #Fore ARIMA migration simulations
 library(parallel)
 library(uuid) #To add a index to each batch
-####If the prelimnary data needs to be reloaded run the supplied bash script to download, process, and generate all needed data sets for the Monte Carlo population Simulation. Otherwise skip this step to save time
+####If the preliminary data needs to be reloaded run the supplied bash script to download, process, and generate all needed data sets for the Monte Carlo population Simulation. Otherwise skip this step to save time
 RELOAD_DATA <- FALSE 
 if(RELOAD_DATA){system("bash Prelim_Process.sh")}
 
@@ -102,15 +102,44 @@ OPERATOR_MIGRATION <- OPERATOR_LIN_MIGRATION
 SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL,OPERATOR_MIGRATION,CONSTRUCTION_MIGRATION,MIGRATION_MULTIPLIERS ){
 	
 	TERRA_POWER_EFFECT <-	rep(0,YEARS_AHEAD)
+	OTHER_PROJECT_EFFECTS_PERM <-	rep(0,YEARS_AHEAD)
+	OTHER_PROJECT_EFFECTS_TEMP <-	rep(0,YEARS_AHEAD)
+	OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED <-	rep(0,YEARS_AHEAD)
 	OPERATOR_MIGRATION <- LOCAL_WORK_ADJ(OPERATOR_MIGRATION ,0.85) #Assume between 85%-100% operators live in Kemmerer
 	CONSTRUCTION_MIGRATION <- LOCAL_WORK_ADJ(CONSTRUCTION_MIGRATION,0.41) #Assume between 41%-100% operators live in Kemmerer
 	CONSTRUCTION_MIGRATION[7] <- CONSTRUCTION_MIGRATION[7] -       sum(CONSTRUCTION_MIGRATION )
 	CONSTRUCTION_POPULATION_ADDED <- cumsum(CONSTRUCTION_MIGRATION)
-	PERMANENT_TERRAPOWER_MIGRATION <- INDUCED_SIMULATION(CONSTRUCTION_MIGRATION,OPERATOR_MIGRATION,MIGRATION_MULTIPLIERS)+OPERATOR_MIGRATION
+	ADDED_FROM_BASELINE <- runif(1) #The percentage of the possible growth in industries like restaurants to add compared to the Kemmerer IMPLAN model which understates possible structural growth 
-		
+		PERMANENT_TERRAPOWER_MIGRATION <- INDUCED_SIMULATION(CONSTRUCTION_MIGRATION,OPERATOR_MIGRATION,MIGRATION_MULTIPLIERS,ADDED_FROM_BASELINE)+OPERATOR_MIGRATION
+
 ###############NOTE NEED TO USE THIS AT END TO ADJUST THE RESULTS WHILE LEAVING THE DEMOGRAPHIC MATRIX
 	TERRA_POWER_EFFECT[1:7]	<- 	TERRA_POWER_EFFECT[1:7]+CONSTRUCTION_MIGRATION 
-	MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL))+TERRA_POWER_EFFECT)
+	#Dry Pine Migration
+	DRY_PINE <- round(runif(1,0.15,0.45)*c(126,176,-302)*1.15)
+	DRY_PINE[3] <- DRY_PINE[3]-sum(DRY_PINE)
+	OTHER_PROJECT_EFFECTS_TEMP[2:4] <- DRY_PINE
+	OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[2:4]  <-  cumsum(DRY_PINE)
+	#Data Center Simulation
+	if(runif(1)>0.67){
+	  MW_CAPACITY <- runif(1,20,120)
+	  PERM_WORKERS <- MW_CAPACITY/8+30 #Seems to be a good ration 32.5 workers at 20 capacity and 45 at max this is a heurstic I came up with
+	  OTHER_PROJECT_EFFECTS_PERM
+	  DATA_CENTER_START_YEAR <- which.max(runif(YEARS_AHEAD-6))+3
+	  PERM_DATA_CENTER_OP <- PERM_WORKERS*(1+(2.515334+19.52286*ADDED_FROM_BASELINE)/40) #Add induced
+	 	 #Estimate construction workers on data center
+	  AVG_CONST_COST_PER_MW <- (10.4+13.2)/2
+	 RELATIVE_COST_DRAW <- runif(1,10.4*0.75,13.2*1.25)/AVG_CONST_COST_PER_MW #Draw from the range of reported construction costs per MW of capacity in Cheyene. See "DATA CENTER DEVELOPMENT COST GUIDE"
+	 PERCENT_IN_KEMMERER <- runif(1,0.41,1) #What percent of the in migration will live in Kemmerer vs other areas. Taken from TerraPower assumptions
+	 POP_PER_YEAR <- (1.15*46.02*MW_CAPACITY)/3 #MW capacity ranges from 18 to 120 based on historic Sabre projects. IMPLAN estimates 46.02 total jobs, per MW of capacity (using average cost per MWh). 1.15 accounts for 5% of workers bringing families
+	 TEMP_DATA_MIGRATION <- round(PERCENT_IN_KEMMERER*RELATIVE_COST_DRAW*(POP_PER_YEAR))
+
+
+	  OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+3] <- OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+3]+TEMP_DATA_MIGRATION #Enter
+	  OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+6] <- OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+6]-TEMP_DATA_MIGRATION #Leave
+	  OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]<-  OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]+TEMP_DATA_MIGRATION
+	  OTHER_PROJECT_EFFECTS_PERM[DATA_CENTER_START_YEAR+6] <- OTHER_PROJECT_EFFECTS_PERM[DATA_CENTER_START_YEAR+3]+PERM_DATA_CENTER_OP
+	}
+	MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL))+TERRA_POWER_EFFECT+OTHER_PROJECT_EFFECTS_PERM)
 	#The runif applies a downshift ranging from the historic decline rate all the way to the Lincoln rate applied in the model
 	
 	FINAL_REPORT_VALUES <- matrix(NA,ncol=6,nrow=YEARS_AHEAD)
@@ -125,6 +154,9 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL
 	}
 	FINAL_REPORT_VALUES[1:7,3] <- as.numeric(FINAL_REPORT_VALUES[1:7,3])+CONSTRUCTION_POPULATION_ADDED 
 	FINAL_REPORT_VALUES[1:7,6] <- as.numeric(FINAL_REPORT_VALUES[1:7,6]) +CONSTRUCTION_MIGRATION
+	
+	FINAL_REPORT_VALUES[,6] <- as.numeric(FINAL_REPORT_VALUES[,6]) +OTHER_PROJECT_EFFECTS_TEMP[1:nrow(FINAL_REPORT_VALUES)]
+	FINAL_REPORT_VALUES[,3] <- as.numeric(FINAL_REPORT_VALUES[,3]) +OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[1:nrow(FINAL_REPORT_VALUES)]
 	return(FINAL_REPORT_VALUES)
 }
 
@@ -135,10 +167,10 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL
 
 	NCORES <- detectCores()-1
 	BATCH_SIZE <- NCORES*10
-	TOTAL_SIMULATIONS <- 10^5
+	TOTAL_SIMULATIONS <- 10^6
 	N_RUNS <-ceiling(TOTAL_SIMULATIONS/BATCH_SIZE )
 
-	SIM_RES_FILE <- paste0(RES_SIM_DIR,"Kemmerer_2024_Simulation_County_Migration_Rate.csv")
+	SIM_RES_FILE <- paste0(RES_SIM_DIR,"Kemmerer_2024_Upper_Bound_With_Data_Center_Simulation.csv")
 NEW_RES_FILE <- !file.exists(SIM_RES_FILE)
 OPERATOR_LIN_MIGRATION <- OPERATORS %>% pull("Operator_Emp_Migrated")
 CONSTRUCTION_LIN_MIGRATION <- CONSTRUCTION %>% pull("Construction_Emp_Migrated")

1D_Run_Simulation_Lower_Bound.r

@@ -4,7 +4,7 @@ library(fixest) #Estimating a model of birth rates, to provide variance in the b
 library(forecast) #Fore ARIMA migration simulations
 library(parallel)
 library(uuid) #To add a index to each batch
-####If the prelimnary data needs to be reloaded run the supplied bash script to download, process, and generate all needed data sets for the Monte Carlo population Simulation. Otherwise skip this step to save time
+####If the preliminary data needs to be reloaded run the supplied bash script to download, process, and generate all needed data sets for the Monte Carlo population Simulation. Otherwise skip this step to save time
 RELOAD_DATA <- FALSE 
 if(RELOAD_DATA){system("bash Prelim_Process.sh")}
 
@@ -102,15 +102,44 @@ OPERATOR_MIGRATION <- OPERATOR_LIN_MIGRATION
 SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL,OPERATOR_MIGRATION,CONSTRUCTION_MIGRATION,MIGRATION_MULTIPLIERS ){
 	
 	TERRA_POWER_EFFECT <-	rep(0,YEARS_AHEAD)
+	OTHER_PROJECT_EFFECTS_PERM <-	rep(0,YEARS_AHEAD)
+	OTHER_PROJECT_EFFECTS_TEMP <-	rep(0,YEARS_AHEAD)
+	OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED <-	rep(0,YEARS_AHEAD)
 	OPERATOR_MIGRATION <- LOCAL_WORK_ADJ(OPERATOR_MIGRATION ,0.85) #Assume between 85%-100% operators live in Kemmerer
 	CONSTRUCTION_MIGRATION <- LOCAL_WORK_ADJ(CONSTRUCTION_MIGRATION,0.41) #Assume between 41%-100% operators live in Kemmerer
 	CONSTRUCTION_MIGRATION[7] <- CONSTRUCTION_MIGRATION[7] -       sum(CONSTRUCTION_MIGRATION )
 	CONSTRUCTION_POPULATION_ADDED <- cumsum(CONSTRUCTION_MIGRATION)
-	PERMANENT_TERRAPOWER_MIGRATION <- INDUCED_SIMULATION(CONSTRUCTION_MIGRATION,OPERATOR_MIGRATION,MIGRATION_MULTIPLIERS)+OPERATOR_MIGRATION
+	ADDED_FROM_BASELINE <- runif(1) #The percentage of the possible growth in industries like restaurants to add compared to the Kemmerer IMPLAN model which understates possible structural growth 
-		
+		PERMANENT_TERRAPOWER_MIGRATION <- INDUCED_SIMULATION(CONSTRUCTION_MIGRATION,OPERATOR_MIGRATION,MIGRATION_MULTIPLIERS,ADDED_FROM_BASELINE)+OPERATOR_MIGRATION
+
 ###############NOTE NEED TO USE THIS AT END TO ADJUST THE RESULTS WHILE LEAVING THE DEMOGRAPHIC MATRIX
 	TERRA_POWER_EFFECT[1:7]	<- 	TERRA_POWER_EFFECT[1:7]+CONSTRUCTION_MIGRATION 
-	MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL))-55+TERRA_POWER_EFFECT)
+	#Dry Pine Migration
+	DRY_PINE <- round(runif(1,0.15,0.45)*c(126,176,-302)*1.15)
+	DRY_PINE[3] <- DRY_PINE[3]-sum(DRY_PINE)
+	OTHER_PROJECT_EFFECTS_TEMP[2:4] <- DRY_PINE
+	OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[2:4]  <-  cumsum(DRY_PINE)
+	#Data Center Simulation
+	if(runif(1)>0.67){
+	  MW_CAPACITY <- runif(1,20,120)
+	  PERM_WORKERS <- MW_CAPACITY/8+30 #Seems to be a good ration 32.5 workers at 20 capacity and 45 at max this is a heurstic I came up with
+	  OTHER_PROJECT_EFFECTS_PERM
+	  DATA_CENTER_START_YEAR <- which.max(runif(YEARS_AHEAD-6))+3
+	  PERM_DATA_CENTER_OP <- PERM_WORKERS*(1+(2.515334+19.52286*ADDED_FROM_BASELINE)/40) #Add induced
+	 	 #Estimate construction workers on data center
+	  AVG_CONST_COST_PER_MW <- (10.4+13.2)/2
+	 RELATIVE_COST_DRAW <- runif(1,10.4*0.75,13.2*1.25)/AVG_CONST_COST_PER_MW #Draw from the range of reported construction costs per MW of capacity in Cheyene. See "DATA CENTER DEVELOPMENT COST GUIDE"
+	 PERCENT_IN_KEMMERER <- runif(1,0.41,1) #What percent of the in migration will live in Kemmerer vs other areas. Taken from TerraPower assumptions
+	 POP_PER_YEAR <- (1.15*46.02*MW_CAPACITY)/3 #MW capacity ranges from 18 to 120 based on historic Sabre projects. IMPLAN estimates 46.02 total jobs, per MW of capacity (using average cost per MWh). 1.15 accounts for 5% of workers bringing families
+	 TEMP_DATA_MIGRATION <- round(PERCENT_IN_KEMMERER*RELATIVE_COST_DRAW*(POP_PER_YEAR))
+
+
+	  OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+3] <- OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+3]+TEMP_DATA_MIGRATION #Enter
+	  OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+6] <- OTHER_PROJECT_EFFECTS_TEMP[DATA_CENTER_START_YEAR+6]-TEMP_DATA_MIGRATION #Leave
+	  OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]<-  OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]+TEMP_DATA_MIGRATION
+	  OTHER_PROJECT_EFFECTS_PERM[DATA_CENTER_START_YEAR+6] <- OTHER_PROJECT_EFFECTS_PERM[DATA_CENTER_START_YEAR+3]+PERM_DATA_CENTER_OP
+	}
+	MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL))-55+TERRA_POWER_EFFECT+OTHER_PROJECT_EFFECTS_PERM)
 	#The runif applies a downshift ranging from the historic decline rate all the way to the Lincoln rate applied in the model
 	
 	FINAL_REPORT_VALUES <- matrix(NA,ncol=6,nrow=YEARS_AHEAD)
@@ -125,6 +154,9 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL
 	}
 	FINAL_REPORT_VALUES[1:7,3] <- as.numeric(FINAL_REPORT_VALUES[1:7,3])+CONSTRUCTION_POPULATION_ADDED 
 	FINAL_REPORT_VALUES[1:7,6] <- as.numeric(FINAL_REPORT_VALUES[1:7,6]) +CONSTRUCTION_MIGRATION
+	
+	FINAL_REPORT_VALUES[,6] <- as.numeric(FINAL_REPORT_VALUES[,6]) +OTHER_PROJECT_EFFECTS_TEMP[1:nrow(FINAL_REPORT_VALUES)]
+	FINAL_REPORT_VALUES[,3] <- as.numeric(FINAL_REPORT_VALUES[,3]) +OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[1:nrow(FINAL_REPORT_VALUES)]
 	return(FINAL_REPORT_VALUES)
 }
 
@@ -135,10 +167,10 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL
 
 	NCORES <- detectCores()-1
 	BATCH_SIZE <- NCORES*10
-	TOTAL_SIMULATIONS <- 10^5
+	TOTAL_SIMULATIONS <- 10^6
 	N_RUNS <-ceiling(TOTAL_SIMULATIONS/BATCH_SIZE )
 
-	SIM_RES_FILE <- paste0(RES_SIM_DIR,"Kemmerer_2024_Simulation_Historic_Migration_Rate.csv")
+	SIM_RES_FILE <- paste0(RES_SIM_DIR,"Kemmerer_2024_With_Lower_Bound_With_Data_Center_Simulation.csv")
 NEW_RES_FILE <- !file.exists(SIM_RES_FILE)
 OPERATOR_LIN_MIGRATION <- OPERATORS %>% pull("Operator_Emp_Migrated")
 CONSTRUCTION_LIN_MIGRATION <- CONSTRUCTION %>% pull("Construction_Emp_Migrated")

2A_2016_Result_Analysis_2016_and_Historic_Trend.r

@@ -0,0 +1,68 @@
+library(tidyverse)
+	if(!exists("SAVE_PRELIM_LOC")){SAVE_PRELIM_LOC <- "./Results/Primary_Simulation_Results/Preliminary_Results/"}
+	dir.create(SAVE_PRELIM_LOC, recursive = TRUE, showWarnings = FALSE)
+
+	RES <- read_csv("Results/Simulations/Kemmerer_2016_Simulation.csv")
+source("Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r")
+	HIST <-	readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>=1940)
+	POP_DATA <- GET_DATA(RES,3,2017)
+	POP_PLOT <- MAKE_GRAPH(POP_DATA,COLOR="orchid3")
+	POP_PLOT <- POP_PLOT+geom_line(data=HIST,aes(x=Year,y=Population),color='orange',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2060, by = 10),2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 500))+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Population")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
+
+	POP_PLOT_ZOOM <- MAKE_GRAPH(POP_DATA,COLOR="orchid3")
+	HIST_ZOOM <-HIST %>% filter(Year>=2010)
+	POP_PLOT_ZOOM <- POP_PLOT_ZOOM+geom_line(data=HIST_ZOOM,aes(x=Year,y=Population),color='orange1',linewidth=0.75)+ scale_x_continuous(breaks = seq(2000, 2025, by = 1))+ scale_y_continuous(breaks = seq(0, 35000, by = 100))+ggtitle("Kemmerer & Diamondville, Population Forecast (2016)")+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Population")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
+png(paste0(SAVE_PRELIM_LOC,"Population_Fan_Chart_2016_Results.png"), width = 12, height = 8, units = "in", res = 600)
+POP_PLOT 
+dev.off()
+
+png(paste0(SAVE_PRELIM_LOC,"Population_Fan_Chart_2016_Results_Zoomed_In.png"), width = 6, height = 6, units = "in", res = 600)
+POP_PLOT_ZOOM
+dev.off()
+#####################
+BIRTH_DATA <- GET_DATA(RES,4,2017)
+BIRTH_PLOT <- MAKE_GRAPH(BIRTH_DATA,COLOR="orchid3")
+
+BIRTH_PLOT <- BIRTH_PLOT+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Births),color='orange1',linewidth=0.75)+ scale_x_continuous(breaks = seq(2000, 2025, by = 1))+ scale_y_continuous(breaks = seq(0, 100, by = 5))+ expand_limits( y = 0)+ggtitle("Kemmerer & Diamondville, Birth Forecast (2016)")+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Births")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
+png(paste0(SAVE_PRELIM_LOC,"Birth_Fan_Chart_2016.png"), width = 12, height = 8, units = "in", res = 600)
+	BIRTH_PLOT 
+dev.off()
+#################
+DEATH_DATA <- GET_DATA(RES,5,2017) %>% filter(!is.na(MIN))
+DEATH_PLOT <- MAKE_GRAPH(DEATH_DATA,COLOR="orchid3")+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Deaths),color='orange1',linewidth=0.75)+ scale_x_continuous(breaks = seq(2000, 2025, by = 1))+ scale_y_continuous(breaks = seq(0, 35000, by = 5))+ expand_limits( y = 0)+ggtitle("Kemmerer & Diamondville, Mortality Forecast (2016)")+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Deaths")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
+
+png(paste0(SAVE_PRELIM_LOC,"Mortality_Fan_Chart_2016.png"), width = 12, height = 8, units = "in", res = 600)
+	DEATH_PLOT
+dev.off()
+
+#################
+MIGRATION_DATA <- GET_DATA(RES,6,2017) %>% filter(!is.na(MIN))
+MIGRATION_PLOT <- MAKE_GRAPH(MIGRATION_DATA,COLOR="orchid3")+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Migration),color='orange1',linewidth=0.75)+ scale_x_continuous(breaks = seq(2000, 2025, by = 1))+ scale_y_continuous(breaks = seq(-1000, 1000, by = 25))+ expand_limits( y = 0)+ggtitle("Kemmerer & Diamondville, Migration Forecast (2016)")+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Migration")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
+
+png(paste0(SAVE_PRELIM_LOC,"Migration_Fan_Chart_2016.png"), width = 12, height = 8, units = "in", res = 600)
+	MIGRATION_PLOT	
+dev.off()
+
+######################
+OTHER <-	readRDS("Data/Cleaned_Data/Population_Data/RDS/Full_Lincoln_County_Population_Data.Rds") %>% mutate(Region='Rest of Lincoln County') %>% select(Year,Population,Region)
+
+KEM <-	readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% select(Year,Population,Region) %>% select(Year,Population,Region)
+TEMP <- KEM %>% rename(KEM_POP=Population) %>% select(Year,KEM_POP)
+OTHER <- OTHER %>% left_join(TEMP) %>% filter(!is.na(KEM_POP)) %>% mutate(Population=Population-KEM_POP) %>% select(Year,Population,Region)
+
+PLOT_DATA <- rbind(OTHER,KEM)
+#%>% group_by(Year) %>% mutate(Population=ifelse(Region=='Rest of Lincoln County',Population-min(Population),Population)) %>% ungroup
+
+#PLOT_DATA <- rbind(PLOT_DATA[,c(1:3,4)] %>% rename(Value=Population) %>% mutate(Category="Population"),
+#PLOT_DATA[,c(1:3,5)] %>% rename(Value=Births) %>% mutate(Category="Births"),
+#PLOT_DATA[,c(1:3,6)] %>% rename(Value=Deaths) %>% mutate(Category="Mortality"),
+#PLOT_DATA[,c(1:3,7)] %>% rename(Value=Migration) %>% mutate(Category="Net Migration"))
+
+
+
+POP_DIFF_PLOT <- ggplot(PLOT_DATA ,aes(x=Year,y=Population,group=Region,color=Region))+geom_line(linewidth=2)+scale_color_manual(values=c("cadetblue","darkred"))+ scale_x_continuous(breaks = c(seq(1910, 2025, by = 10)))+ scale_y_continuous(breaks = seq(0, 35000, by = 1000))+ theme_bw()+ theme(legend.position = "bottom")
+png(paste0(SAVE_PRELIM_LOC,"Population_Historic_Plot.png"), width = 12, height = 8, units = "in", res = 600)
+POP_DIFF_PLOT 
+dev.off()
+
+

2A_Result_Analysis_2016.r

@@ -1,65 +0,0 @@
-library(tidyverse)
-###Process the simulations and save the main percentile results by year
-	RES <- read_csv("Results/Simulations/Kemmerer_2016_Simulation.csv")
-	HIST <-	readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>=1940)
-	YEARS <- min(RES$Year):max(RES$Year)
-	######Population
-	START_POP <- HIST %>% filter(Year==2016) %>% pull(Population)
-	GRAPH_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(RES %>% filter(Year==x) %>% pull(Population),c(0.025,0.05,0.1,0.25,0.4,0.5,0.6,0.75,0.9,0.95,0.975))}))  %>% as_tibble %>% mutate(Year=YEARS)
-	#Add 2016 as a starting value
-	GRAPH_DATA <- 	rbind(GRAPH_DATA[1,],GRAPH_DATA)
-	GRAPH_DATA[1,]  <-t(c(rep(START_POP,ncol(GRAPH_DATA)-1),min(GRAPH_DATA$Year)-1))
-	FAN_DATA <- GRAPH_DATA
-	GRAPH_DATA <- GRAPH_DATA %>% pivot_longer(cols=!Year,names_to=c("Percentile"),values_to="Population")
-	GRAPH_DATA$Percentile <-	factor(GRAPH_DATA$Percentile,levels=rev(c('2.5%','5%','10%','25%','40%','50%','60%','75%','90%','95%','97.5%')))
-	START_POP <- 	HIST %>% filter(Year==2016) %>% pull(Population)
-	GRAPH_DATA <-	rbind(GRAPH_DATA,GRAPH_DATA  %>% filter(Year==2017) %>% mutate(Population=START_POP,Year=2016))
-	MEDIAN_PRED <- 	GRAPH_DATA %>% filter(Percentile=='50%')
-	######Migration
-	GRAPH_DATA_MIGRATION <- do.call(rbind,lapply(YEARS,function(x){quantile(RES %>% filter(Year==x) %>% pull(Net_Migration),c(0.025,0.05,0.1,0.25,0.4,0.5,0.6,0.75,0.9,0.95,0.975))}))  %>% as_tibble %>% mutate(Year=YEARS)
-	FAN_DATA_MIGRATION <- GRAPH_DATA_MIGRATION
-	GRAPH_DATA_MIGRATION <- GRAPH_DATA_MIGRATION %>% pivot_longer(cols=!Year,names_to=c("Percentile"),values_to="Migration")
-	GRAPH_DATA_MIGRATION$Percentile <-	factor(GRAPH_DATA_MIGRATION$Percentile,levels=rev(c('2.5%','5%','10%','25%','40%','60%','75%','90%','95%','97.5%')))
-	MEDIAN_PRED_MIGRATION <- GRAPH_DATA_MIGRATION %>% filter(Percentile=='50%')
-	######Mortalities
-	GRAPH_DATA_MORTALITY <- do.call(rbind,lapply(YEARS,function(x){quantile(RES %>% filter(Year==x) %>% pull(Deaths),c(0.025,0.05,0.1,0.25,0.4,0.5,0.6,0.75,0.9,0.95,0.975))}))  %>% as_tibble %>% mutate(Year=YEARS)
-	FAN_DATA_MORTALITY <- GRAPH_DATA_MORTALITY
-	GRAPH_DATA_MORTALITY <- GRAPH_DATA_MORTALITY %>% pivot_longer(cols=!Year,names_to=c("Percentile"),values_to="Deaths")
-	GRAPH_DATA_MORTALITY$Percentile <-	factor(GRAPH_DATA_MORTALITY$Percentile,levels=rev(c('2.5%','5%','10%','25%','40%','60%','75%','90%','95%','97.5%')))
-	MEDIAN_PRED_MORTALITY<- GRAPH_DATA_MORTALITY %>% filter(Percentile=='50%')
-	######Births
-	START_BIRTHS <- HIST %>% filter(Year==2016) %>% pull(Births)
-
-	GRAPH_DATA_BIRTHS <- do.call(rbind,lapply(YEARS,function(x){quantile(RES %>% filter(Year==x) %>% pull(Births),c(0.025,0.05,0.1,0.25,0.4,0.5,0.6,0.75,0.9,0.95,0.975))}))  %>% as_tibble %>% mutate(Year=YEARS)
-	GRAPH_DATA_BIRTHS <- 	rbind(GRAPH_DATA_BIRTHS[1,],GRAPH_DATA_BIRTHS)
-	GRAPH_DATA_BIRTHS[1,]  <-t(c(rep(START_BIRTHS,ncol(GRAPH_DATA_BIRTHS)-1),min(GRAPH_DATA_BIRTHS$Year)-1))
-
-	FAN_DATA_BIRTHS <- GRAPH_DATA_BIRTHS
-	GRAPH_DATA_BIRTHS<- GRAPH_DATA_BIRTHS %>% pivot_longer(cols=!Year,names_to=c("Percentile"),values_to="Births")
-	GRAPH_DATA_BIRTHS$Percentile <-	factor(GRAPH_DATA_BIRTHS$Percentile,levels=rev(c('2.5%','5%','10%','25%','40%','60%','75%','90%','95%','97.5%')))
-	MEDIAN_PRED_BIRTHS<- GRAPH_DATA_BIRTHS %>% filter(Percentile=='50%')
-
-
-
-	#write_csv(GRAPH_DATA,PERCENTILE_DATA)
-	#Add historic
-	MEDIAN_PRED <- 	GRAPH_DATA %>% filter(Percentile=='50%')
-	GRAPH_DATA <- 	GRAPH_DATA %>% filter(Percentile!='50%')
-
-	ALPHA=0.2
-	COLOR <- 'black'
-#GRAPH <-
-	nrow(RES)/10
- ggplot(data=GRAPH_DATA)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Population,group=Percentile,color=Percentile))+geom_line(data=HIST,aes(x=Year,y=Population),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED,aes(x=Year,y=Population),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2030, by = 5)))+ scale_y_continuous(breaks = seq(0, 35000, by = 250))+theme_bw()+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0)
-
- ggplot(data=GRAPH_DATA_MIGRATION)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MIGRATION,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Migration,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Migration),color='black',linewidth=0.75)
- +geom_line(data=MEDIAN_PRED_MIGRATION,aes(x=Year,y=Migration),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2030, by = 5)))+ scale_y_continuous(breaks = seq(0, 35000, by = 250))+theme_bw()+ggtitle("Kemmerer & Diamondville, Migration Forecast")+ expand_limits( y = 0)
-
-
-
-
- ggplot(data=GRAPH_DATA_MORTALITY)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_MORTALITY,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Deaths,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Deaths),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED_MORTALITY,aes(x=Year,y=Deaths),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2000, 2030, by = 5)))+ scale_y_continuous(breaks = seq(0, 50, by = 5))+theme_bw()+ggtitle("Kemmerer, Wyoming Death Forecast")+ expand_limits( y = 0)
-
-
-ggplot(data=GRAPH_DATA_BIRTHS)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA_BIRTHS,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Births,group=Percentile,color=Percentile))+geom_line(data=HIST %>% filter(Year>=2009),aes(x=Year,y=Births),color='black',linewidth=0.75)+geom_line(data=MEDIAN_PRED_BIRTHS,aes(x=Year,y=Births),color='black',linetype=4,linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2000, 2030, by = 5)))+ scale_y_continuous(breaks = seq(0, 200, by = 5))+theme_bw()+ggtitle("Kemmerer, Wyoming Birth Forecast")+ expand_limits( y = 0)
-

2B_Result_Analysis.r

@@ -20,7 +20,7 @@ POP_PLOT
 dev.off()
 
 #Population Summary table for the key years
-KEY_YEARS <- c(2029,2030,2035,2045,2055,2065)
+KEY_YEARS <- c(2027,2030,2035,2045,2055,2065)
 KEY_YEAR_SUMMARY <- POP_DATA %>% filter(Year %in% KEY_YEARS,Group==5) %>% select(-Group)  %>% rename(CI_95_Lower=MIN,CI_95_Upper=MAX) %>% left_join(POP_DATA %>% filter(Year %in% KEY_YEARS,Group==0) %>% select(-Group)  %>% select(Year,Median=MIN)) %>% select(Year,CI_95_Lower,Median,CI_95_Upper)
 KEY_YEAR_SUMMARY <- round(KEY_YEAR_SUMMARY )
 KEY_YEAR_SUMMARY <- t(KEY_YEAR_SUMMARY )
@@ -49,21 +49,23 @@ png(paste0(SAVE_RES_LOC,"Migration_Fan_Chart_Main_Results.png"), width = 12, hei
 dev.off()
 
 #####Key year table
-KEY <- RES %>% filter(Year %in% c(2029,2030,2035,2045,2055,2065))
+
- AVG_VALUES <- KEY %>% group_by(Year) %>% summarize(MED=median(Population),MEAN=mean(Population))
+YEARS <- c(2027,2030,2035,2045,2055,2065)
+KEY_YEAR_DATA <- RES %>% filter(Year %in% YEARS)
+ AVG_VALUES <- KEY_YEAR_DATA %>% group_by(Year) %>% summarize(MED=median(Population),MEAN=mean(Population))
 AVG_VALUES <- rbind(AVG_VALUES[,1:2]%>% rename(Value=MED) %>% mutate('Summary Stat.'="Median"),AVG_VALUES[,c(1,3)] %>% rename(Value=MEAN) %>% mutate('Summary Stat.'="Mean"))
 
-HISTOGRAM <- ggplot(KEY, aes(x = Population,group=-Year,Color=Year,fill=Year)) + geom_histogram(alpha=0.3,bins=800)+geom_vline(data = AVG_VALUES, aes(xintercept = Value,group=`Summary Stat.`,color = `Summary Stat.`), size = 0.75)+scale_color_manual(values=c("red","black","black"))+ facet_grid(rows=vars(Year))+ scale_x_continuous(breaks = c(seq(0, 10000, by = 500)))+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())+ylab("Number of Simulation")+guides(fill= guide_legend(nrow = 1))
+HISTOGRAM <- ggplot(KEY_YEAR_DATA, aes(x = Population,group=-Year,Color=Year,fill=Year)) + geom_histogram(alpha=0.3,bins=800)+geom_vline(data = AVG_VALUES, aes(xintercept = Value,group=`Summary Stat.`,color = `Summary Stat.`), linewidth= 0.75)+scale_color_manual(values=c("red","black","black"))+ facet_grid(rows=vars(Year))+ scale_x_continuous(breaks = c(seq(0, 10000, by = 500)))+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())+ylab("Number of Simulation")+guides(fill= guide_legend(nrow = 1))
 png(paste0(SAVE_RES_LOC,"Population_Histogram_Main_Results.png"), width = 8, height = 12, units = "in", res = 600)
 	HISTOGRAM
 dev.off()
-
 #rm(KEY_YEARS)
 POP_LEVELS <- seq(2000,6000,100)
-YEARS <- c(2030,2035,2045,2055,2065)
+YEARS
+if(exists("KEY_YEARS")){rm(KEY_YEARS)}
 for(i in YEARS ){
-	KEY <- RES %>% filter(Year==i )  %>% pull(Population)
+	CURRENT_YEAR_DATA <- RES %>% filter(Year==i )  %>% pull(Population)
-	ECDF <- ecdf(KEY)
+	ECDF <- ecdf(CURRENT_YEAR_DATA)
 	ECDF_VALUES <- ECDF(POP_LEVELS)
 	if(!exists("KEY_YEARS")){KEY_YEARS<- ECDF_VALUES} else{KEY_YEARS<- cbind(KEY_YEARS,ECDF_VALUES)}
 
@@ -76,5 +78,5 @@ PLOT_RED <- "red"
 KEY_YEARS <- KEY_YEARS%>% as.data.frame
 Capacity_Risk <- KEY_YEARS%>%  gt(rownames_to_stub = TRUE,caption="Year")  %>%  data_color(    fn = scales::col_numeric( palette = c(PLOT_RED, PLOT_YELLOW, PLOT_GREEN),  domain = c(0, 1) ) ) %>%  fmt_percent( decimals = 1,  drop_trailing_zeros = FALSE) %>% tab_stubhead(label =c("Capacity"))
 gtsave( data = Capacity_Risk ,  filename = "./Results/Primary_Simulation_Results/Main_Results/Capacity_Table_Main_Results.html")
-system("wkhtmltopdf --disable-smart-shrinking --no-stop-slow-scripts --enable-local-file-access --page-width 85mm --page-height 328mm ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table_Main_Results.html ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table_Main_Results.pdf") 
+system("wkhtmltopdf --disable-smart-shrinking --no-stop-slow-scripts --enable-local-file-access --page-width 99mm --page-height 328mm ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table_Main_Results.html ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table_Main_Results.pdf") 
 

2C_Upper_Bound_Result_Analysis.r

@@ -6,7 +6,7 @@ source("Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r") #Function
 
 
 ###Process the simulations and save the main percentile results by year
-	RES <- read_csv("Results/Simulations/Kemmerer_2024_Simulation_County_Migration_Rate.csv")
+	RES <- read_csv("Results/Simulations/Kemmerer_2024_Upper_Bound_With_Data_Center_Simulation.csv")
 	RES[,"Year"] <- RES[,"Year"]
 	RES<-	RES %>% filter(!is.na(Year))
 	RES <-	RES %>% filter(!is.na(Population))
@@ -61,6 +61,7 @@ dev.off()
 #rm(KEY_YEARS)
 POP_LEVELS <- seq(2000,6000,100)
 YEARS <- c(2030,2035,2045,2055,2065)
+if(exists("KEY_YEARS")){rm(KEY_YEARS)}
 for(i in YEARS ){
 	KEY <- RES %>% filter(Year==i )  %>% pull(Population)
 	ECDF <- ecdf(KEY)

2D_Lower_Bound_Result_Analysis.r

@@ -6,7 +6,7 @@ source("Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r") #Function
 
 
 ###Process the simulations and save the main percentile results by year
-	RES <- read_csv("Results/Simulations/Kemmerer_2024_Simulation_Historic_Migration_Rate.csv")
+	RES <- read_csv("Results/Simulations/Kemmerer_2024_With_Lower_Bound_With_Data_Center_Simulation.csv")
 	RES[,"Year"] <- RES[,"Year"]
 	RES<-	RES %>% filter(!is.na(Year))
 	RES <-	RES %>% filter(!is.na(Population))
@@ -62,6 +62,7 @@ dev.off()
 #rm(KEY_YEARS)
 POP_LEVELS <- seq(2000,6000,100)
 YEARS <- c(2030,2035,2045,2055,2065)
+if(exists("KEY_YEARS")){rm(KEY_YEARS)}
 for(i in YEARS ){
 	KEY <- RES %>% filter(Year==i )  %>% pull(Population)
 	ECDF <- ecdf(KEY)

Create_Result_Figures.sh

@@ -1,4 +1,4 @@
-#Rscript 2A_Result_Analysis_2016.r
+Rscript 2A_2016_Result_Analysis_2016_and_Historic_Trend.r 
 Rscript 2B_Result_Analysis.r
 Rscript 2C_Upper_Bound_Result_Analysis.r
 Rscript 2D_Lower_Bound_Result_Analysis.r

Run_Simulations.sh

@@ -37,7 +37,7 @@ printf "Execution time: %02d hours, %02d minutes, %02d seconds\n" "$hours" "$min
 
 
 start_time=$(date +%s)
-Rscript 1B_Run_Full_Simulation.r
+Rscript 1B_Run_Full_Simulation.r 
 echo "Main results with variable migration rates complete! 1 million simulations total"
 end_time=$(date +%s)
 elapsed_seconds=$((end_time - start_time))

Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r

@@ -1,17 +1,18 @@
 #Functions to create consistent fan functions from simulations
 	###Create a data set which pairs of the upper and lower confidence intervals. This allows the layers to be plotted
-GET_DATA <- function(RES,COL_NUM){ 
+GET_DATA <- function(RES,COL_NUM,SIM_START=2025){ 
+	LAST_HISTORIC <- SIM_START-1
 	YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE)
 	LEVELS <- seq(0.01,0.99,by=0.005)
 	GROUPS <- floor(length(LEVELS)/2)
 	FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((RES %>% filter(Year==x))[,COL_NUM])),LEVELS)}))  %>% as_tibble %>% mutate(Year=YEARS)
 	FAN_DATA <- 	rbind(FAN_DATA[1,],FAN_DATA)
-	START_VALUE <- (HIST %>% filter(Year==2024))[,COL_NUM+1] %>% as.numeric
+	START_VALUE <- (HIST %>% filter(Year==LAST_HISTORIC ))[,COL_NUM+1] %>% as.numeric
-	FAN_DATA <- 	FAN_DATA  %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>2024) %>% unique
+	FAN_DATA <- 	FAN_DATA  %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>LAST_HISTORIC) %>% unique
 	NUM_YEARS <- length(unique(FAN_DATA$Year) )
 	FAN_DATA$Group <- rep(c(1:GROUPS,0,rev(1:GROUPS)),NUM_YEARS)
 	FAN_DATA <- FAN_DATA %>% group_by(Year,Group)  %>% summarize(MIN=min(value),MAX=max(value))
-	TEMP <- FAN_DATA %>% filter(Year==2025) %>% mutate(Year=2024) %>% ungroup
+	TEMP <- FAN_DATA %>% filter(Year==SIM_START) %>% mutate(Year=LAST_HISTORIC) %>% ungroup
 	TEMP[,3:4] <- START_VALUE 
 	FAN_DATA <- rbind(TEMP,FAN_DATA %>% ungroup) %>% as_tibble
 	return(FAN_DATA)