Building ARIMA model

library(tseries)

package 㤼㸱tseries㤼㸲 was built under R version 3.4.3
    㤼㸱tseries㤼㸲 version: 0.10-42

    㤼㸱tseries㤼㸲 is a package for time series analysis and computational finance.

    See 㤼㸱library(help="tseries")㤼㸲 for details.

library(readr)
library(ggplot2)
library(forecast)

package 㤼㸱forecast㤼㸲 was built under R version 3.4.3This is forecast 8.2 
  Want to meet other forecasters? Join the International Institute of Forecasters:
  http://forecasters.org/

beer1 <- read_csv("D:\\PG Business Analytics\\TSF\\Group Assignment\\beer.csv")

Parsed with column specification:
cols(
  OzBeer = col_double()
)

head(beer1)

ggplot(beer1,aes(x=beer1$OzBeer,y=c(1:72)))+geom_line()

ts<-ts(beer1$OzBeer,frequency = 4)
plot(ts,col="red")

cnt_ma2=ma(beer1$OzBeer,2)
cnt_ma4=ma(beer1$OzBeer,4)
cnt_ma8=ma(beer1$OzBeer,8)
ggplot() +
  geom_line(data = beer1, aes(x = c(1:72), y = ts, colour = "Quarterly Actual Sales")) +
  geom_line(data = beer1, aes(x = c(1:72), y = cnt_ma2, colour = "Hal-yearly Moving Average"))  +
  geom_line(data = beer1, aes(x = c(1:72), y = cnt_ma4, colour = "Yearly Moving Average"))  +
  geom_line(data = beer1, aes(x = c(1:72), y = cnt_ma8, colour = "Two Year Moving Average"))  +
  ylab('Beer Sales')

kpss.test(ts)

p-value smaller than printed p-value

    KPSS Test for Level Stationarity

data:  ts
KPSS Level = 3.0456, Truncation lag parameter = 1, p-value = 0.01

So the series is not stationary

ndiffs(ts)

[1] 1

dif = diff(ts,differences = 1)
kpss.test(dif)

p-value greater than printed p-value

    KPSS Test for Level Stationarity

data:  dif
KPSS Level = 0.036853, Truncation lag parameter = 1, p-value = 0.1

After differenciating it once, we see that the series is stationary.

Another test for stationarity, Dickey-Fuller test

adf.test(ts)

    Augmented Dickey-Fuller Test

data:  ts
Dickey-Fuller = -0.73068, Lag order = 4, p-value = 0.9633
alternative hypothesis: stationary

hence we accept NULL hypothesis

adf.test(dif)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  dif
Dickey-Fuller = -6.2957, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

this proves that d = 1.

Part A

Holt-Winters method

fit = hw(ts,seasonal="multiplicative")
summary(fit)

Forecast method: Holt-Winters' multiplicative method

Model Information:
Holt-Winters' multiplicative method 

Call:
 hw(y = ts, seasonal = "multiplicative") 

  Smoothing parameters:
    alpha = 0.0643 
    beta  = 0.0406 
    gamma = 1e-04 

  Initial states:
    l = 256.4184 
    b = 0.022 
    s=1.1734 0.9247 0.8768 1.0251

  sigma:  0.0298

     AIC     AICc      BIC 
650.6460 653.5492 671.1360 

Error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
Training set 1.699405 9.183581 7.517402 0.4673523 2.382933 0.5551514 -0.1745673

Forecasts:
      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
19 Q1       467.4689 449.6426 485.2953 440.2059 494.7320
19 Q2       404.1306 388.6368 419.6244 380.4349 427.8264
19 Q3       430.8064 414.1248 447.4881 405.2941 456.3188
19 Q4       552.4171 530.6915 574.1427 519.1906 585.6435
20 Q1       487.6920 468.0895 507.2945 457.7126 517.6715
20 Q2       421.4266 404.0029 438.8502 394.7794 448.0738
20 Q3       449.0489 429.8322 468.2656 419.6595 478.4384
20 Q4       575.5641 549.9233 601.2049 536.3498 614.7783

So we get a AIC of 650.

plot(ts,col="blue")
lines(fitted(fit),col="red")

forecast(fit,4*2)

      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
19 Q1       467.4689 449.6426 485.2953 440.2059 494.7320
19 Q2       404.1306 388.6368 419.6244 380.4349 427.8264
19 Q3       430.8064 414.1248 447.4881 405.2941 456.3188
19 Q4       552.4171 530.6915 574.1427 519.1906 585.6435
20 Q1       487.6920 468.0895 507.2945 457.7126 517.6715
20 Q2       421.4266 404.0029 438.8502 394.7794 448.0738
20 Q3       449.0489 429.8322 468.2656 419.6595 478.4384
20 Q4       575.5641 549.9233 601.2049 536.3498 614.7783

plot(forecast(fit,4*2))

This concludes part A

Part B

ARIMA model

plot(dif,col = "blue")
abline(lm(dif~time(dif)),col = "red")

decom<-decompose(ts)
plot(decom,col = "blue")

#ma4<-ma(dif,4)
par(mfrow=c(2,1))
Acf(diff(ts), main='')
Pacf(diff(ts), main='')

from the above charts we can conclude that p = 3 and q = 2. but we must test the same with our model

auto.arima((ts), seasonal=TRUE)

Series: (ts) 
ARIMA(2,1,1)(0,1,1)[4] 

Coefficients:
          ar1      ar2      ma1     sma1
      -0.5280  -0.2791  -0.5445  -0.6170
s.e.   0.1922   0.1700   0.1852   0.1044

sigma^2 estimated as 109.6:  log likelihood=-252.08
AIC=514.16   AICc=515.14   BIC=525.18

fit<-arima(ts,order =c(2,1,1),seasonal = list(order =c(0,1,1),period =4))
fit

Call:
arima(x = ts, order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))

Coefficients:
          ar1      ar2      ma1     sma1
      -0.5280  -0.2791  -0.5445  -0.6170
s.e.   0.1922   0.1700   0.1852   0.1044

sigma^2 estimated as 103.1:  log likelihood = -252.08,  aic = 514.16

tsdisplay(residuals(fit), lag.max=15, main='Seasonal Model Residuals')

Box.test(residuals(fit),lag=4,type="Ljung")

    Box-Ljung test

data:  residuals(fit)
X-squared = 0.66143, df = 4, p-value = 0.956

pred <- predict(fit, n.ahead = 4*2)
ts.plot(ts,pred$pred, log = "y", lty = c(1,3),ylabel ="Passenger",main="Forecasted Series",col="red")

NAs introduced by coercion