#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