# These are the spectral analysis programs I demonstrated in class
# on 10/07/08.

# As with other test programs, you may have to adjust the "scan"
# statements to include the full path to the file

# Example 1: Lake Huron data

y1<-scan("lake.tsm")
plot(y1,type='l',xlab='Year',ylab='Level of Lake')

# Fit and display linear trend
x1<-1:length(y1)
reg1<-lm(y1~x1)
lines(x1,reg1$fitted)
summary(reg1)

# calculate residuals, display acf and pacf
y2<-reg1$resid
acf(y2)
pacf(y2)

# calculate spectrum by AR method (the program automatically
# chooses AR(2) as best

spectrum(y2,method='ar')

# periodogram method - without smoothing

spectrum(y2,method='pgram')

# or (this gives the same result)

spec.pgram(y2)

# Periodogram with various smoothing

spec.pgram(y2,spans=c(3,3))
spec.pgram(y2,spans=c(7,7))
spec.pgram(y2,spans=c(11,11))


# Example 2: Mauna Loa CO2 data

y3<-scan("maunaloa.txt")
plot(y3,type='l')

spec.pgram(y3)
spec.pgram(y3,spans=c(5,5))
spec.pgram(y3,spans=c(13,13))
spec.pgram(y3,spans=c(21,21))

# remove both linear and quadratic trends, plot residuals

x3<-1:length(y3)
x32<-x3^2
y4<-lm(y3~x3)$resid
plot(y4,type='l')
y4<-lm(y3~x3+x32)$resid
plot(y4,type='l')

# spectral densities for residuals from quadratic trend

spec.pgram(y4)
spec.pgram(y4,spans=c(5,5))
spec.pgram(y4,spans=c(13,13))
spec.pgram(y4,spans=c(21,21))

# repeat with tapering

par(mfrow=c(2,1))
spec.pgram(y4,spans=c(5,5),taper=0)
spec.pgram(y4,spans=c(5,5),taper=0.5)

par(mfrow=c(1,1))

# Global climate data

y5<-scan("uea1.txt")
plot(y5,type='l')
y6<-y5[609:1632]
plot(y6,type='l')

spec.pgram(y6)
spec.pgram(y6,spans=c(5,5))
spec.pgram(y6,spans=c(13,13))
spec.pgram(y6,spans=c(21,21))