#1. AR spectral density

load("/u/buhlmann/teaching/zeitreihen/demonstrations/ocwave.RData")
ocw <- ocwave
plot(ocw)

fit <- ar.yw(ocw)
fit

#AR-Spektrum for ocean wave data
sp <- spec.ar(fit)

peaks <- identify(sp$freq,sp$spec)
1/sp$freq[peaks]  
[1] 26.26316  #pseudo-period


data(sunspots)
plot(sunspots)

fit <- ar.yw(ts(sunspots))
fit

#AR-spectrum for sunspots data
sp <- spec.ar(fit)

peaks <- identify(sp$freq,sp$spec)
1/sp$freq[peaks]  
[1] 142.571429  12.632911   8.316667   5.117949 #in years



#2. Smoothing of the periodogram

par(mfrow=c(2,1))
spec.pgram(ocw)
spec.pgram(ocw,demean=F,detrend=F)

spec.pgram(ocw)
spec.pgram(ocw,demean=T,detrend=F)

#Tapering for bias reduction but leading to higher variance
spec.pgram(ocw,ylim=c(10^(-3),10^7))
spec.pgram(ocw,taper=0.4,ylim=c(10^(-3),10^7),plot=F)

par(mfrow=c(1,1))
sp  <- spec.pgram(ocw,spans=5)
#modify argument "spans" (even number) 
sp$kernel$coef

data("lynx")
spec.pgram(log(lynx))
sp <- spec.pgram(log(lynx),spans=3)
spar <- spec.ar(log(lynx))

#identification of periods which correspond to peaks in spectral density
peaks <- identify(sp$freq,sp$spec)
peaksar <- identify(spar$freq,spar$spec)

1/sp$freq[peaks]  #pseudo-periods in years
[1] 40.000000 10.000000  5.000000  3.333333  2.448980
1/spar$freq[peaksar]
[1] 27.722222  9.689320  4.965174  3.326667  2.526582


data("sunspots")
spec.pgram(sunspots)
sp <- spec.pgram(sunspots,spans=5)
#identification of periods which correspond to peaks in spectral density
peaks <- identify(sp$freq,sp$spec)

spar <- spec.ar(sunspots)
peaksar <- identify(spar$freq,spar$spec)

1/sp$freq[peaks] #pseudo-periods in years
[1] 10.9090909  1.1428571  0.7121662
1/spar$freq[peaksar] 
[1] 11.8809524  1.0527426  0.6930556  0.4264957

#3. Spektrum via overlapping segments

#10 segments for ocwave
sp <- spec.pgram(ocw[1:102],detrend=F,ylim=c(0.01,10^8))
averaged <- sp$spec
for(i in 2:10)
  {
        spi <- spec.pgram(ocw[((i-1)*102 + 1):(i*102)],detrend=F,plot=F)
            lines(sp$freq,spi$spec,col=i)
            averaged <- averaged + spi$spec
      }

averaged <- averaged/10

sp <- spec.pgram(ocw[1:102],detrend=F,ylim=c(0.01,10^8))
points(sp$freq,averaged,cex=1.5,type="b",col=2)


#with artificial trend
ocw1 <- 0.05*(1:1024) + 0.000005*(1:1024)^3 + ocw

sp <- spec.pgram(ocw1[1:102],detrend=F,ylim=c(0.01,10^8))
averaged <- sp$spec
for(i in 2:10)
  {
        spi <- spec.pgram(ocw1[((i-1)*102 + 1):(i*102)],detrend=F,plot=F)
            lines(sp$freq,spi$spec,col=i)
            averaged <- averaged + spi$spec
      }

averaged <- averaged/10

sp <- spec.pgram(ocw1[1:102],detrend=F,ylim=c(0.01,10^8))
points(sp$freq,averaged,cex=1.5,type="b",col=2)

#periodogram of trend
spec.pgram(0.05*(1:1024) + 0.000005*(1:1024)^3)
       #Trends "induce" spectral density with "pole" at zero