# transfer function of (seasonal) differences
f.transfer <- function(coef,freq)
{
  ## Purpose: computing \sum coeff(u)*exp(-i*freq*u) 
  ##          (so-called transfer function)
  ## -------------------------------------------------------------------------
  ## Arguments: coef=impulse resoponse coefficients, freq=frequency
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  2 Dec 96, 11:11
  p <- length(coef)-1
  z <- (cos(freq)+sin(freq)*1i)
  a  <- outer(z,0:p,FUN="^")
  a%*%coef}
nu <- seq(-0.5,0.5,length=1001)
coef <- c(1,-1) #first difference
coef <- c(1,-2,1) #second difference
coef <- c(1,rep(0,11),-1) # seasonal difference
A <- f.transfer(coef,nu*2*pi)
par(mfrow=c(1,2))
plot(nu,Mod(A),type="l",xlab="frequency",ylab="amplitude modulation")
plot(nu,Arg(A),type="l",xlab="frequency",ylab="phase shift")
# phase shift is an odd function, amplitude modulation an even function

# Properties of the periodogram
# Functions: fft for Fast Fourier transform spec.pgram for periodogram

help(fft)
help(spec.pgrm)

f.fourier <- function(x)
{
  ## Purpose: Discrete Fourier transform given as amplitude and phase
  ## -------------------------------------------------------------------------
  ## 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))}

# White noise: Checking of exponential distribution and independence
x <- rnorm(512)
pgram <- (f.fourier(x)$ampl)^2/length(x)
par(mfrow=c(2,2))
k <- length(pgram)-1
hist(pgram[2:k],nclass=20,xlab="Periodogram values",main="")
plot(-log(1-(1:(k-1))/k),sort(pgram[2:k]),xlab="Exponential quantiles",
     ylab="Ordered periodogram values")
abline(a=0,b=1)
plot(pgram[2:(k-1)],pgram[3:k],xlab="Periodogram",
     ylab="Periodogram shifted by 1")
plot(pgram[2:(k-2)],pgram[4:k],xlab="Periodogram",
     ylab="Periodogram shifted by 2")
# In logarithmic scale
lpgram <- log(pgram)
par(mfrow=c(3,1))
freq <- 0:(length(x)%/%2)/length(x)
plot(freq,pgram,type="l",xlab="frequency", ylab="squared amplitude")
plot(freq,lpgram,type="l",xlab="frequency", ylab="log(Amplitude^2)")
spec.pgram(x,taper=0,demean=TRUE,detrend=FALSE)
    # same thing except scale on y axis and value for frequency zero omitted
par(mfrow=c(2,2))
hist(lpgram[2:k],nclass=20,xlab="log Periodogram values",main="")
plot(log(-log(1-(1:(k-1))/k)),sort(lpgram[2:k]),xlab="Gumbel quantiles",
     ylab="Ordered log periodogram values")
abline(a=0,b=1)
plot(lpgram[2:(k-1)],lpgram[3:k],xlab="log Periodogram",
     ylab="log eriodogram shifted by 1")
plot(lpgram[2:(k-2)],lpgram[4:k],xlab="log Periodogram",
     ylab="log periodogram shifted by 2")

# periodogram of log(lynx)
spec.pgram(log(lynx),taper=0,demean=TRUE,detrend=FALSE)
# Can you guess what the true spectral might be ?
spec.ar(ar.mle(log(lynx)),add=TRUE,col=2)
# spectral density from a fitted AR model added. Agreement
# with periodogram seems reasonable (notice the confidence band
# given by the blue segments at the right boundary).

# Moving average: Consider the ratio periodogram:spectrum which should
#                 still be i.i.d. exponential

coef <- cos(0.3*pi*(0:7)) # choose some coefficients 
coef <- coef/sqrt(sum(coef^2)) # normalisation
x <- filter(rnorm(520),filter=coef,sides=1)[9:520] #generate MA
pgram <- (f.fourier(x)$ampl)^2/length(x)
freq <- 2*pi*(0:(length(x)%/%2))/length(x)
pgram <- pgram/Mod(f.transfer(coef,freq))^2
# Now repeat the plots that were done with white noise.

#AR(4): Example which shows the benefit of tapering:
ar.coef <- c(2.7607,-3.8106,2.6535,-0.9238)
Mod(polyroot(c(1,-ar.coef))) #causality condition is satisfied
2*pi/Arg(polyroot(c(1,-ar.coef))) # periods of associated damped harmonics
x <- arima.sim(n=1024,model=list(ar=ar.coef))
par(mfrow=c(2,1))
plot(x)
acf(x)
par(mfrow=c(1,1)) #various amounts of tapering
spec.pgram(x,taper=0,demean=TRUE,detrend=FALSE)
spec.pgram(x,taper=0.1,demean=TRUE,detrend=FALSE,add=TRUE,col=2)
spec.pgram(x,taper=0.3,demean=TRUE,detrend=FALSE,add=TRUE,col=3)
# How do tapered and non-tapered estimates differ ?
# Can you see the two peaks of the spectrum ?
abline(v=c(1/7.14,1/9.09))
# computing and plotting the true spectrum
nu <- seq(0,0.5,length=501)
A <- f.transfer(c(1,-ar.coef),nu*2*pi)
lines(nu,1/Mod(A)^2,col=4)
# fitting an AR(4) and plotting its spectrum
spec.ar(ar.mle(x,order.max=4),add=TRUE,col=5)
# almost the same thing as the true spectrum


# Fejer and other kernels that determine the expected value of
# the (tapered) periodogram
n <- 10 #series length
p <- 0.4 #percentage of tapering on both sides (split cosine taper)
h <- rep(1,n)
r <- floor(0.5+n*p)
h[1:r] <- 0.5*(1-cos(pi*(2*(1:r)-1)/(2*n*p)))
h[n+1-(1:r)] <- h[1:r]
plot(h) #shows the taper
# Left: Fejer kernel (no tapering),
# Right: Kernel corresponding to 40% tapering on both sides
# Sample sizes 10 (upper) and 100 (lower)
# Tapering reduces the side peaks, but increases the width of the main peak
par(mfrow=c(2,2))
for (n in c(10,100)) {
  for (p in c(0,0.4)) {
    h <- rep(1,n)
    r <- floor(0.5+n*p)
    if (r > 0) {
      h[1:r] <- 0.5*(1-cos(pi*(2*(1:r)-1)/(2*n*p)))
      h[n+1-(1:r)] <- h[1:r]
    }
    N <- 4096 # determines the resolution for plotting 
    h <- c(h,rep(0,N-n))
    H <- Mod(fft(h))^2
    H <- c(H[(N/2+1):2],H[1:(N/2+1)])/sum(h^2)
    plot(seq(-0.5,0.5,length=N+1,),H,type="l",
         ylim=c(-0.5,n),xlab="frequency",ylab="kernel")
  }}