# 
# Computes the maximum likelihood ARIMA coefficients for Series 5. Monthly Traffic Fatalities in Ontario, January 1960 to December 1W4
# 
# 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(
    61,    65,    55,    56,    91,    80,   135,   129,   129,   130,   109,   126,
    73,    68,    74,    95,   105,   108,   127,   108,   126,   154,   127,   103,
    95,    59,    68,    82,    92,   124,   139,   167,   138,   146,   128,   145,
    91,    66,    89,    98,   113,   130,   127,   157,   157,   136,   145,   112,
    71,    95,    95,   105,   116,   104,   128,   181,   130,   124,   123,   152,
    85,    86,   101,   105,   135,   142,   138,   190,   129,   175,   146,   179,
   103,   104,    80,   108,   125,   134,   150,   149,   148,   167,   160,   168,
    99,   102,    98,   104,   114,   165,   184,   161,   205,   189,   152,   146,
    72,    83,    97,   105,   134,   145,   172,   175,   149,   166,   142,   146,
   111,    84,   116,   102,   153,   176,   145,   175,   147,   174,   169,   131,
    81,    85,   107,   111,   135,   114,   150,   156,   146,   166,   154,   130,
    92,    80,    87,   101,   130,   151,   195,   224,   191,   188,   169,   161,
   121,    89,   105,   134,   155,   203,   192,   195,   225,   186,   172,   157,
   130,    84,   122,   116,   147,   183,   211,   256,   207,   196,   153,   154,
    94,    89,   118,   101,   150,   150,   191,   214,   173,   170,   175,   123
)
n <- length(Y)

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

# consider the autocorrelation function directly 
postscript("../../WriteUp/Graphics/Chapter6/prob_6_14_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 traffic fatalities") 

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

# based on the sample autocorrelation above lets assume an AR(1) and a SMA(1) model of Yd
#
fit <- arima( Y, order=c(1,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_14_data_d12_sacf_model_1.eps", onefile=FALSE, horizontal=FALSE)
acf(residuals(fit))
dev.off()

# based on the sample autocorrelation of these residuals above lets add in a MA(1) term to the original model
#
fit <- arima( Y, order=c(1,0,1), 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_14_data_res_sacf.eps", onefile=FALSE, horizontal=FALSE)
acf(residuals(fit))
dev.off()