76 lines
3.6 KiB
R
76 lines
3.6 KiB
R
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"))
|
|
|