# 
# Computes the maximum likelihood SARIMA coefficients for the monthly percentage changes in Canadian wages 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(
     0.71, 1.47, 1.61, 2.97, 3.58, 0.94,  0.68,  0.32, -0.67,  0.02, -2.56, 
    -0.84, 0.38, 1.18, 1.95, 4.08, 2.41,  0.55,  1.16,  1.70,  0.13,  0.63, -2.19 ,
    -0.32, 1.48, 0.84, 1.45, 3.00, 3.24, -0.62,  0.35,  3.26, -0.53,  0.90, -1.31 ,
    -0.91, 0.46, 0.50, 0.96, 2.23, 3.27, -1.95,  0.65,  3.05,  0.23, -0.45, -1.78 ,
    -1.44, 1.49, 1.22, 2.22, 4.76, 3.04, -2.39,  0.89,  3.33, -0.14, -0.47, -1.79 ,
    -0.49, 0.58, 2.11, 0.69, 3.81, 3.32, -2.26,  1.46,  3.81,  0.96, -0.03,  0.50 ,
    -1.70, 2.51, 0.75, 2.21, 3.15, 3.05, -2.36,  0.26,  4.94,  1.44, -0.12, -1.08 ,
     0.19, 1.03, 2.12, 1.58, 4.51, 2.96,  0.30,  1.66,  4.65, -0.43, -0.43,  1.19 ,
    -2.93, 0.21, 2.09, 1.60, 4.64, 2.96,  0.86, -1.88,  5.01 
)
n <- length(Y)

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

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

# take the first seasonal difference:
Yd <- diff(Y,lag=12,differences=1) 
nd <- length(Yd)
plot(1:nd, Yd, type="l", xlab="time", ylab="first seasonal difference of monthly arrivals of U.S. citizens") 

# plot the autocorrelation function of this first seasonal difference 
postscript("../../WriteUp/Graphics/Chapter6/prob_6_12_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)_{12} model on Y 
#

fit <- arima( Yd, order=c(0,0,0), seasonal=list(order=c(0,0,1),period=12), include.mean=TRUE, method="ML" ) # <- is the mean \theta_0 significant?  No.
print(fit) 

fit <- arima( Y, order=c(0,0,0), seasonal=list(order=c(0,1,1),period=12), include.mean=TRUE, method="ML" )
print(fit) 

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

predict( fit, n.ahead=12 )