tvt<-function(y,span,index,prob){

# Compute time-varying threshold for a given monthly time series
#
#  Description:
#
#       Returns smooth long term mean, time-varying threshold and fitted values
#       given by long term mean plus mean seasonal cycle
#
# Usage: svt(y,span,index,prob)
#
# Input:
#
#       y: vector containing a monthly time series of data
#
#    span: fraction of the total number of points of 'y' to be used to compute
#          the long term mean
# 
#   index: index vector composed of logicals (TRUE and FALSE) of the same 
#          length as 'y' defining the months to be considered when computing the
#          time-varying threshold
#
#    prob: a value between 0 and 1 that informs the percentage of point
#          that will be left above the time-varying threshold 
#
# Outputs
#
#    $mfit: fitted data (i.e. long term mean + mean annual cycle)
#
#     $ltm: long term mean
#
#$theshold: time-varying threshold
#
# Authors:
#
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 28 Feb 2006
#      Chris Ferro <c.a.t.ferro@reading.ac.uk>


x<-1:length(y)
lo<-loess(y~x,span=span)

resid<-rep(NA,length(y))
resid[!is.na(y)]<-lo$residuals

meanresid<-rep(NA,12)
for(i in 1:12){
meanresid[i]<-mean(resid[seq(i,length(y),12)],na.rm=T)
}

fitted<-rep(NA,length(y))
fitted[!is.na(y)]<-lo$fitted

mfit<-fitted+rep(meanresid,length(y)/12)

step <- sort(y[index] - mfit[index])
step <- unique(step[step > 0])
th <- mfit[index]
p <- mean(y[index] > th, na.rm = TRUE)
k <- 0
while(p > prob) {
  k <- k + 1
  th <- mfit[index] + step[k]
  p <- mean(y[index] > th, na.rm = TRUE)
}

list(threshold=th,ltm=fitted,mfit=mfit)
}