From 034f69924b93737e1b5f4d029a1cec0cd4b39ea6 Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Thu, 23 Oct 2025 16:56:50 -0600 Subject: [PATCH] Made loop, working on data clean --- 1_Run_Full_Simulation.r | 31 +++++++----- Migration_Regression.r | 87 +++++++++++++++++++++++++++++++++ Scripts/Data_Load.r | 2 +- Scripts/Monte_Carlo_Functions.r | 3 +- 4 files changed, 109 insertions(+), 14 deletions(-) create mode 100644 Migration_Regression.r diff --git a/1_Run_Full_Simulation.r b/1_Run_Full_Simulation.r index ab4d373..0b96c0c 100644 --- a/1_Run_Full_Simulation.r +++ b/1_Run_Full_Simulation.r @@ -6,7 +6,7 @@ source("Scripts/Birth_Simulation_Functions.r") source("Scripts/Monte_Carlo_Functions.r") #####User Configuration Values - NUM_SIMULATIONS <- 10^6 #Number of Monte Carlo Simulations to run + NUM_SIMULATIONS <- 10^4 #Number of Monte Carlo Simulations to run RERUN_MORTALITY_SIMULATION <- TRUE #Rerun the Monte Carlo simulation of future mortality rates (not actual deaths) even if a Rds file of a mortality rates exists. This can be used to speed up runs when FALSE NUM_YEARS_PROJECTED <- 50 #How many years into the future should each Monte Carlo run project to. For example 25 years if starting from 2025 and ending in 2050. BIRTH_RATE_REG_RESULTS <- "Data/Regression_Results/Birth_Rate_Model.Rds" #Location of the regression used to model variance in birth rates in each Monte Carlo simulation. @@ -22,17 +22,24 @@ START_DEM_DATA <- readRDS("Data/Cleaned_Data/Lincoln_Demographic_Data.Rds") %>% MORTALITY_SIMULATION <- readRDS("./Data/Simulated_Data_Sets/MORTALITY_MONTE_CARLO.Rds") #Load the Mortality simulation to speed up simulation #First run -C_RES <- RUN_SINGLE_SIM(MOD_BIRTHS,FIRST_PREDICT_YEAR_POPULATION_DATA,START_DEM_DATA,MORTALITY_SIMULATION,SIM_NUMBER=1,START_OF_SIM=2023) -C_RES[[3]] +C_RES <-RUN_SINGLE_SIM(MOD_BIRTHS,FIRST_PREDICT_YEAR_POPULATION_DATA,START_DEM_DATA,MORTALITY_SIMULATION,SIM_NUMBER=1,START_OF_SIM=2023) +# Year County Population Births Deaths Migration Min_Birth_Group PREV_BIRTH PREV_TWO_BIRTH Male_Birth_Group #Second run, work on making into loop and saving the output to file (check if that will slow the results). Maybe use a predifined matrix, so that the results can be stored quirckly -RUN_SINGLE_SIM(MOD_BIRTHS,C_RES[[3]],C_RES[[2]],MORTALITY_SIMULATION,SIM_NUMBER=1,START_OF_SIM=2023) -#Test -sapply(1:1000,function(x){RUN_SINGLE_SIM(MOD_BIRTHS,C_RES[[3]],C_RES[[2]],MORTALITY_SIMULATION,SIM_NUMBER=1,START_OF_SIM=2023)}) - - - -#Run the full simulation in the current year (1), across all simulations x, passing in the demographic, and mortality data. -#TEST <- mclapply(1:10^3,function(x){MORTALITY_SIM(1,x,LIN_CURRENT_DEM,FALSE,Mortality_Rate_Sim )},mc.cores = detectCores()-1) -#TEST[[1]] +SINGLE_PATH_SIM <- function(j){ + C_RES <- RUN_SINGLE_SIM(MOD_BIRTHS,FIRST_PREDICT_YEAR_POPULATION_DATA,START_DEM_DATA,MORTALITY_SIMULATION,SIM_NUMBER=j,START_OF_SIM=2023) + RES <- C_RES[[1]] + for(i in 1:NUM_YEARS_PROJECTED){ + C_RES <- RUN_SINGLE_SIM(MOD_BIRTHS,C_RES[[3]],C_RES[[2]],MORTALITY_SIMULATION,SIM_NUMBER=j,START_OF_SIM=2023) + RES <- rbind(RES,C_RES[[1]]) + } + return(RES) +} +#Run the full simulation across all simulations simulating changes in demographic, and mortality data. +system.time({ +FULL_RESULTS <- mclapply(1:NUM_SIMULATIONS,SINGLE_PATH_SIM,mc.cores = detectCores()-1) +}) +plot(FULL_RESULTS[[2000]]$Population) +TEMP <- sapply(1:length(FULL_RESULTS),function(x){(FULL_RESULTS[[x]] %>% pull(Population))[25] }) +hist(TEMP) diff --git a/Migration_Regression.r b/Migration_Regression.r new file mode 100644 index 0000000..9624c23 --- /dev/null +++ b/Migration_Regression.r @@ -0,0 +1,87 @@ +##########################Model Migration Trends +library(tidyverse) +library(fixest) +library(corrplot) +##Run Regression + DEMOGRAPHIC_DATA <- readRDS("Data/Cleaned_Data/Wyoming_County_Demographic_Data.Rds") + + #Extract the population trend data to connect with demographics (Population,births,deaths) + POP_DATA <- readRDS("Data/Cleaned_Data/Wyoming_County_Population.Rds") +#Identify births, deaths an migration from existing data. +C_YEAR <- 1983 +C_COUNTY <- 'Albany' +POP_DATA %>% filter(Year==C_YEAR,County==C_COUNTY) +sum((DEMOGRAPHIC_DATA %>% filter(Year==C_YEAR,County==C_COUNTY))[,4:5])+34 + +sum((DEMOGRAPHIC_DATA %>% filter(Year==C_YEAR-1,County==C_COUNTY,Age==0))[,4:5]) +sum((DEMOGRAPHIC_DATA %>% filter(Year==C_YEAR,County==C_COUNTY,Age==1))[,4:5]) + +sum((DEMOGRAPHIC_DATA %>% filter(Year==C_YEAR,County==C_COUNTY,Age==0))[,4:5]) + + + + + + + + + + +#############################OTHER TESTING +DATA <- POP_DATA %>% left_join(DEMOGRAPHIC_DATA) %>% filter(!is.na(Births)) +DATA$Age_Group <- NA +DATA <- DATA %>% mutate(Age_Group=ifelse(Age<=5,"Infant",Age_Group)) +DATA <- DATA %>% mutate(Age_Group=ifelse(Age>5 & Age<18,"Child",Age_Group)) +DATA <- DATA %>% mutate(Age_Group=ifelse(Age>=18 & Age<25,"Young_Adult",Age_Group)) +DATA <- DATA %>% mutate(Age_Group=ifelse(Age>=25 & Age<35,"Young_Working_Adult",Age_Group)) +DATA <- DATA %>% mutate(Age_Group=ifelse(Age>=35 & Age<60,"Mid_Adult",Age_Group)) +DATA <- DATA %>% mutate(Age_Group=ifelse(Age>=60,"Retired_Adult",Age_Group)) +DATA %>% filter(Age_Group=="Retired_Adult") +DATA <- DATA %>% ungroup %>% group_by(Year,County,Population,Births,Deaths,Migration,Age_Group) %>% summarize(Num_Male=sum(Num_Male,na.omit=TRUE),Num_Female=sum(Num_Female,na.omit=TRUE)) %>% ungroup +TEMP <- DATA %>% select(-County) %>% pivot_wider(values_from=c(Num_Male,Num_Female),names_from=Age_Group) + +corrplot(cor(TEMP,use="pairwise.complete.obs")) + +REG_TEMP <- DATA %>% pivot_wider(values_from=c(Num_Male,Num_Female),names_from=Age_Group) %>% mutate(Population=Population-Births+Deaths) +REG_TEMP %>% arrange(County,Year) %>% filter(County!='Albany',Year>2015) +#############Looks like Births deaths and migration should be shifted back (or population forward) +POP_DATA %>% group_by(County) %>% arrange(Year) %>% mutate(PREV=Population-Births+Deaths-Migration) %>% arrange(County,Year) %>% filter(Year>2018) +(26500)-501+166+266 +35836+541-184+1137-36209 +(11831-13324)-259+83 +DIFF <- 26519-26165 +DIFF-501+166 +(27380-26633)-413+146 +C_YEAR <-1980 +REG_TEMP %>% filter(Year==C_YEAR-1) +TEMP <- DEMOGRAPHIC_DATA %>% filter(County=='Albany', Year==C_YEAR) +sum(TEMP[1,4:5] ) +TEMP[,4:5] <-DEMOGRAPHIC_DATA %>% filter(County=='Albany', Year==C_YEAR) %>% select(Num_Male,Num_Female)-DEMOGRAPHIC_DATA %>% filter(County=='Albany', Year==C_YEAR-1) %>% select(Num_Male,Num_Female) +TEMP + + +REG_TEMP +REG_TEMP$UPWARD <- ifelse(REG_TEMP$Migration>0,1,0) +REG_TEMP[,5:16] <- log(((REG_TEMP[,5:16]))) +REG_TEMP$Migration <- log(abs(REG_TEMP$Migration)) + +summary(feols(Migration~UPWARD*(Num_Male_Infant+Num_Male_Child+Num_Male_Young_Adult+Num_Male_Young_Working_Adult+Num_Male_Retired_Adult+Num_Female_Infant+Num_Female_Child+Num_Female_Young_Adult+Num_Female_Young_Working_Adult+Num_Female_Retired_Adult)+Population+Population+Year|County,data=REG_TEMP)) +summary(feols(Migration~UPWARD*(Num_Male_Infant+Num_Male_Child+Num_Male_Young_Adult+Num_Male_Young_Working_Adult+Num_Male_Retired_Adult+Num_Female_Infant+Num_Female_Child+Num_Female_Young_Adult+Num_Female_Young_Working_Adult+Num_Female_Retired_Adult)+Population+Population+Year|County,data=REG_TEMP)) + +summary(lm(Migration~.,data=REG_TEMP)) + + + ,Young_Adult=Age>=18,"Child",Age_Group)) +%>% mutate(Child=Age<18,Young_Adult=Age>=18 & Age<35,Mid_Adult=Age>=35 & Age<=60,Retired_Adult=Age>60) %>% group_by(Year,County,Population,Births,Deaths,Migration,Child,Young_Adult,Mid_Adult,Retired_Adult) %>% summarize(Num_Male=sum(Num_Male),Num_Female =sum(Num_Female)) +TEST <- POP_DATA %>% left_join(DEMOGRAPHIC_DATA) %>% filter(!is.na(Births)) %>% pivot_wider(names_from=Age,values_from=c(Num_Male,Num_Female)) +TEST +head(colnames(TEST)) +TEST <- TEST +corrplot(cor(TEST,use="pairwise.complete.obs")) + #Merger the two data sets and drop any records that cannot be used in the regression (this makes the "predict" function output the right number of records) + REG_DATA <- POP_DATA %>% left_join(DEMOGRAPHIC_DATA) %>% filter(!is.na(Births)) + + REG_DATA <- REG_DATA %>% group_by(County) %>% mutate(PREV_MIG=lag(Migration),PREV_TWO_MIG=lag(Migration,2),PREV_POP=lag(Population),PREV_BIRTHS=lag(Births)) %>% ungroup + REG_DATA$County <- factor(REG_DATA$County) +feols((Migration)~(PREV_MIG)+(PREV_TWO_MIG)+PREV_BIRTHS+PREV_POP|Year+County,data=REG_DATA) +REG_DATA %>% filter(!is.na(Births)) diff --git a/Scripts/Data_Load.r b/Scripts/Data_Load.r index 999c8f5..ff3e726 100644 --- a/Scripts/Data_Load.r +++ b/Scripts/Data_Load.r @@ -244,7 +244,7 @@ DEM_DATA <- rbind(DEM_2020,DEM_DATA) %>% ungroup %>% arrange(Year,Age) %>% uniq #########################################Mortality Rate GET_MORTALITY_DATA <- function(FILE,SEX,LOWER_AGE,UPPER_AGE){ - #Create clean moratlity rate data + #Create clean mortality rate data #Data gathered from https://hdpulse.nimhd.nih.gov/data-portal/mortality/table?cod=247&cod_options=cod_15&ratetype=aa&ratetype_options=ratetype_2&race=00&race_options=race_6&sex=2&sex_options=sex_3&age=177&age_options=age_11&ruralurban=0&ruralurban_options=ruralurban_3&yeargroup=5&yeargroup_options=year5yearmort_1&statefips=56&statefips_options=area_states&county=56000&county_options=counties_wyoming&comparison=counties_to_us&comparison_options=comparison_counties&radio_comparison=areas&radio_comparison_options=cods_or_areas NAMES <- c("County","FIPS","Death_Rate","Lower_Rate","Upper_Rate","Deaths","Trend_Category","Trend","Lower_Trend","Upper_Trend") diff --git a/Scripts/Monte_Carlo_Functions.r b/Scripts/Monte_Carlo_Functions.r index 76c6d24..f3260fe 100644 --- a/Scripts/Monte_Carlo_Functions.r +++ b/Scripts/Monte_Carlo_Functions.r @@ -1,5 +1,5 @@ #Uncomment to check the function line by line - #REG_BIRTH_MODEL=MOD_BIRTHS;START_BASIC_DATA=FIRST_PREDICT_YEAR_POPULATION_DATA;START_DETAILED_DATA=START_DEM_DATA;Mortality_Rate_Sim=MORTALITY_SIM;SIM_NUMBER=1;START_OF_SIM=2023 + #REG_BIRTH_MODEL=MOD_BIRTHS;START_BASIC_DATA=FIRST_PREDICT_YEAR_POPULATION_DATA;START_DETAILED_DATA=START_DEM_DATA;Mortality_Rate_Sim=MORTALITY_SIMULATION;SIM_NUMBER=1;START_OF_SIM_YEAR=2023 RUN_SINGLE_SIM <- function(REG_BIRTH_MODEL,START_BASIC_DATA,START_DETAILED_DATA,Mortality_Rate_Sim,SIM_NUMBER,START_OF_SIM_YEAR=2023){ #REG_BIRTH_MODEL: Feols regression object of population model. @@ -37,5 +37,6 @@ RUN_SINGLE_SIM <- function(REG_BIRTH_MODEL,START_BASIC_DATA,START_DETAILED_DATA, NEXT_BASIC_DATA[,"Male_Birth_Group"] <- sum(NEXT_DETAILED_DATA[NEXT_DETAILED_DATA$Age>=18 & NEXT_DETAILED_DATA$Age<=30,"Num_Male"]) NEXT_BASIC_DATA[,"Female_Birth_Group"] <- sum(NEXT_DETAILED_DATA[NEXT_DETAILED_DATA$Age>=18 & NEXT_DETAILED_DATA$Age<=28,"Num_Female"]) NEXT_BASIC_DATA[,"Min_Birth_Group"] <- min(NEXT_BASIC_DATA[,c("Female_Birth_Group","Male_Birth_Group")]) + NEXT_BASIC_DATA <- NEXT_BASIC_DATA[,-10:-11] return(list(START_BASIC_DATA,NEXT_DETAILED_DATA,NEXT_BASIC_DATA)) }