# 
# Computes the maximum likelihood ARIMA coefficients for the gas usage data set.
# 
# 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: 
Y <- c(
    302, 262, 218, 175, 100,  77, 43, 47, 49, 69, 152, 205,
    246, 294, 242, 181, 107, 56, 49, 47, 47, 71, 151, 244, 
    280, 230, 185, 148, 98, 61, 46, 45, 55, 48, 115, 185, 
    276, 220, 181, 151, 83, 55, 49, 42, 46, 74, 103, 200, 
    237, 247, 215, 182, 80, 46, 65, 40, 44, 63, 85, 185, 
    247, 231, 167, 117, 79, 45, 40, 38, 41, 69, 152, 232, 
    282, 255, 161, 107, 53, 40, 39, 34, 35, 56, 97, 210, 
    260, 257, 210, 125, 80, 42, 35, 31, 32, 50, 92, 189
) 
#    256, 250, 198, 136, 73, 39, 32, 30, 31, 45 <- excluded data ... held out 
n <- length(Y)

# take the first seasonal difference to make the data stationary
Yd <- diff(Y,lag=12,differences=1)

# construct a SARIMA (0,0,1)(0,1,1)_{12} model on Y
#        or a SARIMA (0,0,1)(0,0,1)_{12} model on Yd ... this call matches the results from the book quite nicely ... 
#
fit <- arima( Yd, order=c(0,0,1), seasonal=list(order=c(0,0,1),period=12), include.mean=TRUE, method="ML" )

print(fit)

predict( fit, n.ahead=10 ) 

# Note: the call below does not match the book as well as the previous call even though they should be equivalant.
#
fit <- arima( Y, order=c(0,0,1), seasonal=list(order=c(0,1,1),period=12), include.mean=TRUE, method="ML" )

print(fit)

predict( fit, n.ahead=10 ) 

# lets look at the residuals 
#postscript("test.eps", onefile=FALSE, horizontal=FALSE)
#acf(residuals(fit))
#dev.off()
#pacf(residuals(fit))
# diffinv(x,lag=,difference=)