Description

I chose to evaluate the percentage of all US deaths due to Pneumonia and Influenza. The data can be retrieved from the following CDC website:

Mortality Data

The dataset begins at week 40 in 2015 and continues until week 41 in 2019, for a total of 210 observations.

Motivation

I have a background in Finance and spent many years as an equity analyst. I primarily covered the healthcare sector, with a focus on healthcare services and medical devices. Trading in healthcare stocks sometimes requires an estimate of patient visits for the current period, but since I did not have a background in Statistics I relied mostly on conversations with hospital administrators. I often wondered if there was a better way to forecast, but until now I have never attempted to do so with time series models.

Unfortunately there is no weekly or monthly data available for healthcare utilization, and so I chose this dataset as a proxy.

Transformation

I began by viewing the time series plot of the dataset, as shown below. We can see that there is a clear seasonal pattern, with peaks during the winter months and troughs during the summer. We can also see that the variance from year to year is not constant, with a particularly strong flu season 2 years ago. To determine if a transformation would help stabilize the variance I ran a Box-Cox test, which suggested that an inverse transformation would be useful (value of -0.9999). We can see from the inverse plot below that some of the variance is lessened.

library(astsa)
library(forecast)
library(tseries)
library(TSA)
library(knitr)
library(kableExtra)

#Import Dataset
PIdths <- read.csv("~/Baruch/Baruch Class Projects/9701 Time Series/PIdeaths.csv")
PIdths = ts(PIdths) #convert to time series

#Plot the timeseries
plot.ts(PIdths,main="Percentage of all deaths due to pneumonia and influenza",ylab="% of Deaths")

#Determine whether transformation is necessary
BoxCox.lambda(PIdths)

## [1] -0.9999242

#result of -.9999 indicates inverse transformation 

InvPIdths = PIdths^-1 # Convert data to inverse

#Plot the inverse timeseries
plot.ts(InvPIdths,main="Inverse % of all deaths due to pneumonia and influenza",ylab="% of Deaths")

Next, I viewed both the ACF and PACF of the inversed data. We can see the seasonal pattern in the ACF, and a spike in the PACF at 52 weeks.

#Plot the sample ACF and PACF
acf(PIdths,100, main="")

pacf(PIdths,100)

Stationarity

To normalize the seasonal pattern I differenced the weekly data by 52. I then evaluated the Augmented Dickey-Fuller test along with the ndiffs function, which suggested that another difference would help make the data stationary. We can see from the resulting plot that the time series does appear to have a constant mean and relatively stable variance. I again evaluated the Augmented Dickey-Fuller test, this time on the differenced data (52 week difference followed by 1 week difference), and with a p-value of <0.01 we can reject the null hypothesis that the data is not stationary.

# 52 week data.  Clear seasonal trend, so take 52 week difference
dInvPIdths = diff(PIdths,52)

#Plot the 52 week differenced timeseries
plot.ts(dInvPIdths,main="Inverse % of all deaths due to P&I, 52 week diff",ylab="% of Deaths")

#Test whether data is stationary or if differencing is necessary
adf.test(dInvPIdths) # Augmented Dickey-Fuller test

    Augmented Dickey-Fuller Test

data:  dInvPIdths
Dickey-Fuller = -3.5047, Lag order = 5, p-value = 0.04429
alternative hypothesis: stationary

ndiffs(dInvPIdths)  # confirm with ndiffs.  Recommends 1 difference

[1] 1

ddInvPIdths = diff(dInvPIdths,1) #differece the data again

#Plot the double differenced timeseries
plot.ts(ddInvPIdths,main="Inverse % of all deaths due to P&I, diff 52 and 1",ylab="% of Deaths")

# plot looks stable with constant variance and mean

#Test whether data is stationary or if differencing is necessary
adf.test(ddInvPIdths) # Augmented Dickey-Fuller test

    Augmented Dickey-Fuller Test

data:  ddInvPIdths
Dickey-Fuller = -5.4507, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary

