####################################################################

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

####################################################################

options(echo = FALSE)
library(fSeries)

#--------------------------------------------------------------------
# All plots will be named plot1.png, plot2.png... from each
# successive call to a plotting function below.

png(file="plots%d.png", width = 680, height = 450, pointsize = 12, bg="transparent", res = 200)

#--------------------------------------------------------------------
# Read in series data. Store in time series object "deaths". This data is
# from R and can be alternatively loaded using "data(deaths, package="MASS")".
# But use this general method so can load any series from file.

x <- scan("P:\\My Documents\\FrankMasci\\Filtering_tests_ARIMAsim\\ARIMA_in_R\\deaths_data.txt", list(0));
data <- x[[1]];

startyear <- 1974;
startperiod <- 1;  # will start at JAN.
freq <- 12;        # for monthly.
deaths <- ts( data, start=c(startyear, startperiod), frequency=freq );

# plot the timeseries:
plot(deaths);

#--------------------------------------------------------------------
# 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, 30);
pacf(deaths, 30);

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

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

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

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

# one year ahead forecasts:
write("\nForecasts:",file="");
predict(fit, 12);

#--------------------------------------------------------------------
# Simulation. Can only simulate non-seasonal components!
# can reduce a (p,d,q)(P,D,Q)_s model to a pure (p,d,q) model but
# algebra is difficult. Need general method.

var <- 0.001;
p <- 3;
d <- 1;
q <- 1;
ar_coefficients <- c(0.05,0.7,0.12);
ma_coefficients <- c(0.3);

ts.sim <- arima.sim(n = 500, list(order = c(p, d, q),
                    ar = ar_coefficients,  # ar coeffs.
                    ma = ma_coefficients), # ma coeffs.
                    sd = sqrt(var));

ts.plot(ts.sim)
acf(ts.sim, 36);
pacf(ts.sim, 36);

# other diagnostics for above:
write("\nAR and MA Roots:",file="");
armaRoots(ar_coefficients);
armaRoots(ma_coefficients);

write("\nTrue ACF of model:",file="");
armaTrueacf(list(ar = ar_coefficients, ma = ma_coefficients), lag.max = 30)

dev.off()
q()