### Simulate from the raw spectrum ############

### check on the raw periodogram..

library(tseries)

lf = matrix(scan("LeesFerry-monflows-1906-2006.txt"), ncol=13, byrow=T)
lf = lf[,2:13]
lf = lf*0.0004690502	# convert to cms

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
#####

######################3

#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
specsimrawar = matrix(0,ncol=nsim,nrow=nfreq)
specsimrawsur = matrix(0,ncol=nsim,nrow=nfreq)
for(isim in 1:nsim){

## simulate from AR1 model

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

## simulate from spectrum
ysims = surrogate(y,fft=TRUE,amplitude=TRUE)


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


pgramsim=1:nfreq
pgramsims=1:nfreq

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

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


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

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



specsimrawar[,isim]=fspecraw
specsimrawsur[,isim]=fspecrawsur


}

## spectrum of the original data


#### raw spectrum of original data

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
}

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


### Plot the raw spectrum from AR simulations

plot(plotfreqs,pgram,xlab="Freq. cy/yr",ylab="Spectrum", ylim=range(specsimrawar,pgram),col="red",lwd=2)

for(i in 1:nsim){
lines(plotfreqs, specsimrawar[,i],col="grey")
}

lines(plotfreqs,pgram,col="red",lwd=3)

## plot the raw spectrum from surrogate data - i..e, simulations from spectrum

plot(plotfreqs,pgram,xlab="Freq. cy/yr",ylab="Spectrum", ylim=range(specsimrawsur,pgram),col="red",lwd=3,
type="l")

for(i in 1:nsim){
lines(plotfreqs, specsimrawsur[,i],col="grey")
}

lines(plotfreqs,pgram,col="red",lwd=3)