#Load the data and make the plot
data=scan("c:\\river-flows.txt")
plot(data,type="l")
#Compute seasonal components
data=ts(data,frequency=365)
pr1=cycle(data)
seasmean=tapply(data,pr1,mean)
seastd=tapply(data,pr1,var)
seastd=seastd^0.5
#Compute deseasonalised series
datad=(data-rep(seasmean,length=length(data)))/rep(seastd,length=length(data))
#Plot deseasonalised time series
plot(datad,type="l")
#Estimate the autoregressive model
datadd=c(0,datad[1:(length(data)-1)])
fi=0.1
res=function(fi)
{
sumsq=sum((datad-fi*datadd)^2)
return(sumsq)
}
pr2=optim(fi,res)$par
#Compute the residuals
epsilon=datad-pr2*datadd
#Plot autocorrelation function of epsilon
acf(epsilon)
#and check that the acf is lying within the band delimited by the 2 blue lines
#Apply the Kolmogorov-Smirnov test
#Compute cumulative frequency of not exceedance
cumfreq=ppoints(epsilon)
#Compute cumulative probability of the Gaussian distribution
cumprob=pnorm(sort(epsilon),mean(epsilon),sqrt(var(epsilon)))
statest=max(abs(cumfreq-cumprob))
#Check that states is lower or equal to maximum acceptable stats found on google
#Data generation
#Direct generation with R
sddata=arima.sim((365*50),model=list(ar=pr2),innov=rnorm((365*50),mean(epsilon),sqrt(var(epsilon))))
#Add seasonal components back
sdata=sddata*rep(seastd,50)+rep(seasmean,50)

#Computation of the annual flow duration curves

fdc=array(0,dim=c(50,365))
for(i in 1:50)
{
pr3=rev(sort(sdata[((i-1)*365+1):(i*365)]))
fdc[i,]=pr3
}
#Plotting the group of the 50 fdc
plot(seq(1,365,1),fdc[1,],type="l",ylim=c(0,max(fdc)))
for (i in 1:50)
lines(seq(1,365,1),fdc[i,])

#You may check if negative values are there.
sdata[sdata<0]=0