Preample

library("forecast")

Coke Revenue

Coca-Cola Net Revenue, Quarterly, 1994-2009.

data <- read.csv("http://ptrckprry.com/course/forecasting/data/cokerev.csv")
date <- as.Date(data$date)
time <- 1:length(date)
revenue <- data$revenue

Time Series Plots

plot(date, revenue, type="l", col=2)

ACF, PACF

Original Series

par(mfrow=c(1,2))
Acf(revenue)
Pacf(revenue)

First Difference

par(mfrow=c(1,2))
Acf(diff(revenue))
Pacf(diff(revenue))

Seasonal Difference of First Difference

par(mfrow=c(1,2))
Acf(diff(diff(revenue), 4))
Pacf(diff(diff(revenue), 4))

Seasonal Difference of Second Difference

par(mfrow=c(1,2))
Acf(diff(diff(diff(revenue)), 4))
Pacf(diff(diff(diff(revenue)), 4))

Second Seasonal Difference of First Difference

par(mfrow=c(1,2))
Acf(diff(diff(diff(revenue), 4), 4))
Pacf(diff(diff(diff(revenue), 4), 4))

AICC

Based on ACFs, PACFs, choose d=1, D=1. Next, we need (p,q) and (P,Q)

d <- 1
D <- 1

for (p in 0:2) {
    for (q in 0:2) {
        for (P in 0:2) {
            for (Q in 0:2) {
                fit <- list(aicc=NA)
                try({
                    fit <- Arima(revenue, c(p,d,q),
                                 seasonal=list(order=c(P,D,Q), period=4),
                                 method="ML")
                })
                cat("ARIMA(", p, ",", d, ",", q, ")",
                    "(", P, ",", D, ",", Q, ")[4]",
                    " : ", fit$aicc, "\n", sep="")
            }
        }
    }
}
ARIMA(0,1,0)(0,1,0)[4] : 907.8622
ARIMA(0,1,0)(0,1,1)[4] : 885.4825
ARIMA(0,1,0)(0,1,2)[4] : 887.6426
ARIMA(0,1,0)(1,1,0)[4] : 897.3474
ARIMA(0,1,0)(1,1,1)[4] : 887.6569
ARIMA(0,1,0)(1,1,2)[4] : 888.9501
ARIMA(0,1,0)(2,1,0)[4] : 891.3722
ARIMA(0,1,0)(2,1,1)[4] : 888.948
ARIMA(0,1,0)(2,1,2)[4] : 890.8344
ARIMA(0,1,1)(0,1,0)[4] : 909.2511
ARIMA(0,1,1)(0,1,1)[4] : 887.2063
ARIMA(0,1,1)(0,1,2)[4] : 889.4555
ARIMA(0,1,1)(1,1,0)[4] : 898.7736
ARIMA(0,1,1)(1,1,1)[4] : 889.4649
ARIMA(0,1,1)(1,1,2)[4] : 891.0353
ARIMA(0,1,1)(2,1,0)[4] : 893.601
ARIMA(0,1,1)(2,1,1)[4] : 891.0112
ARIMA(0,1,1)(2,1,2)[4] : 893.124
ARIMA(0,1,2)(0,1,0)[4] : 911.2511
ARIMA(0,1,2)(0,1,1)[4] : 888.5856
ARIMA(0,1,2)(0,1,2)[4] : 890.8874
ARIMA(0,1,2)(1,1,0)[4] : 900.8786
ARIMA(0,1,2)(1,1,1)[4] : 890.9047
ARIMA(0,1,2)(1,1,2)[4] : 892.4356
ARIMA(0,1,2)(2,1,0)[4] : 895.6441
ARIMA(0,1,2)(2,1,1)[4] : 892.5115
ARIMA(0,1,2)(2,1,2)[4] : 894.6675
ARIMA(1,1,0)(0,1,0)[4] : 909.2912
ARIMA(1,1,0)(0,1,1)[4] : 887.3139
ARIMA(1,1,0)(0,1,2)[4] : 889.5596
ARIMA(1,1,0)(1,1,0)[4] : 898.8511
ARIMA(1,1,0)(1,1,1)[4] : 889.5702
ARIMA(1,1,0)(1,1,2)[4] : 891.0995
ARIMA(1,1,0)(2,1,0)[4] : 893.6083
ARIMA(1,1,0)(2,1,1)[4] : 891.0788
ARIMA(1,1,0)(2,1,2)[4] : 893.1617
ARIMA(1,1,1)(0,1,0)[4] : 905.0479
ARIMA(1,1,1)(0,1,1)[4] : 886.9345
ARIMA(1,1,1)(0,1,2)[4] : 889.2028
ARIMA(1,1,1)(1,1,0)[4] : 896.7257
ARIMA(1,1,1)(1,1,1)[4] : 889.2308
ARIMA(1,1,1)(1,1,2)[4] : 890.7414
ARIMA(1,1,1)(2,1,0)[4] : 892.3508
ARIMA(1,1,1)(2,1,1)[4] : 890.8749
ARIMA(1,1,1)(2,1,2)[4] : 892.9558
ARIMA(1,1,2)(0,1,0)[4] : 906.0811
ARIMA(1,1,2)(0,1,1)[4] : 889.279
ARIMA(1,1,2)(0,1,2)[4] : 891.6079
ARIMA(1,1,2)(1,1,0)[4] : 899.112
ARIMA(1,1,2)(1,1,1)[4] : 891.6447
ARIMA(1,1,2)(1,1,2)[4] : 893.1306
ARIMA(1,1,2)(2,1,0)[4] : 894.4997
ARIMA(1,1,2)(2,1,1)[4] : 893.3428
ARIMA(1,1,2)(2,1,2)[4] : 895.4136
ARIMA(2,1,0)(0,1,0)[4] : 911.4095
ARIMA(2,1,0)(0,1,1)[4] : 888.6688
ARIMA(2,1,0)(0,1,2)[4] : 891.0068
ARIMA(2,1,0)(1,1,0)[4] : 900.9236
ARIMA(2,1,0)(1,1,1)[4] : 891.0137
ARIMA(2,1,0)(1,1,2)[4] : 892.633
ARIMA(2,1,0)(2,1,0)[4] : 895.7048
ARIMA(2,1,0)(2,1,1)[4] : 892.6311
ARIMA(2,1,0)(2,1,2)[4] : 894.7932
ARIMA(2,1,1)(0,1,0)[4] : 905.1105
ARIMA(2,1,1)(0,1,1)[4] : 889.28
ARIMA(2,1,1)(0,1,2)[4] : 891.6094
ARIMA(2,1,1)(1,1,0)[4] : 898.8752
ARIMA(2,1,1)(1,1,1)[4] : 891.6461
ARIMA(2,1,1)(1,1,2)[4] : 893.1331
ARIMA(2,1,1)(2,1,0)[4] : 894.4648
ARIMA(2,1,1)(2,1,1)[4] : 893.3457
ARIMA(2,1,1)(2,1,2)[4] : 895.4183
ARIMA(2,1,2)(0,1,0)[4] : 902.8535
ARIMA(2,1,2)(0,1,1)[4] : 890.469
ARIMA(2,1,2)(0,1,2)[4] : 892.8883
ARIMA(2,1,2)(1,1,0)[4] : 897.7937
ARIMA(2,1,2)(1,1,1)[4] : 892.9219
ARIMA(2,1,2)(1,1,2)[4] : 894.6571
ARIMA(2,1,2)(2,1,0)[4] : 894.8784
ARIMA(2,1,2)(2,1,1)[4] : 894.9
ARIMA(2,1,2)(2,1,2)[4] : 897.1286

