# 
# Computes the maximum likelihood SARIMA coefficients for the quarterly GM stock earnings per share
# 
# 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:
X <- c( 1951:1979 ) 
Y <- c(
    0.5300,    0.5200,    0.3400,    0.5000,
    0.4700,    0.5300,    0.4400,    0.6400,
    0.5700,    0.6000,    0.5200,    0.5400,
    0.7100,    0.8900,    0.6000,    0.8300,
    1.1400,    1.2700,    0.9000,    0.9900,
    1.0100,    0.7900,    0.4800,    0.7400,
    0.9300,    0.7800,    0.4300,    0.8500,
    0.6500,    0.5200,    0.2200,    0.8300,
    1.0300,    1.0500,    0.4700,    0.5100,
    1.1400,    1.0100,    0.3000,    0.9000,
    0.6500,    0.8800,    0.3000,    1.2000,
    1.2300,    1.4100,    0.6400,    1.5500,
    1.4400,    1.6200,    0.7200,    1.7700,
    1.8700,    2.1000,    0.7700,    1.3000,
    2.2200,    2.2300,    0.9100,    2.0500,
    2.0800,    1.9100,    0.3500,    1.9000,
    1.3500,    1.8200,    0.5100,    1.9700,
    1.4600,    1.8800,    0.6200,    2.0500,
    1.8200,    1.5600,    0.7900,    1.7700,
    1.2100,    1.6400,   -0.2800,   -0.4900,
    2.1200,    1.9700,    0.7500,    1.8800,
    2.2600,    2.5200,    0.4100,    2.3200,
    2.8400,    2.7800,    0.9200,    1.8000,
    0.4100,    1.0500,    0.0500,    1.7600,
    0.2000,    1.1400,    6.8400,    2.1400,
    2.7800,    3.1600,    1.3700,    2.7700,
    3.1400,    3.8200,    1.4000,    3.2600,
    3.0300,    3.8300,    1.8400,    3.5400,
    4.3900,    4.1300,    0.0600,    1.4600
) 
n <- length(Y)

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

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

# take the first seasonal difference to make the data stationary (this looks like it has a period of four) 
Yd <- diff(Y,lag=4,differences=1)
nd <- length(Yd)
plot(1:nd, Yd, type="l", xlab="time", ylab="first seasonal difference of quarterly GM stock earnings") 

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

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

# an alternative arima call in which arima does all of the seasonal differences itself (which seems to give more reasonable values):
#
fit <- arima( Y, order=c(0,0,0), seasonal=list(order=c(0,1,1),period=4), include.mean=TRUE, method="ML" )
print(fit) 

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

# lets predict the next 4 observations:
predict(fit, n.ahead=4 )