#install.packages("astsa") require(astsa) #help(package = astsa) ma = filter(soi, sides=2, c(.5, rep(1,11), .5)/12) # moving average ts.plot(soi, ma, lty=2:1, col=1:2, lwd=1:2, main="Southern Oscillation Index with 12 month moving average") #p = 15 #arfit = sarima(soi, p,0,0, diagnostics = F, details = F) #z = c(1, -arfit$fit$coef[1:p]) #sort(round(abs(polyroot(z)),2)) # Ordered, rounded moduli of zeros of the char. eq'n. ############################### win.graph() par(mfrow=c(2,1)) soi1.per = spec.pgram(soi, log="no", main = "SOI spectrum", xlab = "cycles/unit time") mat = cbind(soi1.per$freq/12, soi1.per$spec) # cycles/month colnames(mat) = c("freq", "power") mat abline(v= soi1.per$freq[40], lty=2) abline(v= soi1.per$freq[8], lty=2) soi2.per = spec.pgram(soi, log="yes", main = "SOI log-spectrum", xlab = "cycles/year") abline(v= soi2.per$freq[40], lty=2) abline(v= soi2.per$freq[8], lty=2) ############################### win.graph() par(mfrow=c(2,1)) k = kernel("daniell",1) # L = 3, m = 1 soi.ave = spec.pgram(soi, k, main = "smoothed log-spectrum; L=3") soi.ave$df k = kernel("daniell",5) # L = 11, m = 5 soi.ave = spec.pgram(soi, k, main = "smoothed log-spectrum; L=11") soi.ave$df ############################### # Add 95% confidence limits to the periodogram with L=9: win.graph() k = kernel("daniell",4) # L = 9, m = 4 soi.ave = spec.pgram(soi, k, plot = F) # Don't plot it yet Upper = qchisq(.975, soi.ave$df) # Upper chisquare quantile Lower = qchisq(.025, soi.ave$df) # Lower chisquare quantile up = log(soi.ave$df/Lower) # Add this to log(power) to get the upper conf. limit down = log(Upper/soi.ave$df) # Subtract this from log(power) to get the lower conf. limit upper.limit = ts(log(soi.ave$spec) + up) lower.limit = ts(log(soi.ave$spec) - down) plot(soi.ave$freq,log(soi.ave$spec), type = "l", lty=1, col = 2, ylim = c(min(lower.limit), max(upper.limit)), xlab = "cycles/year", ylab = "log(power)", main = "log-spectrum for SOI, with 95% CIs") lines(soi.ave$freq,upper.limit, lty=1,col=1, lwd=2) lines(soi.ave$freq,lower.limit, lty=1, col=1, lwd=2) ############################### win.graph() x = ts(cbind(soi,rec), frequency = 12) L = 19 s = spec.pgram(x, kernel("daniell",(L-1)/2), plot=F) plot(s, plot.type = "coh", ci.lty = 2, xlab = "cycles/year", main = "Coherence: SOI & Recruits") # These last two lines could be combined: s = spec.pgram(x, kernel("daniell",(L-1)/2), plot.type = "coh", ci.lty = 2, xlab = "cycles/year", main = "Coherence: SOI & Recruits") # Add a horizontal line corresponding to a p-value of .001: s$df # df = 32.12546 f = qf(.999, 2, s$df-2) # f = 8.764348 c = f/(.5*s$df-1+f) # c = 0.3678314 abline(h = c) ###############################