# 
# Computes the maximum likelihood ARIMA coefficients for the demand for parts 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(
    954,  765,  867,  940,  913, 1014,  801,  990,  712,  959,  828,  965 ,
    915,  891,  991,  971, 1129, 1091, 1195, 1295, 1046, 1121, 1033, 1222 ,
   1199, 1012, 1404, 1137, 1421, 1162, 1639, 1545, 1420, 1916, 1491, 1295 ,
   1764, 1727, 1654, 1811, 1520, 1635, 1984, 1898, 1853, 2015, 1709, 1667 ,
   1625, 1562, 2121, 1783, 1474, 1657, 1746, 1763, 1517, 1457, 1388, 1501 ,
   1227, 1342, 1666, 2091, 1629, 1451, 1727, 1736, 1952, 1420, 1345,  842 ,
   1576, 1485, 1928, 2072, 1887, 2175, 2199, 1961, 2091, 1993, 1595, 1372 
) 
n <- length(Y)

# take the log transformation to stabalize the variance: 
Y <- log(Y) 

# we plot the original data:
postscript("../../WriteUp/Graphics/Chapter6/demand_data_plt.eps", onefile=FALSE, horizontal=FALSE)
plot(1:n, Y, type="l", xlab="time", ylab="log of demand for repair parts") 
dev.off()

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

# take the first 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 demand for repair parts") 

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

# based on the sample autocorrelation above lets assume a MA(1) model of Yd
#
fit <- arima( Yd, order=c(0,0,1), include.mean=TRUE, method="ML" )
print(fit)

# drop the insignificant mean component and refitting 
#
fit <- arima( Yd, order=c(0,0,1), include.mean=FALSE, method="ML" )
print(fit)

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

# these residuals indicate we might need a SMA(1) component to be added into the model
#
fit <- arima( Yd, order=c(0,0,1), seasonal=list(order=c(0,0,1),period=12), include.mean=TRUE, method="ML" )
print(fit) 

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