Estimation

(fit <- Arima(revenue, c(0,1,0), seasonal=list(order=c(0,1,1), period=4)))
Series: revenue 
ARIMA(0,1,0)(0,1,1)[4]                    

Coefficients:
         sma1
      -0.7489
s.e.   0.0979

sigma^2 estimated as 67137:  log likelihood=-440.64
AIC=885.28   AICc=885.48   BIC=889.57

Diagnostic Checking

resid <- residuals(fit)

Time Series Plot of Residuals

plot(date, resid, type="l", col=2)

Should we have taken logs?

ACF, PACF of Residuals

par(mfrow=c(1,2))
Acf(resid)
Pacf(resid)

Ljung-Box Tests

Box.test(resid, lag=12, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 5.4142, df = 7, p-value = 0.6096
Box.test(resid, lag=24, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 16.657, df = 19, p-value = 0.6131
Box.test(resid, lag=36, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 32.212, df = 31, p-value = 0.4065
Box.test(resid, lag=48, type="Ljung-Box", fitdf=5)

    Box-Ljung test

data:  resid
X-squared = 39.39, df = 43, p-value = 0.6286

Forecasting

Gaussian

plot(forecast(fit, h=20))

Bootstrap

plot(forecast(fit, h=20, boot=TRUE))

Comparison With Non-Seasonal Model

fit0 <- Arima(revenue, c(0,1,1))

Gaussian

plot(forecast(fit0, h=20))

Bootstrap

plot(forecast(fit0, h=20, boot=TRUE))

hist(resid(fit0), 20)