# 
# Computes the maximum likelihood ARIMA coefficients for the growth rate 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(
      601, 604, 620, 626, 641, 642, 645, 655, 
      682, 678, 692, 707, 736, 753, 763, 775, 
      775, 783, 794, 813, 823, 826, 829, 831, 
      830, 838, 854, 872, 882, 903, 919, 937,
      927, 962, 975, 995, 1001, 1013, 1021, 1028,
      1027, 1048, 1070, 1095, 1113, 1143, 1154, 1173,
      1178, 1183, 1205, 1208, 1209, 1223, 1238, 1245,
      1258, 1278, 1294, 1314, 1323, 1336, 1355, 1377,
      1416, 1430, 1455, 1480, 1514, 1545, 1589, 1634,
      1669, 1715, 1760, 1812, 1809, 1828, 1871, 1892,
      1946, 1983, 2013, 2045, 2048, 2097, 2140, 2171,
      2208, 2272, 2311, 2349, 2362, 2442, 2479, 2528,
      2571, 2634, 2684, 2790, 2890, 2964, 3085, 3159,
      3237, 3358, 3489, 3588, 3624, 3719, 3821, 3934,
      4028, 4129, 4205, 4349, 4463, 4598, 4725, 4827,
      4939, 5067, 5231, 5408, 5492, 5653, 5828, 5965
)
n <- length(Y)

# compute the growth rate from the raw data
X <- c()
for( ii in 1:n-1 ){
  X[ii] <- (Y[ii+1]-Y[ii])/Y[ii]
}

# take the first difference to make the data stationary
Xd <- diff(X)

# construct a ARIMA(0,1,1) model on this data
#       or an ARIMA(0,0,1) model on the first difference of this data
fit <- arima( Xd, order=c(0,0,1), include.mean=TRUE, method="ML" )

print(fit)

# restimate without including the mean:
fit <- arima( Xd, order=c(0,0,1), include.mean=FALSE, method="ML" )

print(fit)

###predict( fit, n.ahead=1 )