#  Spectral estimation

# SS code

source("http://www.stat.pitt.edu/stoffer/tsa2/Rcode/itall.R")


# simulated AR(4) with two peaks

n <- 1024
phi <- c(2.7607, -3.8106, 2.6535, -0.9238)
true.f <- spec.arma(phi)

ar4.data <- arima.sim(model=list(ar=phi), n=1024)

# raw periodogram shows leakage effects
pgram <- spec.pgram(ar4.data, taper=0.0, demean=FALSE, detrend=FALSE, plot=FALSE) 
lines(pgram$freq, pgram$spec)

# slight tapering improves, but still variable
pgram.tapered <- spec.pgram(ar4.data, taper=0.1, demean=FALSE, detrend=FALSE, plot=FALSE)
lines(pgram.tapered$freq, pgram.tapered$spec, col="red")


################################################################
# read data (variable star)

xt <- scan("http://www-stat.wharton.upenn.edu/~stine/stat910/data/varstar.dat")
plot(xt, type="o")

# spectral estimate

p2 <- spec.pgram(xt, demean=TRUE, taper=0.1, plot=FALSE); 
lines(p2$freq,p2$spec, col="lightblue")

plot(p2$freq,p2$spec,log="y", type="l")

p3 <- spec.pgram(xt, demean=TRUE, taper=0.1, plot=FALSE, spans=c(5,3)); 
lines(p3$freq,p3$spec, col="red")


################################################################
#
# multitaper methods  (sapa package)
#
# library(sapa)
#
#

n    <- 1024
time <- 0:(n-1)

# plot taper functions  (discrete prolate spheroid)

h.k <- taper(type="dpss", n.sample=n, n.taper=8, bandwidth=4)   # bw = number multiples of FF
dim(h.k)                                                        # each row is taper
round(h.k %*% t(h.k), digits=4)                                 # orthogonal

stackPlot(time, t(h.k[1:4,]))            # stackPlot loaded with sapa
stackPlot(time, t(h.k[5:8,]))                       


# FT of data tapers 

ft.k <- Mod(t(mvfft(t(h.k))))

freq <- 1:32
plot (time[freq],ft.k[1, freq],log="y",type="l")
lines(time[freq],ft.k[2, freq],      col="red")
lines(time[freq],ft.k[3, freq],      col="green")
lines(time[freq],ft.k[4, freq],      col="blue")


# FT of tapered data; note eventual presence of leakage

ar4.data <- arima.sim(model=list(ar=phi), n=1024)

true.f <- spec.arma(phi)
s1 <- spec.pgram(h.k[1,] * ar4.data, taper=0.0, plot=FALSE); lines(s1$freq, n * s1$spec)            
s2 <- spec.pgram(h.k[2,] * ar4.data, taper=0.0, plot=FALSE); lines(s2$freq, n * s2$spec, col="red")
s3 <- spec.pgram(h.k[3,] * ar4.data, taper=0.0, plot=FALSE); lines(s3$freq, n * s3$spec, col="green")
s4 <- spec.pgram(h.k[4,] * ar4.data, taper=0.0, plot=FALSE); lines(s4$freq, n * s4$spec, col="blue")

true.f <- spec.arma(phi)
s5 <- spec.pgram(h.k[5,] * ar4.data, taper=0.0, plot=FALSE); lines(s5$freq, n * s5$spec)            
s6 <- spec.pgram(h.k[6,] * ar4.data, taper=0.0, plot=FALSE); lines(s6$freq, n * s6$spec, col="red")
s7 <- spec.pgram(h.k[7,] * ar4.data, taper=0.0, plot=FALSE); lines(s7$freq, n * s7$spec, col="green")            
s8 <- spec.pgram(h.k[8,] * ar4.data, taper=0.0, plot=FALSE); lines(s8$freq, n * s8$spec, col="blue")


# Estimator obtained by averaging

true.f <- spec.arma(phi)
s.bar = (s1$spec+s2$spec+s3$spec+s4$spec+s5$spec+s6$spec)/6
lines(s1$freq, n*s.bar, col="seagreen")
lines(pgram.tapered$freq, pgram.tapered$spec, col="red")


# built-in routines  (see docs with sapa package)
# You can see the increasing variance due to leakage at higher frequency in the plot

ff <- SDF(as.vector(ar4.data), method="multitaper", taper.=h.k);  plot(ff)