#install.packages("astsa")
require(astsa)

spaic = spec.ar(soi, log="no")
cbind(spaic$freq, spaic$spec)[1:40,] # inspect, to find the El Nino frequency
en =  spaic$freq[19]
abline(v=en, lty="dotted")


soi.ar = ar(soi, order.max=30)
win.graph()
plot(1:30, soi.ar$aic[-1], type="o", ylab = "AIC")

soi.ar = ar(soi, order.max=60)
win.graph()
plot(1:60, soi.ar$aic[-1], type="o", ylab = "AIC")

# Compare with the periodogram
k = kernel("daniell", m=2)
win.graph()
soi.ave = spec.pgram(soi, k, log = "no")

#plot both
win.graph()
par(mfrow=c(1,2))
spaic = spec.ar(soi, log="no")
soi.ave = spec.pgram(soi, k, log = "no")