#Fitting seasonal ARIMA model to US Accidental Deaths series

plot(USAccDeaths)

dif1 <- diff(USAccDeaths)
dif12 <- diff(USAccDeaths,lag=12)
difdif12 <- diff(dif12)

par(mfrow=c(2,3))
acf(dif1)
acf(dif12)
acf(difdif12)
acf(dif1,type="partial")
acf(dif12,type="partial")
acf(difdif12,type="partial")

#Fit ARIMA(0,1,1) x (0,1,1)12 

fit <- arima(USAccDeaths, order=c(0,1,1),
seasonal=list(order=c(0,1,1),period=12))

pred <- predict(fit,n.ahead=12)

# graph the forecasts 
par(mfrow=c(1,1))
plot(c(USAccDeaths, pred$pred), type="n")
lines(1:72, USAccDeaths, lty=1)
lines(73:84, pred$pred + 2*pred$se, lty=2)
lines(73:84, pred$pred - 2*pred$se, lty=2)
points(73:84, pred$pred, pch=3)
title("95% Forecast intervals")