# 
# Computes the maximum likelihood ARIMA coefficients for Series 8. Beer Shipments-Four-Week Totals (in thousands of units)
# 
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

# get the data into an R vector:
# Read downwards, left to right.
# Year 1 Year 2 Year 3 Year 4
# 
Y <- c(
   18.7050,   19.5980,   20.8160,   21.2600,
   20.2320,   21.4630,   23.7430,   24.1090,
   20.4670,   23.2870,   25.1520,   26.3200,
   22.1230,   24.0650,   28.8040,   27.7010,
   25.0360,   27.4470,   31.1580,   34.5020,
   26.8390,   30.4130,   31.5400,   33.2970,
   29.6400,   32.3070,   32.8490,   31.2520,
   30.9350,   32.9740,   33.7480,   35.1730,
   28.2780,   29.9730,   31.9100,   36.2070,
   24.2350,   23.9860,   27.6090,   31.5110,
   22.3700,   26.9530,   25.1700,   28.5600,
   21.2240,   24.2500,   24.0400,   26.8280,
   21.0610,   23.5180,   25.3680,   26.6600
) 
n <- length(Y)

# we plot the original data:
postscript("../../WriteUp/Graphics/Chapter6/prob_6_13_data_plt.eps", onefile=FALSE, horizontal=FALSE)
plot(1:n, Y, type="l", xlab="time", ylab="beer shipments") 
dev.off()

# consider the autocorrelation function directly 
postscript("../../WriteUp/Graphics/Chapter6/prob_6_13_data_direct_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(Y)
dev.off()

# take a seasonal difference to make the data stationary
Yd <- diff(Y,lag=1,differences=1)
nd <- length(Yd)
plot(1:nd, Yd, type="l", xlab="time", ylab="first seasonal difference of beer shipments") 

# plot the autocorrelation function of this first difference 
postscript("../../WriteUp/Graphics/Chapter6/prob_6_13_data_d1_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(Yd)
dev.off()

# based on the sample autocorrelation above lets assume a SAR(1) model of Yd
# this is a SARIMA (0,1,0)(1,0,0)_{12} model on Y 
#
fit <- arima( Y, order=c(0,1,0), seasonal=list(order=c(1,0,0),period=4), include.mean=TRUE, method="ML" )
print(fit) 

# lets look at these residuals 
postscript("../../WriteUp/Graphics/Chapter6/prob_6_13_data_res_sacf_model_1.eps", onefile=FALSE, horizontal=FALSE)
acf(residuals(fit))
dev.off()

# based on the sample autocorrelation of the residuals above add an AR(1) term to the above model
# this is a SARIMA (0,1,1)(1,0,0)_{12} model on Y 
#
fit <- arima( Y, order=c(0,1,1), seasonal=list(order=c(1,0,0),period=4), include.mean=TRUE, method="ML" )
print(fit) 

# lets look at the new residuals 
postscript("../../WriteUp/Graphics/Chapter6/prob_6_13_data_res_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(residuals(fit))
dev.off()