diff --git a/Scripts/2F_Create_Age_Sex_Specfic_Mortality_Trends.r b/Scripts/2F_Create_Age_Sex_Specfic_Mortality_Trends.r index 9a71ee6..7610da9 100644 --- a/Scripts/2F_Create_Age_Sex_Specfic_Mortality_Trends.r +++ b/Scripts/2F_Create_Age_Sex_Specfic_Mortality_Trends.r @@ -17,6 +17,15 @@ REG_DATA_2016 <- REG_DATA %>% filter(Year>=2016) MOD <- feols(Age_.[0:85]~US_Adj_Death_Rate+Sex*Year+WUPI+L_WUPI,REG_DATA, data.save = TRUE) MOD_2016 <- feols(Age_.[0:85]~US_Adj_Death_Rate+Sex*Year+WUPI+L_WUPI,REG_DATA_2016, data.save = TRUE) +###Save results + if(!exists("SAVE_FIG_LOC")){SAVE_FIG_LOC <- "./Results/Age_Mortality_Regression"} + dir.create(SAVE_FIG_LOC , recursive = TRUE, showWarnings = FALSE) + REG_TABLE_LOC <- paste0(SAVE_FIG_LOC,"/Morality_by_Age_Regressions.png") + REG_TABLE_TEX_LOC <- paste0(SAVE_FIG_LOC,"/Mortality_by_Age_Regressions.tex") + +DICT <- c("US_Adj_Death_Rate"="US Mortality Rate","SexMale"="Male","WUPI"="Pandemic Index","L_WUPI"="Pandemic Index (1 Year)") +etable(MOD[[1]],MOD[[19]],MOD[[31]],MOD[[61]],MOD[[81]],dict=DICT,style.tex=style.tex(yesNo="$\\checkmark$"),file=REG_TABLE_TEX_LOC ,export=REG_TABLE_LOC,replace=TRUE) + ###Simulate each age-sex death rate over time with the models #########When project far into the future some death rate values become negative. Make bounds to limit the forecast to a reasonable range. In this case I select half of the historic minimum, or double the historic maximum as upper an lower bounds in the study period. diff --git a/Scripts/2G_Single_Age_Sex_ARIMA_Models.r b/Scripts/2G_Single_Age_Sex_ARIMA_Models.r index cc41c37..f23094b 100644 --- a/Scripts/2G_Single_Age_Sex_ARIMA_Models.r +++ b/Scripts/2G_Single_Age_Sex_ARIMA_Models.r @@ -1,6 +1,8 @@ library(tidyverse) library(forecast) library(lmtest) +library(texreg) + #################### DATA_WOMEN <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Mortality_Rate_and_Pandemic_Data_for_Regression.Rds") %>% filter(Sex=='Female') DATA_MEN <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Mortality_Rate_and_Pandemic_Data_for_Regression.Rds") %>% filter(Sex=='Male') @@ -39,11 +41,13 @@ FORECAST_XREG_2016[,] <- 0 MOD_US_WOMEN <- auto.arima(TS_WOMEN_US,lambda=0,biasadj=TRUE,xreg=TS_PANDEMIC) +MOD_US_WOMEN MOD_US_WOMEN_2016 <- auto.arima(TS_WOMEN_US_2016,lambda=0,biasadj=TRUE,xreg=TS_PANDEMIC_2016) #checkresiduals(MOD_US_WOMEN) MEN_XREG <- cbind(TS_WOMEN_US,TS_PANDEMIC) MEN_XREG_2016 <- cbind(TS_WOMEN_US_2016,TS_PANDEMIC_2016) +MOD_US_MEN <- auto.arima(TS_MEN_US,lambda=0,biasadj=TRUE,xreg=MEN_XREG) MOD_US_MEN <- auto.arima(TS_MEN_US,lambda=0,biasadj=TRUE,xreg=MEN_XREG) MOD_US_MEN_2016 <- auto.arima(TS_MEN_US_2016,lambda=0,biasadj=TRUE,xreg=MEN_XREG_2016) @@ -73,3 +77,40 @@ saveRDS(MOD_US_MEN_2016,paste0(SAVE_LOC_MOD,"ARIMA_US_Men_Mortality_by_Age_2016. saveRDS(MOD_LIN_WOMEN_2016,paste0(SAVE_LOC_MOD,"ARIMA_Lincoln_Women_Mortality_by_Age_2016.Rds")) saveRDS(MOD_LIN_MEN_2016,paste0(SAVE_LOC_MOD,"ARIMA_Lincoln_Men_Mortality_by_Age_2016.Rds")) +###Mortality ARIMA results save + + if(!exists("SAVE_LOC_ARIMA_FIGURES")){SAVE_LOC_ARIMA_FIGURES <-"./Results/Mortality_ARIMA/"} + dir.create(SAVE_LOC_ARIMA_FIGURES, recursive = TRUE, showWarnings = FALSE) +################Figures +png(paste0(SAVE_LOC_ARIMA_FIGURES,"US_Men_Mortality_ARIMA_Residual_Checks.png"), res = 600, height = 12, width=16, units = "in") + checkresiduals(MOD_US_MEN) +dev.off() + +png(paste0(SAVE_LOC_ARIMA_FIGURES,"US_Women_Mortality_ARIMA_Residual_Checks.png"), res = 600, height = 12, width=16, units = "in") + checkresiduals(MOD_US_WOMEN) +dev.off() + +png(paste0(SAVE_LOC_ARIMA_FIGURES,"Lincoln_County_Men_Mortality_ARIMA_Residual_Checks.png"), res = 600, height = 12, width=16, units = "in") + checkresiduals(MOD_LIN_MEN) +dev.off() + +png(paste0(SAVE_LOC_ARIMA_FIGURES,"Lincoln_County_Women_Mortality_ARIMA_Residual_Checks.png"), res = 600, height = 12, width=16, units = "in") + checkresiduals(MOD_LIN_WOMEN) +dev.off() + + + +##################Tables +DICT <- list("ar1"="Autoregressive (1 Year)","ar2"="Autoregressive (2 Year)","ma1"="Moving Average (1 Year)","ma2"="Moving Average (2 Year)","WUPI"="Pandemic Index","L_WUPI"="Pandemic Index (1 Year)","intercept"="Average Mortality","Women Mortality Rate"="TS\\_WOMEN\\_US" ) +NOTE <- list("The US Mortality rate variable is tied to the sex of the given model.","Lambda set to zero for all model which log-transforms the results.") +FILE_NAME <- "Mortality_ARIMA_Tables" +REG_TABLE_LOC <- paste0(SAVE_LOC_ARIMA_FIGURES,FILE_NAME,".tex") +sink(REG_TABLE_LOC) +cat("\\documentclass[border=0pt]{article}","\n","\\pagestyle{empty}","\n","\\usepackage{booktabs,dcolumn}","\n","\\begin{document}") +texreg(l=list(MOD_US_WOMEN,MOD_US_MEN,MOD_LIN_WOMEN,MOD_LIN_MEN),digits=4,custom.model.names=c("\\textbf{U.S. Women}","\\textbf{U.S. Men}","\\textbf{Lincoln Women}","\\textbf{Lincoln Men}"),table=FALSE,use.packages=FALSE,booktabs=TRUE,dcolumn=TRUE,caption.above=TRUE) +cat("\n","\\end{document}") +sink() +system(paste0("pdflatex ",REG_TABLE_LOC)) +system(paste0("pdfcrop --margins '5 5 5 5' ",FILE_NAME,".pdf ",SAVE_LOC_ARIMA_FIGURES,FILE_NAME,".pdf")) +file.remove(list.files(pattern=FILE_NAME)) +