library(tidyverse)
library(forecast)
library(lmtest)
####################
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_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_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"))