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

library(tseries)
source("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/rawspec.R")
source("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/smoothspec.R")

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

## Select the May monthly flow..
y=lf[,5]
deltaT = 1
t = 1:length(y)

#####  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 a best AR model to the data
    arbest = ar((y-mean(y)),order.max=25)
	p = order(arbest$aic)[1]

## AR-1 model
#armodel = ar((y-mean(y)),order.max=1)

armodel = arima((y-mean(y)), order = c(p,0,0), include.mean=TRUE, method="ML")

#2. Simulate from  the fitted model of the same length

nsim = 250

## matrix to hold the raw spectrum..
specsimrawar = matrix(0,ncol=nsim,nrow=nfreq)
specsimrawsur = matrix(0,ncol=nsim,nrow=nfreq)

### matrix to hold the smoothed spectrum
specsimrawars = matrix(0,ncol=nsim,nrow=nfreq)
specsimrawsurs = 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)

    ysim=arima.sim(n=N,list(ar=c(armodel$coef[1:p])), sd=sqrt(armodel$sigma2)) + mean(y)

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


#3. Compute the raw spectrum  of the simulations

specsimrawar[,isim]=rawspec(ysim,deltaT)
specsimrawsur[,isim]=rawspec(ysims,deltaT)

### Compute the smooth spectrum of the simulations..

specsimrawars[,isim]=smoothspec(ysim,deltaT)
specsimrawsurs[,isim]=smoothspec(ysims,deltaT)

#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
pgram=rawspec(y,deltaT)

#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

#### smooth spectrum of the data
pgrams = smoothspec(y,deltaT)

### 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)
############## Plot the smooth spectrum of the AR simulations
plot(plotfreqs,pgrams,xlab="Freq. cy/yr",ylab="Spectrum", ylim=range(specsimrawars,pgrams),col="red",lwd=3,
type="l")

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

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

############## Plot the smooth spectrum of the surrogate data
plot(plotfreqs,pgrams,xlab="Freq. cy/yr",ylab="Spectrum", ylim=range(specsimrawsurs,pgrams),col="red",lwd=3,
type="l")

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

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