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

options(echo = FALSE)

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

#png(file="plot%03d.png", width = 680, height = 450, pointsize = 12, bg="transparent", res = 200)
jpeg(filename="plot%03d.jpg", width=700, height=480, pointsize=12, quality=75, bg="white", res=NA)
#bitmap(file="plot%03d.png", width=5, height=5, res=400)

#--------------------------------------------------------------------
# Read in series data.

x1 <- scan("RetailPG98_D8_WithTDcorr.txt", list(0,0,0));
xd1 <- x1[[3]];
x2 <- scan("RetailPG98_D8_NoTDcorr.txt", list(0,0,0));
xd2 <- x2[[3]];

startyear <- 1962;
startperiod <- 4;
# All data:
ts1 <- ts( xd1, start=c(startyear, startperiod), frequency=12 );
ts2 <- ts( xd2, start=c(startyear, startperiod), frequency=12 );

# plot full timeseries:
ts.plot(ts1, ts2, gpars=list(ylab="With TD (BLACK); Without TD (RED)", col=c(1:2)),
        main="S*I (D8) with/without static TD");

#--------------------------------------------------------------------
# Spectra.

# frequency scale is fraction of a cycle per Delta_t sampling interval
# in input series. Max=0.5 cycle per Delta_t => 1 cycle=2*Delta_t=2 Mnths:
# Period (time represented by 1 complete cycle) = "1/f" mnths. (= 2*pi/omega where
# omega = 2*pi*f).

spectrum_xd <- spectrum(cbind(xd1,xd2), method="pgram", plot=TRUE);
# monthly TD freqs:
abline(v=(1/2.87), col="green", lty=1, lwd=1);
abline(v=(1/2.32), col="green", lty=1, lwd=1);
abline(v=(1/3.29), col="green", lty=1, lwd=1);
abline(v=(1/4.67), col="green", lty=1, lwd=1);
abline(v=(1/12), col="blue", lty=1, lwd=1);

spectrum_xd1 <- spectrum(xd1, method="pgram", plot=FALSE);
spectrum_xd2 <- spectrum(xd2, method="pgram", plot=FALSE);

f <- 1/(spectrum_xd1$freq);
s1 <- spectrum_xd1$spec;
s2 <- spectrum_xd2$spec;

#--------------------------------------------------------------------
# Customized plot with two spectra and legend:

plot(f, s1, main="Spectra of S*I (D8): \"Retail (PG98)\"", xlab="Period (Months)",
     ylab="Spectral Density", pch=" ", frame=TRUE, xlim=range(2,6.5),
     ylim=range(0,0.04));
lines(f, s1, col="black",lty=1)
lines(f, s2, col="red",lty=1)
abline(v=4.67, col="green", lty=1, lwd=1);
abline(v=3.29, col="green", lty=1, lwd=1);
abline(v=2.32, col="green", lty=1, lwd=1);
abline(v=2.87, col="green", lty=1, lwd=1);
abline(v=12, col="blue", lty=1, lwd=1);
abline(v=6, col="blue", lty=1, lwd=1);
abline(v=4,col="blue", lty=1, lwd=1);
abline(v=3, col="blue", lty=1, lwd=1);
abline(v=2.4, col="blue", lty=1, lwd=1);

legend(4.7, 0.04, c("TD corrected",
                  "No TD correction",
                  "QTD cycles",
                  "Seasonal cycles"),
           col = c("black",
                   "red",
                   "green",
                   "blue"),
           text.col="black",
           lty = c(1, 1, 1, 1),
           lwd = c(1, 1, 1, 1),
           pch = c(" ", " ", " ", " "),
           merge=FALSE, bty='n');

dev.off()
q()