117 lines
5.8 KiB
R
117 lines
5.8 KiB
R
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')
|
|
|
|
DATA_WOMEN_2016 <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Mortality_Rate_and_Pandemic_Data_for_Regression.Rds") %>% filter(Sex=='Female',Year>=2016)
|
|
DATA_MEN_2016 <- readRDS("Data/Cleaned_Data/Mortality_Data/RDS/Mortality_Rate_and_Pandemic_Data_for_Regression.Rds") %>% filter(Sex=='Male',Year>=2016)
|
|
|
|
|
|
#Create time series data
|
|
ST_YEAR <- DATA_MEN %>% pull(Year) %>% min
|
|
|
|
|
|
TS_MEN_US <- DATA_MEN %>% select(Mort_Rate_US) %>% ts(start=ST_YEAR,frequency=1)
|
|
TS_MEN_US_2016 <- DATA_MEN_2016 %>% select(Mort_Rate_US) %>% ts(start=ST_YEAR,frequency=1)
|
|
|
|
TS_MEN_LIN <- DATA_MEN %>% select(Mort_Rate_Lincoln) %>% ts(start=ST_YEAR,frequency=1)
|
|
TS_MEN_LIN_2016 <- DATA_MEN_2016 %>% select(Mort_Rate_Lincoln) %>% ts(start=ST_YEAR,frequency=1)
|
|
|
|
|
|
TS_WOMEN_US <- DATA_WOMEN %>% select(Mort_Rate_US) %>% ts(start=ST_YEAR,frequency=1)
|
|
TS_WOMEN_US_2016 <- DATA_WOMEN_2016 %>% select(Mort_Rate_US) %>% ts(start=ST_YEAR,frequency=1)
|
|
|
|
TS_WOMEN_LIN <- DATA_WOMEN %>% select(Mort_Rate_Lincoln) %>% ts(start=ST_YEAR,frequency=1)
|
|
TS_WOMEN_LIN_2016 <- DATA_WOMEN_2016 %>% select(Mort_Rate_Lincoln) %>% ts(start=ST_YEAR,frequency=1)
|
|
|
|
|
|
TS_PANDEMIC <- DATA_MEN %>% select(WUPI,L_WUPI) %>% ts(start=ST_YEAR,frequency=1)
|
|
TS_PANDEMIC_2016 <- DATA_MEN_2016 %>% select(WUPI,L_WUPI) %>% ts(start=ST_YEAR,frequency=1)
|
|
|
|
|
|
FORECAST_XREG <- TS_PANDEMIC
|
|
FORECAST_XREG_2016 <- TS_PANDEMIC_2016
|
|
|
|
FORECAST_XREG[,] <- 0
|
|
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)
|
|
|
|
#checkresiduals(MOD_US_MEN)
|
|
#plot(forecast(MOD_US_MEN,xreg=FORECAST_XREG))
|
|
|
|
#plot(forecast(MOD_US_WOMEN,xreg=FORECAST_XREG))
|
|
MOD_LIN_WOMEN <- auto.arima(TS_WOMEN_LIN,biasadj=TRUE,xreg=TS_WOMEN_US)
|
|
MOD_LIN_WOMEN_2016 <- auto.arima(TS_WOMEN_LIN_2016,biasadj=TRUE,xreg=TS_WOMEN_US_2016)
|
|
|
|
MOD_LIN_MEN <- auto.arima(TS_MEN_LIN,biasadj=TRUE,xreg=TS_MEN_US)
|
|
MOD_LIN_MEN_2016 <- auto.arima(TS_MEN_LIN_2016,biasadj=TRUE,xreg=TS_MEN_US_2016)
|
|
|
|
#checkresiduals(MOD_LIN_MEN)
|
|
#coeftest(MOD_LIN_WOMEN)
|
|
###Create Location to save models
|
|
if(!exists("SAVE_LOC_MOD")){SAVE_LOC_MOD <-"./Data/Intermediate_Inputs/Age_Mortality_ARIMA_Models/"}
|
|
dir.create(SAVE_LOC_MOD, recursive = TRUE, showWarnings = FALSE)
|
|
|
|
saveRDS(MOD_US_WOMEN,paste0(SAVE_LOC_MOD,"ARIMA_US_Women_Mortality_by_Age.Rds"))
|
|
saveRDS(MOD_US_MEN,paste0(SAVE_LOC_MOD,"ARIMA_US_Men_Mortality_by_Age.Rds"))
|
|
saveRDS(MOD_LIN_WOMEN,paste0(SAVE_LOC_MOD,"ARIMA_Lincoln_Women_Mortality_by_Age.Rds"))
|
|
saveRDS(MOD_LIN_MEN,paste0(SAVE_LOC_MOD,"ARIMA_Lincoln_Men_Mortality_by_Age.Rds"))
|
|
####2016 censured data for validity check
|
|
saveRDS(MOD_US_WOMEN_2016,paste0(SAVE_LOC_MOD,"ARIMA_US_Women_Mortality_by_Age_2016.Rds"))
|
|
saveRDS(MOD_US_MEN_2016,paste0(SAVE_LOC_MOD,"ARIMA_US_Men_Mortality_by_Age_2016.Rds"))
|
|
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))
|
|
|