####################################################################
## Purpose:
## Simulate seasonal series with noise and a seasonal break in a
## multiplicative seasonal model.
##
## Execute following on dos command prompt (omit the "##"):
## "S:\\R\\rw2001\\bin\\Rterm.exe" -q --no-restore --no-save < simulate_SB.r 
## Unix:
## /usr/bin/R -q --no-restore --no-save < simulate_SB.r
####################################################################

options(echo = FALSE)

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

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

#--------------------------------------------------------------------
# Simulate seasonal seasonal series with SB in multiplicative 
# decomposition model:

# standard deviation of random Gaussian noise:
sd <- 0.00;

# Before break (20 years): 
N1 <- 20;
x <- 0.01*seq(1:(N1*12)) + (0.5*sin((pi/12)*seq(1:(N1*12)))**100) + rnorm(N1*12,sd=sd); 
x <- exp(x);
seas_fact_x <- x / exp(0.01*seq(1:(N1*12)));

# After break (10 years, factor decrease "f" in amplitude of seasonal component):
N2 <- 10;
f <- 0.3;
y <- log(x[[N1*12]]) + 0.01*seq(1:(N2*12)) + f*(0.5*sin((pi/12)*seq(1:(N2*12)))**100) + rnorm(N2*12,sd=sd); 
y <- exp(y);
seas_fact_y <- y / exp(log(x[[N1*12]]) + 0.01*seq(1:(N2*12)));

z <- c(x,y);
z_seas <- c(seas_fact_x, seas_fact_y);

#--------------------------------------------------------------------
# Store into ts objects and plot.

startyear <- 1970;
startperiod <- 1;
freq <- 12;
zts <- ts( z, start=c(startyear, startperiod), frequency=freq ); 
zseasts <- ts( z_seas, start=c(startyear, startperiod), frequency=freq );

ts.plot(zts);
ts.plot(zseasts);

years <- time(zts);
mnths <- rep(seq(1:freq),(N1+N2));

#--------------------------------------------------------------------
# Write out series data to file in SEASABS format.

for(i in 1:((N1+N2)*freq)) {
   write(paste(trunc(years[[i]]), mnths[[i]], zts[[i]]), file="Simseries_withSB_noiseless.txt", append=T);
}