diff --git a/1B_Run_Full_Simulation.r b/1B_Run_Full_Simulation.r index 16aa46a..cce55e0 100644 --- a/1B_Run_Full_Simulation.r +++ b/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 - 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 + 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 diff --git a/1C_Run_Simulation_Upper_Bound.r b/1C_Run_Simulation_Upper_Bound.r index f4c1c1e..c7665e6 100644 --- a/1C_Run_Simulation_Upper_Bound.r +++ b/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") diff --git a/1D_Run_Simulation_Lower_Bound.r b/1D_Run_Simulation_Lower_Bound.r index 2b8a5d7..e00210d 100644 --- a/1D_Run_Simulation_Lower_Bound.r +++ b/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") diff --git a/2A_2016_Result_Analysis_2016_and_Historic_Trend.r b/2A_2016_Result_Analysis_2016_and_Historic_Trend.r new file mode 100644 index 0000000..ebe66f5 --- /dev/null +++ b/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() + + diff --git a/2A_Result_Analysis_2016.r b/2A_Result_Analysis_2016.r deleted file mode 100644 index fb6c507..0000000 --- a/2A_Result_Analysis_2016.r +++ /dev/null @@ -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) - diff --git a/2B_Result_Analysis.r b/2B_Result_Analysis.r index 6ffaef3..aaa6700 100644 --- a/2B_Result_Analysis.r +++ b/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) - ECDF <- ecdf(KEY) + CURRENT_YEAR_DATA <- RES %>% filter(Year==i ) %>% pull(Population) + 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") diff --git a/2C_Upper_Bound_Result_Analysis.r b/2C_Upper_Bound_Result_Analysis.r index 0ae073f..4e0e1fb 100644 --- a/2C_Upper_Bound_Result_Analysis.r +++ b/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) diff --git a/2D_Lower_Bound_Result_Analysis.r b/2D_Lower_Bound_Result_Analysis.r index c3145c5..76a7bb8 100644 --- a/2D_Lower_Bound_Result_Analysis.r +++ b/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) diff --git a/Create_Result_Figures.sh b/Create_Result_Figures.sh index 4447b38..4a33db4 100644 --- a/Create_Result_Figures.sh +++ b/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 diff --git a/Run_Simulations.sh b/Run_Simulations.sh index 118cc7d..ba2340e 100644 --- a/Run_Simulations.sh +++ b/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)) diff --git a/Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r b/Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r index 2fc1724..98cd2e7 100644 --- a/Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r +++ b/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 - FAN_DATA <- FAN_DATA %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>2024) %>% unique + 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>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)