ndiffs(ddInvPIdths)  # confirm with ndiffs

[1] 0

# confirmed stationary

ACF and PACF review

I then viewed the ACF and PACF plots of the differenced data. We can see that the seasonal component (52 week lags) suggest a possible seasonal MA(1) model, since the ACF cuts off after lag 1 and PACF tails off. The ACF and PACF of the non-seasonal component do not have clear cut-off points, which seems to suggest that an ARMA model would be appropriate.

acf(ddInvPIdths,200, main="") # non seasonal ACF tails off, seasonal spike at 52, 54 weeks, cuts off

pacf(ddInvPIdths,200) #non seasonal PACF tails off, #seasonal PACF tails off

#Since the seasonal ACF cuts off and seaonal PACF tails off, use seasonal MA model
#Since non-seasonal ACF and PACF tail off, use non-seaonal ARMA
#Use 1 seasonal difference at 52 weeks and 1 non-seaonal difference

Evaluation

Based on the ACF/PACF analysis I began by evaluating seasonal ARIMA models with 1 non-seasonal difference, 1 seasonal difference, a seasonal MA(1) component, and an ARMA non-seasonal component. I began by evaluating the simplest model first, an ARIMA$(1,1,1)$ x $(0,1,1)_{52}$. I then continued evaluating different possibilities until I was satisfied with the residual analysis and AICc/BIC values.

#try simplest model first
fit1 <- sarima(InvPIdths,1,1,1,0,1,1,52) #diagnostics OK, but ACF of some resids significant

fit1$AICc 

[1] -5.658188

fit1$BIC 

[1] -5.59998

Based on the results, I chose the SARIMA$(1,1,5)$ x $(0,1,1)_{52}$ model.

fit3 <- sarima(InvPIdths,1,1,5,0,1,1,52) #diagnostics perfect

fit3$ttable

     Estimate     SE  t.value p.value
ar1    0.8360 0.0730  11.4454  0.0000
ma1   -1.1995 0.1039 -11.5449  0.0000
ma2    0.4070 0.1308   3.1127  0.0022
ma3   -0.0703 0.1337  -0.5260  0.5997
ma4    0.1283 0.1140   1.1259  0.2620
ma5   -0.2450 0.0794  -3.0851  0.0024
sma1  -0.5311 0.1533  -3.4650  0.0007

Since the ma3 and ma4 coefficients were not significant, I eliminated them from the final model.

Residual Diagnostics

We can see that the final model’s residuals have no clear pattern, have no significant ACF values, and all have Ljung-Box p-values over 0.05.

Final Model

The seasonal autoregressive integrated moving average model is represented by:

$$SARIMA(p,d,q)\quad x \quad (P,D,Q)_s$$

$$\Phi_p(B^s)\phi(B)\nabla_s^D\nabla^dX_t = \delta + \Theta_Q(B^s)\theta(B)w_t$$

Where:
$\Phi_P(B^s)= 1-\Phi_1B^s-\Phi_2B^{2s}-...-\Phi_PB^{Ps}$
$\phi(B)=1-\phi_1B-\phi_2B^2-...-\phi_PB^P$
$\nabla_s^D = (1-B^s)^D$
$\nabla^d=(1-B)^d$
$\Theta_Q(B^s)=1+\Theta_1B^s+\Theta_2B^{2s}+...+\Theta_QB^{Qs}$
$\theta(B)=1+\theta_1B+\theta_2B^2+...+\theta_qB^q$

fit4 <- sarima(InvPIdths,1,1,5,0,1,1,52, fixed=c(NA,NA,NA,0,0,NA,NA))

fit4$ttable

     Estimate     SE  t.value p.value
ar1    0.8463 0.0724  11.6970   0e+00
ma1   -1.2164 0.1073 -11.3325   0e+00
ma2    0.4065 0.0879   4.6228   0e+00
ma5   -0.1755 0.0464  -3.7794   2e-04
sma1  -0.5230 0.1514  -3.4541   7e-04