# 
# Computes the maximum likelihood ARIMA coefficients for the unemployed U.S. civilian labor force (in thousands) for the period 1970-1978 
# 
# 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(
        3406,        3794,        3733,        3552,        3384,        4669,
        4510,        4220,        4292,        4259,        4607,        4636,
        5414,        5442,        5175,        4694,        4394,        5490,
        5330,        5061,        4840,        4570,        4815,        4695,
        5447,        5412,        5215,        4697,        4344,        5426,
        5173,        4857,        4658,        4470,        4266,        4116,
        4675,        4845,        4512,        4174,        3799,        4847,
        4550,        4208,        4165,        3763,        4056,        4058,
        5008,        5140,        4755,        4301,        4144,        5380,
        5260,        4885,        5202,        5044,        5685,        6106,
        8180,        8309,        8359,        7820,        7623,        8569,
        8209,        7696,        7522,        7244,        7231,        7195,
        8174,        8033,        7525,        6890,        6304,        7655,
        7577,        7322,        7026,        6833,        7095,        7022,
        7848,        8109,        7556,        6568,        6151,        7453,
        6941,        6757,        6437,        6221,        6346,        5880,
        6897,        6735,        6479,        5685,        5457,        6326,
        6438,        5931,        5797,        5460,        5629,        5725
) 
n <- length(Y)

# we plot the original data:
postscript("../../WriteUp/Graphics/Chapter6/prob_6_11_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_11_data_direct_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(Y)
dev.off()

# take a 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 monthly arrivals of U.S. citizens") 

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

# take a lag=12 difference to make this data stationary:
Ydd <- diff(Yd,lag=12,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_11_data_dds_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(Ydd,lag.max=50)
dev.off()

# based on this we try to add an AR(2) term:
#
fit <- arima( Y, order=c(2,1,0), seasonal=list(order=c(0,1,0),period=12), include.mean=TRUE, method="ML" )
print(fit) 

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

# based on this looks like we are now missing a SMA(1) term
#
fit <- arima( Y, order=c(2,1,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_11_data_res_sacf_final.eps", onefile=FALSE, horizontal=FALSE)
acf(residuals(fit))
dev.off()