{
#simple periodogram analysis..
#y is the vector with the data.

data=matrix(scan("hwk2-1a.txt"),ncol=13,byrow=T)
data=data[,2:13]/10000                   #remove year and divide
y=data
y=array(t(y))
N=length(y)

c=acf(y,type="covariance",lag.max=(N-1))$acf

freqs=2*pi*(1:round((N/2)))/N

plotfreqs=freqs/(2*pi*(1/12))

tordf=1:round((N/2))

tord=1:N

nfreq=length(freqs)

ap=1:nfreq
bp=1:nfreq
for(i in 1:(nfreq-1)){

ap[i]=2*sum(y*cos(freqs[i]*tord))/N
bp[i]=2*sum(y*sin(freqs[i]*tord))/N

}

ap[nfreq]=sum(y*((-1)^(tord)))/N

bp[nfreq]=0

pgrams=N*((ap*ap) + (bp*bp))/(4*pi)
pgrams[nfreq]=N*ap[nfreq]*ap[nfreq]/pi

par(mfrow=c(3,1))

plot(plotfreqs,pgrams,xlab="freq. cy/yr",ylab="Pgram", type="h")

#ACF based..

N1=length(c)

pgramsacf=1:nfreq
for(i in 1:nfreq){
pgramsacf[i]=(c[1]+(2*sum(c[2:N1]*cos(freqs[i]*(1:(N1-1))))))/pi
}
plot(plotfreqs,pgramsacf,xlab="freq. cy/yr",ylab="Pgram", type="h")


#FFT of the data..
pgrams1=Mod(fft((y/N)))^2
pgrams1=N*pgrams1/pi
plot(plotfreqs,pgrams1[2:(nfreq+1)],xlab="freq. cy/yr",ylab="Pgram", type="h")


#Using equation 7.16 of Chatfield p. 128
pgramseq=1:nfreq

for(i in 1:nfreq){
ap[i]=(sum(y*cos(freqs[i]*tord)))^2
bp[i]=(sum(y*sin(freqs[i]*tord)))^2
pgramseq[i]=(ap[i]+bp[i])/(N*pi)
}

#plot(plotfreqs,pgramseq,xlab="freq. cy/yr",ylab="Pgram", type="h")

}