# 
# Computes the maximum likelihood SARIMA coefficients for the montly arrival of U.S. citizens
# 
# 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(
 550,  444, 517, 563, 573, 594,
 897, 1065, 769, 647, 544, 427,
 655,  579, 618, 765, 704, 749,
1055, 1130, 844, 771, 664, 543,
 663,  589, 713, 780, 775, 790,
 993, 1171, 761, 751, 630, 594,
 620,  601, 720, 767, 706, 726, 
 906, 1054, 753, 599, 571, 518,
 627,  531, 553, 624, 625, 701,
 872, 1003, 653, 658, 606, 514,
 571,  493, 585, 590, 617, 711,
 825,  936, 683, 687, 535, 468,
 588,  511, 618, 645, 643, 710,
 919, 1002, 719, 760, 575, 511,
 633, 1077, 742, 740, 612, 584 
) 
n <- length(Y)

# we plot the original data:
postscript("../../WriteUp/Graphics/Chapter6/prob_6_10_data_plt.eps", onefile=FALSE, horizontal=FALSE)
plot(1:n, Y, type="l", xlab="time", ylab="monthly arrivals of U.S. citizens") 
dev.off()

# consider the autocorrelation function directly 
postscript("../../WriteUp/Graphics/Chapter6/prob_6_10_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 monthly arrivals of U.S. citizens") 

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

# take a lag=1 difference to make this data stationary:
Ydd <- diff(Yd,lag=1,differences=1)
nd  <- length(Ydd)
plot(1:nd, Ydd, type="l", xlab="time", ylab="first difference of first seasonal difference of monthly arrivals of U.S. citizens") 

# plot the autocorrelation function of this combined set of differences:
postscript("../../WriteUp/Graphics/Chapter6/prob_6_10_data_dds_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(Ydd)
dev.off()

# based on these sample autocorrelation above lets assume a MA(1) model of Ydd
# this is a SARIMA (0,0,1)(0,0,0)_{12} model on Ydd
#      or a SARIMA (0,1,1)(0,1,0)_{12} model on Y 
#

fit <- arima( Ydd, order=c(0,0,1), seasonal=list(order=c(0,0,0),period=12), include.mean=TRUE, method="ML" )      # <- see if \theta_0 is significant ... NO
print(fit) 

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

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

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