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

# 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)