####################################################################
##Execute following on dos command prompt (omit the "##"):
##"S:\\R\\rw2001\\bin\\Rterm.exe" -q --no-restore --no-save < ARIMA_fit.r
##Unix:
##/usr/bin/R -q --no-restore --no-save < ARIMA_fit.r
####################################################################

options(echo = FALSE)

#--------------------------------------------------------------------
# All plots will be named plot001.png, plot002.png... from each

# successive call to a plotting function below.

#bitmap(file="Noreg_plot_%03d.png", width=5, height=5, res=400)
png(file="Noreg_plot_%03d.png", width = 680, height = 450, pointsize = 12, bg="transparent", res = 200)

#--------------------------------------------------------------------
# Read in series data. Take log of data like in X12-ARIMA. 

x <- scan("new_nsw_largecorr.dat", list(0,0,0));
data <- log(x[[3]]);

startyear <- 1982;
startperiod <- 4;
freq <- 12;
deaths <- ts( data, start=c(startyear, startperiod), frequency=freq );

# plot the timeseries:
plot(deaths);

#--------------------------------------------------------------------
# Read in regressor. Vector of same dimension as series.
# From "R" help on "arima":
# When regressors are specified, they are orthogonalized prior to 
# fitting unless any of the coefficients is fixed. It can be helpful 
# to roughly scale the regressors to zero mean and unit variance.

r <- scan("nswld.txt", list(0,0,0));
nswid <- r[[3]];

#--------------------------------------------------------------------
# Fit ARIMA(p,d,q)(P,D,Q), epsilon ~ N(0,var) model to series.

# plot acf/pacf of raw series to see if differencing needed:
acf(deaths, 36);
pacf(deaths, 36);

# seasonally difference the series: Y(t) - Y(t-12)  
deaths.diff <- diff(deaths, 12)
deaths.diff <- diff(deaths.diff, 1)
plot(deaths.diff);

# plot acf/pacf of differenced series:
acf(deaths.diff, 36);
pacf(deaths.diff, 36);

# Fit arima model:
(fit <- arima(deaths,                         # must be univariate TS.
             order = c(0,1,1),                # non-seasonal
             seasonal = list(order = c(0,1,1),
             period = 12),                    # seasonal
             xreg = NULL,                     # regressors!
             include.mean = FALSE) )          # FALSE if mean term to be removed by differencing.

# diagnostic plots:
tsdiag(fit, gof.lag = 36)

dev.off()
q()