xindex <- function(y,u,threshold = FALSE) {

  # Description:
  #
  #   Evaluates the intervals estimator for the extremal index, an index for
  #   time clusters, for a given time series and threshold.
  #
  # Usage:
  #
  #   xindex(y,u,threshold=FALSE)
  #
  # Arguments:
  #
  #   y: a time series for which the extremal index will be computed
  #
  #   u: Threshold. Can be either a quantile (e.g. 0.8) or a threshold value
  #      (see definition of theshold below).
  #
  #   threshold: Logical. If FALSE (the default) u must be a quantile value.
  #              If TRUE, then u must be a threshold value.
  #
  # Value:
  # 
  #   Estimate of the extremal index.
  #
  # Details:
  #
  # Warning:
  #
  #   The estimator is not constrained to lie in [0,1] and a
  #   default value of 1 is returned if there are fewer than two
  #   extreme values.
  #
  # Authors:
  #
  #   Christopher Ferro <c.a.t.ferro@reading.ac.uk> 21 Dec 2005
  #   Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 
  # 
  # References:
  #
  #   Ferro CAT & Segers J (2003) Inference for clusters of extreme
  #   values. Journal of the Royal Statistical Society B 65, 545-556.
  #
  # See Also:
  #
  # Examples:
  #
  #   y<-rnorm(1000,25,2)
  #   xindex(y,0.9)
  #   xindex(y,28,threshold=TRUE)


 # generates a logical vector indicating which positions correspond to
  # extreme values.
  if(threshold) z <- (y >= u) 
  else z <- (y >= quantile(y, u,na.rm=T))

  if(sum(z,na.rm=T) <= 1) {
    warning("estimator undefined: too few exceedances")
    return(1)
  } else {
    nz <- length(z)          # length of sequence
    s <- c(1:nz)[z]          # exceedance times
    t <- diff(s)             # interexceedance times
    if(max(t,na.rm=T) <= 2) {
      t1 <- mean(t,na.rm=T)
      t2 <- mean(t^2,na.rm=T)
    } else {
      t1 <- mean(t-1,na.rm=T)
      t2 <- mean((t-1)*(t-2),na.rm=T)
    }
  }
  min(1,2*(t1^2)/t2)  
}