################################################################### ################# Generate synthetic data ######################## ################################################################### # This script can be used to build synthetic data under controlled # conditions for test purposes. #### define working environment #### #setwd("D:/ETC-ACM/Tasks/1022_3/R/finalisedScripts") ######## Function to create time series with three components ######### buildTS <- function(length=1000, meanT=50, slopeT=-0.001, facSin1=0.15, facSin2=5, facSinTot=10, sdNoise=10){ # linear trend data <- data.frame(t=meanT + (1:length)*slopeT) # seasonal pattern data$s<- facSinTot*(sin(facSin1*(1:length))+sin(facSin2*(1:length))) # random noise #set.seed(1100) data$r <- rnorm(length, mean=0, sd=sdNoise) # add components data$rts <- data$r + data$t + data$s return(data) } ################################################################### ######################### I Outliers ########################### ################################################################### I <- buildTS(length=10000, slopeT=0) # plot components par(mfrow=c(4,1)) par(mar=c(2,2,2,2)) ts.plot(I$t, main="Linear trend") ts.plot(I$s, main="Seasonal pattern") ts.plot(I$r, main="Random noise") ts.plot(I$rts, main="Combined time series") # chooses outlier locations I.pos <- round(runif(15,1,10000)) I.neg <- round(runif(15,1,10000)) # choose values outside the 95 % CI (>3*sd) I$r[I.pos]<- 60 I$r[I.neg]<- -60 I$rt <- I$t + I$r I$rs <- I$r + I$s I$rts <- I$r + I$t + I$s # plot outliers par(mfrow=c(1,1)) ts.plot(I$rts, main="Additive outliers") points(I.pos,I$rts[I.pos],col="red") points(I.neg,I$rts[I.neg],col="red") # Save time series save(I, I.pos, I.neg, file="I_outliers.RData") ################################################################### #################### II Structure changes ###################### ################################################################### ######## a) Permanent Change ######### IIa <- buildTS(length=10000, slopeT=0) # plot components par(mfrow=c(4,1)) par(mar=c(2,2,2,2)) ts.plot(IIa$t, main="Linear trend") ts.plot(IIa$s, main="Seasonal pattern") ts.plot(IIa$r, main="Random noise") ts.plot(IIa$rts, main="Combined time series") # Add breaks to trend (size 0.1-0.5 sdNoise) IIa.pos <- round(runif(1,1,10000)) IIa.min <- round(runif(1,1,10000)) # large break IIa$t[IIa.pos:nrow(IIa)] <- IIa$t[IIa.pos:nrow(IIa)] + 20 # change in minimum (detection limit) IIa$rts <- IIa$r + IIa$t + IIa$s IIa$rts[which(IIa$rts<30)[which(IIa$rts<30)