# March 20 2001
# Script for periodogram examples for spectral density course

# library(ts)  # need ts

# make a plot of the raw periodograms for the lynx monthly data
lynx.plt1()  

# euler_E( log(X) ) where X ~ exp,1

euler_-0.577216
var.log.exp_pi^2/6
var.exp_1

lynx.pgram_pgram(lynx.dat$dat)

#*****************************
#  Daniell smooth periodograms

f.hat.plt.eg1_function(m=3){
txt_paste( c("Daniell Smooth, m = ",m),collapse="")
k.daniell_kernel("daniell",m=m)
f.hat_kernapply(lynx.pgram$pgram/(2*pi),k.daniell)
freq_lynx.pgram$freq[m:(56-m)]
f.ci.up_f.hat*(1+2*sqrt(var.exp/(2*m+1)))
f.ci.low_f.hat*(1-2*sqrt(var.exp/(2*m+1)))

ylim_c(0, max(f.hat,f.ci.up))
plot(freq,f.hat,type="l",main=txt,ylim=ylim)
lines(freq,f.ci.up,lty=2)
lines(freq,f.ci.low,lty=2)
}

f.hat.plt.eg1(m=3)
f.hat.plt.eg1(m=5)


#*************
#  Daniell smooth log periodograms


log.f.hat.plt1_function(m=3){
 txt_paste( c("Log Periodogram, Daniell Smooth, m = ",m),collapse="")
 k.daniell_kernel("daniell",m=m)
 freq_lynx.pgram$freq[m:(56-m)]

 log.f.hat_kernapply(log(lynx.pgram$pgram/(2*pi)),k.daniell)
 log.f.ci.up_log.f.hat+2*sqrt(var.log.exp/(2*m+1))
 log.f.ci.low_log.f.hat-2*sqrt(var.log.exp/(2*m+1))

 ylim_c(min(log.f.hat,log.f.ci.low), max(log.f.hat,log.f.ci.up))
 plot(freq,log.f.hat,type="l",main=txt,ylim=ylim)
 lines(freq,log.f.ci.up,lty=2)
 lines(freq,log.f.ci.low,lty=2)
}

log.f.hat.plt1(m=3)
log.f.hat.plt1(m=7)

#**************
# simulate an AR series and make some smooth periodogram plots

n_1000
beta0_.6
x_arma.sim(n=n,ar=c(beta0))
x.pgram_pgram(x)

# true spectral density
spec.dens_function(lam,beta){
 1/((2*pi)*abs(1-beta*exp(-lam*1i))^2)
}

f.hat.plt.eg2_function(m=3,x.pgram){
 par(oma=c(0,0,4,0))
 txt_paste( c("Daniell Smooth, m = ",m),collapse="")
 k.daniell_kernel("daniell",m=m)
 f.hat_kernapply(x.pgram$pgram/(2*pi),k.daniell)
 freq_x.pgram$freq[m:(length(x.pgram$pgram)-(m+1))]
 f.ci.up_f.hat*(1+2*sqrt(var.exp/(2*m+1)))
 f.ci.low_f.hat*(1-2*sqrt(var.exp/(2*m+1)))

 ylim_c(0, max(f.hat,f.ci.up))
 cat(length(freq),", ",length(f.hat),"\n")
 plot(freq,f.hat,type="l",main=txt,ylim=ylim)
 lines(freq,f.ci.up,lty=2)
 lines(freq,f.ci.low,lty=2)
 mtext("AR Simulation Example Smoothed Periodogram", outer=T,side=3,line=3,cex=1.5)
}



f.hat.plt.eg2(3,x.pgram)
lam_seq(0,.5,.01)
lines(lam,spec.dens(2*pi*lam,beta0),lty=1)