source("conflines.txt")
lf = matrix(scan("LeesFerry-monflows-1906-2006.txt"), ncol=13, byrow=T)
lf = lf[,2:13]
lf = lf*0.0004690502	# convert to cms
y = lf
y=array(t(y))
deltaT = 1/12	# for monthly data..

## if you want to perform spectrum on monthly data comment the below three
## lines

lfann = apply(lf,1,sum)
y=lfann
deltaT = 1


#####  Get the ak and bk and it is equivalent to
#####  to the FFT  - see below

# number of frequencies possible..
N = length(y)
freqs=2*pi*(1:round((N/2)))/N
plotfreqs=freqs/(2*pi*(deltaT))
nfreq=length(freqs)
a = 1:nfreq
b = 1:nfreq

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

t=1:length(y)

for(k in 1:nfreq){
a[k]=2*sum(y*cos(2*pi*k*t/N))/N
b[k]=2*sum(y*sin(2*pi*k*t/N))/N
}

fspecraw = sqrt(a*a + b*b) / 2

Ck = fft(y)/N		#fourier coefficients r - does fft without dividing by N
fspecrawfft = Mod(Ck)
fspecrawfft = fspecrawfft[2:(nfreq+1)]   #R does summation 0 to N-1 when you need 1 to N

## fspecrawfft and fspecraw will be the same

######################################################

###### Smooth the raw spectrum using a window function
########

N=length(y)

# compute the auto covariance..
nlag = N-1
c=acf(y,type="covariance",lag.max=nlag)$acf

#lag window width
M = round(N/20)		# M should be even
if(M %% 2 > 0)M=M+1
# try different Window widths and different window types....

#Parzen Window..
lambdas=1:M
lambdas[1:(M/2)]=1-(6*((lambdas[1:(M/2)]/M)^2))+ (6*((lambdas[1:(M/2)]/M)^3))
lambdas[((M/2)+1):M]=2*((1-(lambdas[((M/2)+1):M]/M))^3)
lambda0=1
degf=3.709*N/M

#Tukey window..
lambdas=1:M
lambdas[1:M]=0.5*(1+cos(pi*(1:M)/M))
lambda0=1
degf=8*N/(3*M)

# number of frequencies possible..
freqs=2*pi*(1:round((N/2)))/N
plotfreqs=freqs/(2*pi*(deltaT))
nfreq=length(freqs)

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

# compute the smoothed periodgram/spectrum.. - Equation 12.3.15 Wei

pgrams=1:nfreq

for(i in 1:nfreq){
pgrams[i]=((lambda0*c[1]) + (2*sum(lambdas*c[2:(M+1)]*cos(freqs[i]*(1:M)))))/pi
}



cc=acf(y,type="covariance",lag.max=N-1,plot=FALSE)$acf
NN=length(cc)
NN1=NN-1
crev=rev(cc)
cfull=c(crev[2:NN],cc)
Nfull=length(cfull)
### Spectrum from sample ACF - 12.2.1b Wei
pgramacf = 1:nfreq
NN=length(cc)
for(i in 1:nfreq){
pgramacf[i]=(cc[1]+(2*sum(cc[2:NN]*cos(freqs[i]*(1:(NN-1))))))/(2*pi)
}




#confidence intervals.. based on Chisquare..

lowfac=degf/qchisq(0.975,degf)
upfac=degf/qchisq(0.025,degf)
spfft.u=pgrams*upfac
spfft.l=pgrams*lowfac



plot(plotfreqs,pgrams,xlab="freq. cy/yr",ylab="Spectrum", type="l", ylim=range(c(pgrams,pgramacf,spfft.l,spfft.u)))
lines(plotfreqs,pgramacf,type="h")
conflines(plotfreqs,spfft.u,spfft.l)



### To get the significance of a peak - you can do it couple of ways
## as mentioned in the class

#1. Fit an AR-1 model to the data
armodel = ar((y-mean(y)),order.max=1)
#2. Simulate from  the fitted model of the same length

nsim = 250
specsim = matrix(0,ncol=nsim,nrow=nfreq)

for(isim in 1:nsim){

ysim = arima.sim(n=N, list(ar=c(armodel$ar)), sd=sqrt(armodel$var.pred)) + mean(y)

#3. Compute the spectrum  in the same manner as you did for the data
# above..

c=acf(ysim,type="covariance",lag.max=nlag,plot=FALSE)$acf
pgramsim=1:nfreq
for(i in 1:nfreq){
pgramsim[i]=((lambda0*c[1]) + (2*sum(lambdas*c[2:(M+1)]*cos(freqs[i]*(1:M)))))/pi
}
specsim[,isim]=pgramsim

}




#4. Repeat - obtain the 95th percentile of the spectrum at each
#frequency. plot it.
yconflev = 1:nfreq
for(i in 1:nfreq){yconflev[i] = quantile(specsim[i,],c(0.95))}




plot(plotfreqs,pgrams,xlab="freq. cy/yr",ylab="Spectrum", type="l", ylim=range(c(pgrams,pgramacf,spfft.l,spfft.u,
yconflev)),col="blue",lwd=2)
lines(plotfreqs,pgrams,type="h")
conflines(plotfreqs,spfft.u,spfft.l)
lines(plotfreqs,yconflev,col="red")