# 
# Written by:
# -- 
# John L. Weatherwax                2009-02-24
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
#
# load( '../../Data/C8exercise1.mat' ); 
# plot( i, x, '-o' ); 
#
#-----

Y <- c( 1.1650, 0.6268, 0.3746, 0.5688, -0.0494, 1.3092, 2.0501, 4.1521, 5.3647, 6.6350, 5.5558, 3.9833, 4.1652, 3.4180, 3.4135, 2.8803, 2.3075, 0.5277, -1.7480, -1.9277,
       -1.9750, -2.7012, -3.7307, -3.2487, -2.8774, -2.0158, -1.0445, 0.2982, 1.4223, 0.6280, 0.5348, 1.5819, 2.8074, 0.9636, -0.3724, -1.4116, -1.5946, -2.1266, -2.3150,
       -0.6963, 1.0764, 2.2657, 1.8830, 1.5079, 1.3741, 0.6747, -0.0804, -0.0765, 1.4873, -0.1712, -2.1731, -2.0428, -1.2873, -1.0788, -0.4960, -1.1285, -2.4414, -2.7237,
       -1.2267, 0.4074, 0.7781, 1.8284, 2.0665, 1.3679, 0.4703, -0.6599, -1.4576, -1.2613, -2.1542, -2.1860, -0.9612, -0.0895, 1.0975, 1.6431, 3.8796, 5.8398, 4.8118, 4.0408,
        1.6068, 1.2160, 1.3348, 1.9729, 3.0994, 5.3164, 5.9498, 5.8633, 6.1384, 6.9606, 5.9495, 4.5459, 3.5864, 0.3966, -0.2190, -1.5729, -1.0854, -1.5495, -1.0924, 0.0340,
       0.4758, 1.2575 )

# begin by considering the autocorrelation function directly
# 
#postscript("../../WriteUp/Graphics/Chapter8/prob_1_direct_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(Y)
#dev.off()

# Plot the sample partial autcorrelation function ... note the signficant term at $k=1$
pacf( Y ) 

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

# we have elliminated any autocorrelation: 
pacf( fit$residuals )

# but we still don't have non-significant MA terms:
acf( fit$residuals )

# based on this lets assume a MA(1) term for the residuals or a ARIMA(1,0,1) model:
#
fit <- arima( Y, order=c(1,0,1), include.mean=FALSE, method="ML" )
print(fit)

# we have mostly elliminated any autocorrelation some small amount remains around lag 6: 
pacf( fit$residuals )

# but we have elliminated most of any MA components:
acf( fit$residuals )



# As an alternative model lets begin with the first difference to make the data stationary and model the first difference
# 
Yd <- diff(Y,lag=1,differences=1)
nd <- length(Yd)
plot(1:nd, Yd, type="l", xlab="time", ylab="first difference of given data")

# plot the autocorrelation function of this first difference.  There is a significant lag at k=1
#postscript("../../WriteUp/Graphics/Chapter8/prob_1_lag_one_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(Yd,lag.max=50)
#dev.off()

# this function is noisy with "no" large terms:
pacf(Yd,lag.max=50)

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

# lets look at these residuals and see if we have a sufficient model:
#postscript("../../WriteUp/Graphics/Chapter8/prob_1_data_res_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(residuals(fit))
pacf(residuals(fit))
#dev.off()