############################### # Chapter 2 ############################### # Section 2.1 # Non linear function # Smooth transition trend1 = function(beta0,beta1,n){ x = c(1:n) Trend = (1+exp(beta0*(t-beta1)))**(-1) y = Trend + rnorm(n,sd =0.3) return(y) } Ttrend1 = function(beta0,beta1,n){ x = c(1:n) Trend = (1+exp(beta0*(t-beta1)))**(-1) y = Trend + rnorm(n,sd =0.0) return(y) } example1 = trend1(0.2,60,100) true1 = Ttrend1(0.2,60,100) example2 = trend1(6,60,100) true2 = Ttrend1(6,60,100) par(mfrow=c(1,2)) plot.ts(example1,col="blue",lwd=2) lines(true1,lty=2) plot.ts(example2,col="red",lwd=2) lines(true2,lty=2) ############################################# # Section 2.1 # Non linear function # Smooth transition # Burst function burst = function(A,b,alpha,theta,n){ temp = rep(0,n) for(t in c(1:n)){ temp[t] = A*exp(b*(1-cos(alpha*t)))*cos(theta*t) } return(temp) } signal1 = burst(1,b=2, alpha =1, theta = pi/2,n=100) Ysignal = signal1 + rnorm(n=100, sd = 8) par(mfrow=c(1,2)) plot.ts(signal1,lwd=2,col="blue") plot.ts(Ysignal,lwd=2,col="red") lines(signal1,lty=2) ################################ # Section 2.2 (Taking differences) # Simulating a random walk. Taking first differences will remove the random trend. # Create Random Walk (of length n) and plot n = 200 # make vector of length n, with zeros temp = rep(0,n) # initialize with a drawn from a standard random normal temp[1] = rnorm(1) for(i in c(2:n)){ temp[i] = (temp[(i-1)]+rnorm(1)) } # make time series plot of realisation plot.ts(temp) # Section 2.3.1 # See the Rcode ######################################### # Section 2.5 # periodic function (no noise) oscill = rep(c(rep(1,3),rep(0,2)),20) n=length(oscill) FO = fft(oscill)/n freq = 2*pi*c(1:n)/n par(mfrow=c(1,2)) plot(freq,Re(FO),lwd=2,col="blue",type="l") plot(freq,Im(FO),lwd=2,col="blue",type="l") par(mfrow = c(1,1)) plot.ts(oscill,lwd=2,col="blue",type="l") par(mfrow=c(2,1)) plot.ts(Re(per)) plot.ts(Im(per)) #################################################### # Add noise to system # Signal oscillM = rep(c(rep(1,1),rep(-0.5,3)),25) # centered Signal oscillM = oscillM - mean(oscillM) # Signal plus noise oscill = oscillM + rnorm(100,sd=0.5) ############## n=length(oscill) FO = fft(oscill)/n freq = 2*pi*c(1:n)/n t = c(1:n) ############ # plot of observations par(mfrow=c(1,2)) plot(t,oscill,col="blue",type="l",lwd=2) lines(t,oscillM,col="black",type="l",lwd=2,lty=2) plot(t,oscill,col="blue",type="l",lwd=2) lines(t,sin(2*freq),col="green",lwd=1) lines(t,sin(25*freq),col="red",lwd=1) FO = fft(oscill)/n F00 = fft(oscillM)/n freq = 2*pi*c(1:n)/n ############## # plot of cosine and sine transforms par(mfrow=c(1,2)) plot(freq,Re(FO),lwd=2,col="blue",type="l",ylab="Cosine") lines(freq,Re(F00),lwd=1,col="red",type="l") plot(freq,Im(FO),lwd=2,col="blue",type="l",ylab="Sine") lines(freq,Im(F00),lwd=1,col="red",type="l") ############################## # Section 2.5.5 # Smooth trend function and its fft n1=100 rescaled = c(1:n1)/n1 lambda1 = (6*(-rescaled**2+rescaled)-0.7) y = lambda1 + rnorm(100,sd=0.5) time =c(1:100) freq = 2*pi*c(0:19)/n1 FO = abs(fft(y))**2/n1 F1 = FO[c(1:20)] F2 = abs(fft(lambda1))**2/n1 F3 = F2[c(1:20)] par(mfrow=c(1,2)) plot(time,y,lwd=2,type="l",col="blue",main="Signal plus noise") lines(time,lambda1,lwd=2,type="l",col="red") plot(freq,F1,lwd=2,type="l",col="blue",main="Periodogram") points(freq,F1,lwd=2,col="blue") lines(freq,F3,lwd=2,col="red",lty=2) points(freq,F3,lwd=2,col="red") ################################### # Section 2.5.6 Period detection # generating sin + iid n = 128 # Period P = 8 # Generate periodic signal plus iid noise temp <- rnorm(n) sin1 <- 1.5*sin(2*pi*c(1:n)/P) signal <- 1.5*sin(2*pi*c(1:n)/P) + temp PS <- abs(fft(signal)/n)**2 P1 <- abs(fft(temp)/n)**2 P2 <- abs(fft(sin1)/n)**2 PS <- PS[c(1:n/2)] P1 <- P1[c(1:n/2)] P2 <- P2[c(1:n/2)] frequency = 2*pi*c(0:(n-1))/n # Make plot plot(frequency,PS,type = "l") # Locate maximum in R which.max(PS) ################################################### # Generate periodic signal plus dependent noise # generate an AR(2) with Gaussian innovations temp <- arima.sim(list(order=c(2,0,0), ar = c(1.4*cos(2*pi/5), -(0.7**2))), n) sin1 <- 1.5*sin(2*pi*c(1:n)/P) signal <- 1.5*sin(2*pi*c(1:n)/P) + temp PS <- abs(fft(signal)/n)**2 P1 <- abs(fft(temp)/n)**2 P2 <- abs(fft(sin1)/n)**2 PS <- PS[c(1:n/2)] P1 <- P1[c(1:n/2)] P2 <- P2[c(1:n/2)] frequency = 2*pi*c(0:(n-1))/n # Make plot plot(frequency,PS,type = "l") # Locate maximum in R which.max(PS)