# R function for computing the confidence limits of river flows simulated
# by rainfall-runoff models. The routine applies the meta-Gaussian approach
# by Montanari and Brath, "A stochastic approach for assessing the uncertainty
# or rainfall-runoff simulations", Water Resources Research, 40, 
# doi:10.1029/2003WR002540, 2004
# This program runs under R (see http:\\www.r-project.org). R is public domain
# and is very easy to run this program in it. If you do not have any knowledge
# of R, please contact me (alberto.montanari@unibo.it). I can profide you
# with a quick list of instructions that will allow you to get the confidence
# bands you are looking for.

# Please ignore the warning: Collapsing to unique x values in: approx(....)
# It odes not have any effect

# Please keep in mind that when performing the numerical NQT linear extrapolation
# is NOT performed when xout is external to the range of x. Therefore it is strongly
# recommended that xout is never outside the x range. If you need futher details
# please contact alberto.montanari@unibo.it

# The users needs to put in the command line 2 records of observed (dataobs) and
# simulated (modsim) river flows to calibrate the meta-Gaussian approach. The data
# to compute the confidence bands of are placed in the vector "data". Select
# transf=T (the default) if you want to perform preliminary transformation of the
# data as described in Montanari and Brath (2004) cited above. Select plot=T 
# (the default) if you want a plot of the residuals as the ones reported in the
# same paper.
#PLEASE NOTE: the sotware allows the user to specify a lower threshold value for the
#simulated river flow below which the model error is not analyzed. This option
#allows the user to limit the effect of non-stationarity (see Montanari & Brath, 2004).
#It is highly recommended to specify such a threshold, by also looking at the 
#goodness of fit provided by the meta-Gaussian model

# See the note here below if you wish to perform the probability plot correlation
# coefficient test

# PLEASE NOTE that this code is a beta. I tested it and it seems to work well,
# but of course I cannot be sure that it is not affected by bugs. Please report
# to me any problem you experience with it. Thanks.

uncert=function(dataobs,modsim,data,thr=0,trasf=T,plot=T)
{
dataobs=dataobs[modsim>thr]
modsim=modsim[modsim>thr]
err=dataobs-modsim
kasp=mean(err)
err1=abs(err-kasp)

pr1=ppoints(modsim)
pr2=qnorm(pr1)
pr3=sort(modsim)
pr8=approx(pr3,pr2,xout=modsim,rule=2)$y
Ndata=approx(pr3,pr2,xout=data,rule=2)$y
Nmodsim=pr8

if(trasf==T)
{
pr1=ppoints(err1)
pr5=qnorm(pr1)
pr6=sort(err1)
pr7=approx(pr6,pr5,xout=err1)
Nerr=pr7$y

pre1=lsfit(Nmodsim,Nerr)
residui=pre1$residuals
mcond=pre1$coef[[1]]+pre1$coef[[2]]*Nmodsim

corcof=acf(cbind(Nmodsim,Nerr),lag.max=1,plot=F)$acf[1,2,1]

pr10=approx(pr5,pr6,xout=mcond,rule=2)$y

err1=(err-kasp)/pr10+kasp
}
else
{
err1=err
}

pr101=ppoints(err1)
pr105=qnorm(pr101)
pr106=sort(err1)
pr107=approx(pr106,pr105,xout=err1)
Nerr=pr107$y

mNerr=mean(Nerr)
sNerr=(var(Nerr))^0.5
mNmodsim=mean(Nmodsim)
sNmodsim=(var(Nmodsim))^0.5

corcof=acf(cbind(Nmodsim,Nerr),lag.max=1,plot=F)$acf[1,2,1]

mcond=mNerr+corcof*sNerr/sNmodsim*(Ndata-mNmodsim)
scond=sqrt(sNerr^2*(1-corcof^2))

conflim=0.95

pr9=mcond+qnorm(-(1-conflim)/2+1)*scond
pr9a=mcond-qnorm(-(1-conflim)/2+1)*scond

pr700=approx(pr105,pr106,xout=pr9,rule=2)$y
pr701=approx(pr105,pr106,xout=pr9a,rule=2)$y


if(trasf==T)
{
pr700=(pr700-kasp)*approx(modsim,pr10,xout=data,rule=2)$y+kasp+data
pr701=(pr701-kasp)*approx(modsim,pr10,xout=data,rule=2)$y+kasp+data
}

resw=Nerr-corcof*Nmodsim

#If you want to perform propability plot correlation coefficient test
#please remove the # symbols in the 5 lines here below

#pr50=pnorm(sort(resw),mean(resw),sqrt(var(resw)))
#pr51=ppoints(resw)
#cat('Probability plot correlation coefficient test:',fill=T)
#cat('Test statistic:',max(abs(pr50)-abs(pr51)),fill=T)
#cat('95% upper limit of the test statistic:',+1.36/sqrt(length(resw)),fill=T)
if(plot==T)
{
plot(corcof*Nmodsim,resw)
}

result=list()
result$data=data
result$csup=pr700
result$cinf=pr701
return(result)
}