# 4.12.
library(ts)

#Diskrete Fourier-Tranformation: Wir betrachten die Zerlegung einer
#Reihe x(0), x(1), ... , x(N-1) in harmonische Schwingungen mit den
#Frequenzen 2*pi*j/N, j=0,1,2,...,[N/2].
f.fourier <- function(x)
{
  ## Purpose: Discrete Fourier transform given as amplitude and phase.
  ## The fast Fourier transform (FFT) is used
  ## -------------------------------------------------------------------------
  ## Arguments: x (real) series
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 29 Nov 2000, 14:37
  N <- length(x)
  a <- fft(x)
  freq <- (0:(N/2))*2*pi/N
  a <- a[1:length(freq)]
  list(freq=freq,ampl=Mod(a),phase=Arg(a))
}

#Befehle zur Analyse einer Reihe x
xf <- f.fourier(x)
par(mfrow=c(2,1))
plot(xf$freq,xf$ampl,type="h",xlab="Frequenz", ylab="Amplitude")
abline(h=0)
plot(xf$freq,xf$phase,type="h",xlab="Frequenz", ylab="Phase")
abline(h=0)

#Verschiedene Wahlen von x:
#1.Harmonische Schwingung mit Frequenz verschieden von den Fourierfrequenzen
x <- cos(2*pi*40.5*(0:127)/128)
#Das Verschmieren einer Nicht-Fourierfrequenz kann durch ein sogenanntes
#taper reduziert werden. Ein taper multipliziert die Beobachtungen am
#Anfang und am Ende mit Gewichten, die monoton gegen null gehen.
x <- spec.taper(x)

#2.Zwei harmonische Schwingungen mit fast gleicher Frequenz
x <- cos(2*pi*80.5*(0:255)/256)+0.5*cos(2*pi*81.5*(0:255)/256) 


#3.Weisses Rauschen
x <- rnorm(256)
#Zusaetzlich: Betrachte Histogramme von Amplitude und Phase
hist(xf$ampl,nclass=12)
hist(xf$phase,breaks=seq(-pi,pi,length=12))

#4.Harmonische Schwingung ueberlagert von weissem Rauschen
x <- cos(2*pi*80.2*(0:255)/256) + rnorm(256)
#Nehmen Sie stattdessen auch eine Fourierfrequenz oder zwei Frequenzen

#5.Luchsdaten
data("lynx")
x <- log(lynx)
x <- log(lynx)-mean(log(lynx))
#Periode der am staerksten vertretenen Schwingung
k <- length(xf$freq)
length(x)/(0:(k-1))[xf$ampl == max(xf$ampl)]
plot(log(lynx))

#6.Ozeanwellen
x <- scan.url("http://stat.ethz.ch/Teaching/Datasets/ocwave.dat")
x <- x-mean(x)
x <- spec.taper(x)
#Vergleich von erster und zweiter Haelfte der Daten
x1  <- spec.taper(x[1:512])
xf <- f.fourier(x1)
plot(xf$freq,xf$ampl,type="h",xlab="Frequenz", ylab="Amplitude")
x2  <- spec.taper(x[513:1024])
xf  <- f.fourier(x2)
plot(xf$freq,xf$ampl,type="h",xlab="Frequenz", ylab="Amplitude")

#7.Jaehrliche Minima des Nils von 622 - 1284
x <- ts(scan.url("http://stat.ethz.ch/Teaching/Datasets/nilmin.dat"),start=622)
plot(x)
x <- x-mean(x)
#Darstellung in anderer Skala:
xf <- f.fourier(x)
par(mfrow=c(1,1))
plot(xf$freq[-1],xf$ampl[-1],xlab="log-Frequenz",
     ylab="log-Amplitude",log="xy")

#8. Gleitendes Mittel: Vergleich von Amplitude^2 mit dem wahren Spektrum
f.transfer <- function(coef,lambda)
{
  ## Purpose: Berechnung von \sum coeff(u)*exp(-i*lambda*u) 
  ##          (sogenannte Transferfunktion)
  ## -------------------------------------------------------------------------
  ## Arguments: coef=Koeffizienten, lambda=Frequenz
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  2 Dec 96, 11:11
  p <- length(coef)-1
  z <- (cos(lambda)+sin(lambda)*1i)
  a  <- outer(z,0:p,FUN="^")
  a%*%coef
}
lambda <- seq(0,pi,0.01)
coef <- cos(0.3*pi*(0:7)) #Festlegung der Koeffizienten
coef <- coef/sqrt(sum(coef^2)) #Normierung der Koeffiziente
tr <- f.transfer(coef,lambda)
par(mfrow=c(2,1))
plot(lambda,Mod(tr)^2,type="l")
x <- ts(rnorm(520),start=1,freq=1)
x <- filter(x,filter=coef,sides=1)[9:520]
xf <- f.fourier(x)
plot(xf$freq,xf$ampl^2,type="h",xlab="Frequenz", ylab="Quadr. Amplitude")


#9.ARCH
f.arch.sim <- function(n=500,init=200,a)
{
  ## Purpose: Simulating the Arch(1) process
  ## -------------------------------------------------------------------------
  ## Arguments: n=length of the series, init=length of the initial segment
  ## that is discarded, a= Parameter of the model
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  5 Feb 97, 16:49
  m <- n+init
  x <- rnorm(m)
  for (i in (2:m)) x[i] <- sqrt(1+a*x[i-1]^2)*x[i]
  ts(x[(init+1):m],start=1)
}
x <- f.arch.sim(n=256,a=0.8)

#10.Chaotisches System
f.logist <- function(n,start=runif(1))
{
  ## Purpose: Computing the iterated logistic map x(t+1)=4*x(t)*(1-x(t))
  ## -------------------------------------------------------------------------
  ## Arguments: start=starting value, length=number of iterates
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 11 Nov 96, 14:59
  x <- rep(0,n)
  x[1] <- start
  for (i in (2:n)) {x[i] <- 4*x[i-1]*(1-x[i-1])}
  x
}
x <- f.logist(256)
x <- x-mean(x)


#Randomisieren der Phase und Umkehrung: Veraendern der Phase hat keinen
#Einfluss auf das Spektrum. Man sollte daher so einen Prozess erhalten,
#der die gleichen Autokovarianzen hat, aber andere nichtlineare
#Zusammmenhaenge

f.randphase <- function(x)
{
  ## Purpose: Transforming a time series by randomizing the phase
  ## -------------------------------------------------------------------------
  ## Arguments: 
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 29 Nov 2000, 16:36
  x <- fft(x)
  n <- length(x)
  k <- (n-1)%/%2
  x[2:(k+1)] <- complex(mod=Mod(x[2:(k+1)]),arg=runif(k,-pi,pi))
  x[n:(n-k+1)] <- complex(mod=Mod(x[2:(k+1)]),arg=- Arg(x[2:(k+1)]))
  Re(fft(x,inverse=T))/n
}

x <- f.arch.sim(n=256,a=0.8)
y <- f.randphase(x)
plot(x)
plot(y)
acf(x)
acf(y)
acf(x^2)
acf(y^2)