### periodogram of ARIMA (0,1,1) residuals
## first try
help(spectrum)
n <- length(res1)
y <- spec.pgram(res1, plot=F, spans = NULL)
y2 <- y$freq
y3 <- y$spec
y4 <- 10*log10(y3) ### to get decibels
y$bandwidth  ## .0013
# freq seems to be in cycles/datapoint.  
# multiply by # of datapoints per year (which is 12 here) to get cycles/year.
# here sigma^2 = var(res1) = .037
# To see which frequencies are "stat. sig.", use 95%-CIs for spectrum of 
# white noise.  For WN, f(w) = sigma^2/pi = .037/pi, with the frequencies
# defined in radians,  but here the freq is in cycles/datapoint so 
# f(w) = sigma^2. 
# the periodogram estimates (ordinarily times 2 pi/sigma^2 -- 
# but here just multiplied by 2/sigma^2 since the estimates are already 
# multiplied by pi) are chi-square (nu), see 
# Chatfield p. 113 or 121.  Here nu = 2.
# The critical levels of the chi-square are:
# P(chi-square_2 < .0506) = .025.  
# P(chi-square_2 < 7.378) = .975. 
# So put lower bound at .0506 * .037/2
# and upper at 7.378 * .037/2
# Also draw a line at sigma^2, denoting the spectrum for white noise.
# NOTE that this is different from ordinary 95%-CI for f(w)!
bounds <- c(.0506 * .037/2, 7.378 * .037/2)
y2 <- y$freq * 12
par(mfrow=c(2,1))
plot(range(y2),c(0,max(c(bounds,y3))),type="n",xlab="",ylab="")
lines(y2,y3)
mtext(s=1,l=2,"frequency (cycles/year)")
mtext(s=2,l=2,"amplitude")
mtext(s=3,l=2, "Fig 21:  Periodogram of ARIMA(0,1,1) resids")
abline(h=bounds,lty=3)  
abline(h=.037,lty=2)
spectrum(res1)




### smoothed spectrum
y <- spec.pgram(res1, plot=F,spans=9)
y2 <- y$freq
y3 <- y$spec
y$df	## 15.3
y$bandwidth  ## .011
y2 <- y$freq * 12
# use 95%-CIs for spectrum of white noise, 
# where f(w) = sigma^2 = .037, then use formula 
# Chatfield p. 121.  Here nu = 15.3.
# chi-square (15,.025) = 6.262.  chi-square(15, .975) = 27.49. 
# So put lower bound at 6.262 * .037/(15.3*2) 
# and upper at 27.49 * .037/(15.3*2)
# Also draw a line at sigma^2, denoting the spectrum for white noise.
# NOTE that this is different from ordinary 95%-CI for f(w)!
bounds <- c(6.262 * .037/15.3, 27.49 * .037/15.3)
par(mfrow=c(2,1))
plot(range(y2),c(0,max(c(bounds,y3))),type="n",xlab="",ylab="")
lines(y2,y3)
mtext(s=1,l=2,"frequency (cycles/year)")
mtext(s=2,l=2,"amplitude")
mtext(s=3,l=2, "Fig 22:  Smoothed p-gram of ARIMA(0,1,1) resids")
abline(h=bounds,lty=3)  
abline(h=.037,lty=2)
spec.pgram(res1,spans=9)




### AR spectrum
y <- spectrum(res1, method="ar", plot=F)
y2 <- y$freq
y3 <- y$spec
y2 <- y2 * 12
par(mfrow=c(2,1))
plot(y2,y3,type="n",xlab="frequency (cycles/year)",ylab="amplitude")
lines(y2,y3)
mtext(s=3,l=2, "Fig 23:  AR-spectrum estimate of ARIMA(0,1,1) resids")
spectrum(res1,method="ar")

### TRY THE ABOVE AGAIN, but with the 1st line:  
### y <- spectrum(x,method="ar",plot=F)
### AND THE LAST LINE:
### spectrum(x,method="ar")




### STL
y1 <- ts(res1,start=1,frequency=3)  
## since peak freq was about 4 cycles/yr = 1 cycle per 3 observations.
plot(stl(y1, "periodic"))

help(stl)

y1 <- ts(res1,start=1,frequency=12)
plot(stl(y1, "periodic"))