# 
# Computes the maximum likelihood ARIMA coefficients for the car sales 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(
    6.550,  8.728, 12.026, 14.395, 14.587, 13.791, 9.498,   8.251,  7.049,  9.545,  9.364,  8.456 ,
    7.237,  9.374, 11.837, 13.784, 15.926, 13.821, 11.143,  7.975,  7.610, 10.015, 12.759,  8.816 ,
   10.677, 10.947, 15.200, 17.010, 20.900, 16.205, 12.143,  8.997,  5.568, 11.474, 12.256, 10.583 ,
   10.862, 10.965, 14.405, 20.379, 20.128, 17.816, 12.268,  8.642,  7.962, 13.932, 15.936, 12.628 ,
   12.267, 12.470, 18.944, 21.015, 22.015, 18.581, 15.175, 10.306, 10.792, 14.752, 13.754, 11.738 ,
   12.181, 12.965, 19.990, 23.125, 23.541, 21.247, 15.189, 14.767, 10.895, 17.130, 17.697, 16.611 ,
   12.674, 12.760, 20.249, 22.135, 20.677, 19.933, 15.388, 15.113, 13.401, 16.135, 17.562, 14.720 ,
   12.225, 11.608, 20.985, 19.692, 24.081, 22.114, 14.220, 13.434, 13.598, 17.187, 16.119, 13.713 
) 
Y_hold_out <- c( 13.210, 14.251, 20.139, 21.725, 26.099, 21.084, 18.024, 16.722, 14.385, 21.342, 17.180, 14.577 )
n <- length(Y)

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

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

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

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

# based on the sample autocorrelation above lets try a MA(2) model and a SMA(1) for Yd 
# 
fit <- arima( Yd, order=c(0,0,2), seasonal=list(order=c(0,0,1),period=12), include.mean=TRUE, method="ML" )
print(fit)

# or based on the sample autocorrelation above lets try a AR(2) model and a SMA(1) for Yd
# 
fit <- arima( Yd, order=c(2,0,0), seasonal=list(order=c(0,0,1),period=12), include.mean=TRUE, method="ML" )
print(fit)

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

# to compute the forecasts we'll let R do the differences internally in the arima function:
#
fit <- arima( Y, order=c(2,0,0), seasonal=list(order=c(0,1,1),period=12), include.mean=TRUE, method="ML" )
print(fit)

predict(fit, n.ahead=length(Y_hold_out) )