##########################Model Population Trends
library(fixest)
library(tidyverse)
readRDS("Data/Cleaned_Data/Wyoming_County_Population.Rds")
REG_DATA <- readRDS("Data/Cleaned_Data/Wyoming_County_Population.Rds") %>% left_join(readRDS("Data/Cleaned_Data/Wyoming_County_Demographic_Data.Rds") %>%  group_by(County,Year) %>% summarize(POP2=sum(Num_Male+Num_Female))) %>%left_join( readRDS("Data/Cleaned_Data/Wyoming_County_Demographic_Data.Rds") %>% mutate(BIRTH_POP=ifelse(Age==0,1,0)) %>% group_by(County,Year,BIRTH_POP) %>% summarize(BIRTH_POP=sum(Num_Female))) %>% mutate(LN=ifelse(County=="Lincoln",1,0),Pop=Population-Births)
TEMP <- readRDS("Data/Cleaned_Data/Wyoming_County_Demographic_Data.Rds")  %>% mutate(Num_Male=log(Num_Male),Num_Female=log(Num_Female))
#TEMP <-	     pivot_longer(TEMP,c("Num_Male","Num_Female"),names_to="Sex",values_to="Number") %>% mutate(Sex=ifelse(Sex=="Num_Male","Male","Female"))
TEMP <- pivot_wider(TEMP,names_from=Age,values_from=c(Num_Male,Num_Female),names_prefix="Age_") %>% unique
#Establish good bounds for the birth groups
REG_DATA_PRLIM <- readRDS("Data/Cleaned_Data/Wyoming_County_Population.Rds") %>% mutate(Population=Population-Births+Deaths+Migration) %>% select(-Deaths,-Migration) %>% left_join(TEMP) %>% select(-Population,-Num_Male_Age_0,-Num_Female_Age_0)
  summary(lm(Births~. ,data=REG_DATA_PRLIM))
##Run Regression 
  #Pull in Demographic data, and create categories for key groups, male/female population with high fertility, children under one and two (but not zero)
  DEMOGRAPHIC_DATA <- readRDS("Data/Cleaned_Data/Wyoming_County_Demographic_Data.Rds") %>% mutate(Male_Window=Age>=18 & Age<=30,Female_Window=Age>=18 & Age<=28,Under_Two=Age<=2 & Age!=0,Under_One=Age<=1 & Age!=0) %>% group_by(County,Year) %>% summarize(Female_Birth_Group=sum(Num_Female*Female_Window),Male_Birth_Group=sum(Num_Male*Male_Window),Under_Two=sum(Under_Two*(Num_Male+Num_Female)),Under_One=sum(Under_One*(Num_Male+Num_Female)),Min_Birth_Group=ifelse(Female_Birth_Group<Male_Birth_Group,Female_Birth_Group,Male_Birth_Group))
  POP_DATA <- readRDS("Data/Cleaned_Data/Wyoming_County_Population.Rds") %>% mutate(LN=ifelse(County=="Lincoln",1,0))
  REG_DATA <- POP_DATA %>% full_join(DEMOGRAPHIC_DATA) %>% filter(!is.na(Births))
  REG_DATA <- REG_DATA %>% group_by(County) %>% mutate(PREV_BIRTH=lag(Births),PREV_TWO_BIRTH=lag(Births,2)) %>% ungroup %>% filter(!is.na(PREV_TWO_BIRTH),!is.na(Min_Birth_Group))
REG_DATA
  MOD_BIRTHS<- feols(log(Births)~log(PREV_BIRTH)+log(PREV_TWO_BIRTH)+log(Min_Birth_Group)+Year*County,cluster=~Year, data=REG_DATA )
  #acf(RES_DATA %>% pull(RESID))
  #pacf(RES_DATA %>% pull(RESID))

TEST <- REG_DATA %>% filter(LN==1) %>% filter(Year==2022) %>% mutate(Year=2022)
C_PREDICT <- predict(MOD_BIRTHS,TEST,interval = "prediction",level=0.95)
PRED_MEAN <- C_PREDICT$fit
SE_PRED <- (C_PREDICT$ci_high-C_PREDICT$ci_low)/3.92

YEAR <- 2025
NUM_SIMS <- 10000
BIRTHS <- round(exp(rnorm(NUM_SIMS,mean=PRED_MEAN,sd=SE_PRED)))
MALE <- sapply(1:NUM_SIMS,function(x){ rbinom(1,BIRTHS[x],prob=0.5)})
RES <- cbind(rep(YEAR,NUM_SIMS),rep(0,NUM_SIMS),MALE,BIRTHS-MALE) %>% as_tibble
colnames(RES) <- c("Year","Age","Num_Male","Num_Female")
RES
##NOTE TEST[[1]] Comming from death simulation script as a test in order to feedback into births.
CVAL <- TEST[[1]]
CVAL 
Male_Birth_Group <- sum(CVAL[CVAL$Age>=18 & CVAL$Age<=30,] %>% pull(Num_Male))
Female_Birth_Group <- sum(CVAL[CVAL$Age>=18 & CVAL$Age<=28,] %>% pull(Num_Female))
Min_Birth_Group <- min(Male_Birth_Group,Female_Birth_Group )
YEAR <- YEAR+1