leesflow = read.table("LeesFerry-monflows-1906-2006.txt")
leesflow = leesflow[,-1]
leesflow = leesflow/10^6
nyrs = length(leesflow[,1])
X = leesflow[,5]		#May flow
DT=1
DJ=0.25
pad=1
zz=wavelet(X,DT,pad,DJ,6)
ns = length(zz$scale)
z2=matrix(0,nrow=ns,ncol=nyrs)
for(i in 1:ns){
z2[i,]=(sqrt(DT)*DJ)*Re(zz$wave[i,])/((sqrt(zz$scale[i]))*.776*(pi^(-1/4)))
}
Xrecon = apply(z2,2,sum) + mean(X)
## Xrecon and X should be similar..
## to check that the wavelet decomposition adds up to
## to the original data

#### Simulation...
## phase randomization...
### Easy and straightforwward
wavphas = zz$wave
for(i in 1:ns){
wavphas[i,]= Arg(zz$wave[i,])
}

sm.density(X,ylim=c(0,0.75))

index=1:nyrs

nsim=250

for(isim in 1:nsim){
index1=sample(index)

z2=matrix(0,nrow=ns,ncol=nyrs)
for(i in 1:ns){
z2[i,]=(sqrt(DT)*DJ)*Re(Mod(zz$wave[i,])*cos(wavphas[i,index1]))/((sqrt(zz$scale[i]))*.776*(pi^(-1/4)))
}
xmaysim = apply(z2,2,sum) + mean(X)

sm.density(xmaysim,add=T,col="grey")

}
sm.density(X,ylim=c(0,0.75),add=T,lwd=2,col="red")

### Fit AR to each component and simulate
sm.density(X,ylim=c(0,0.75))

nsim=250
for(isim in 1:nsim){

zs = c()

for(i in 1:ns){
Y = z2[i,] - mean(z2[i,])
  arbest = ar(Y,order.max=25)
	p = order(arbest$aic)[1]
	p=2
#zar = arima(Y, order = c(p,0,0), include.mean=TRUE, method="ML")
#zs1 = arima.sim(n=nyrs,list(ar=c(zar$coef[1:p])), sd=sqrt(zar$sigma2)) + mean(z2[i,])

zar = ar(Y,order=p)
zs1 = arima.sim(n=nyrs,list(ar=c(zar$ar)), sd=sqrt(zar$var.pred)) + mean(z2[i,])

zs=cbind(zs,zs1)
}

xmaysim = apply(zs,1,sum) + mean(X)
sm.density(xmaysim,add=T,col="grey")

}


### You can compute the band filtered series 
### Fit an AR or do phase randomization
### For the rest fit a Normal distribution