R/dependence.R

In extRemes: Extreme Value Analysis

atdf <- function(x, u, lag.max=NULL, type=c("all", "rho", "rhobar"), plot=TRUE, na.action=na.fail, ...) {
   out <- list()
   out$call <- match.call()
   if(!is.numeric(x)) stop("atdf: x must be numeric.")
   if(u <= 0 || u >= 1) stop("atdf: u must be between 0 and 1 (non-inclusive).")
   type <- match.arg(type)
   out$type <- type
   out$series <- deparse(substitute(x))
   x <- na.action(as.ts(x))
   n <- length(x)
   if(is.null(lag.max)) lag.max <- floor(10*(log10(n)))
   lag.max <- min(lag.max, n - 1)
   if(lag.max < 0) stop("atdf: lag.max must be at least 0.")
   lag <- 0:lag.max
   out$lag <- lag
   if(type=="all") {
	res <- matrix(NA, lag.max+1, 2)
	colnames(res) <- c("rho", "rhobar")
   } else res <- numeric(lag.max)+NA
   xun <- sort(x)[floor(n*u)]
   for(h in lag) {
	hold <- cbind(x[(1+h):n], x[1:(n - h)])
	hold <- apply(hold,1,min) > xun
	if(type=="all" || type=="rho") {
	   if(type=="all") res[h+1,1] <- sum(hold)/(n*(1-u))
	   else res[h+1] <- sum(hold)/(n*(1-u))
	}
	if(type=="all" || type=="rhobar") {
	   if(type=="all") res[h+1,2] <- -1 + 2*log(1-u)/log(sum(hold)/n)
	   else res[h+1] <- -1 + 2*log(1-u)/log(sum(hold)/n)
	}
   } # end of for 'h' loop.
   out$atdf <- res
   class(out) <- "atdf"
   if(plot) {
	plot(out, ...)
	invisible(out)
   } else return(out)
} # end of 'atdf' function.

plot.atdf <- function(x, type=NULL, ...) {
   if(is.null(type)) type <- x$type
   if(type=="all") {
	if(x$type != "all") stop("plot.atdf: invalid type argument.")
	x1 <- x$atdf[,1]
	x2 <- x$atdf[,2]
	par(mfrow=c(2,1), mar=c(4.1, 4.1, 2.1, 1.1))
   } else if(x$type=="all") {
	if(type=="rho") x1 <- x2 <- x$atdf[,1]
	else x1 <- x2 <- x$atdf[,2]
   } else x1 <- x2 <- x$atdf
   args <- list(...)
   if(is.null(args$main)) {
      if(is.null(args$ylab) && is.null(args$xlab)) {
         if(type=="all" ) plot(x$lag, x1, type="h", ylim=c(0,1), main=paste("auto-tail dependence function\n", x$series, sep=""), ylab="rho", xlab="", ...)
	 else if(type=="chi") plot(x$lag, x1, type="h", ylim=c(0,1), main=paste("auto-tail dependence function", x$series, sep=""), ylab="rho", xlab="lag", ...)
         if(type=="all") plot(x$lag, x2, type="h", ylim=c(-1,1), ylab="rhobar", xlab="lag", ...)
         else if(type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), main=paste("auto-tail dependence function", x$series, sep=""), ylab="rhobar", xlab="lag", ...)
      } else if(is.null(args$ylab)) {
         if(type=="all" || type=="rho") plot(x$lag, x1, type="h", ylim=c(0,1), main=paste("auto-tail dependence function", x$series, sep=""), ylab="rho", ...)
         if(type=="all") plot(x$lag, x2, type="h", ylim=c(-1,1), ylab="rhobar", ...)
	 else if(type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), main=paste("auto-tail dependence function", x$series, sep=""), ylab="rhobar", ...)
      } else if(is.null(args$xlab)) {
         if(type=="all") plot(x$lag, x1, type="h", main=paste("auto-tail dependence function", x$series, sep=""), ylim=c(0,1), xlab="", ...)
	 else if(type=="rho") plot(x$lag, x1, type="h", main=paste("auto-tail dependence function", x$series, sep=""), ylim=c(0,1), xlab="lag", ...)
         if(type=="all") plot(x$lag, x2, type="h", ylim=c(-1,1), xlab="lag", ...)
	 else if(type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), main=paste("auto-tail dependence function", x$series, sep=""), xlab="lag", ...)
      } else {
         if(type=="all" || type=="rho") plot(x$lag, x1, type="h", ylim=c(0,1), main=paste("auto-tail dependence function", x$series, sep=""), ...)
         if(type=="all" || type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), main=paste("auto-tail dependence function", x$series, sep=""), ...)
      }
   } else {
      if(is.null(args$ylab) && is.null(args$xlab)) {
         if(type=="all") plot(x$lag, x1, type="h", ylim=c(0,1), ylab="rho", xlab="", ...)
	 else if(type=="rho") plot(x$lag, x1, type="h", ylim=c(0,1), ylab="rho", xlab="lag", ...)
         if(type=="all") plot(x$lag, x2, type="h", ylim=c(-1,1), ylab="rhobar", xlab="lag", ...)
         else if(type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), ylab="rhobar", xlab="lag", ...)
      } else if(is.null(args$ylab)) {
         if(type=="all" || type=="rho") plot(x$lag, x1, type="h", ylim=c(0,1), ylab="rho", ...)
         if(type=="all" || type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), ylab="rhobar", ...)
      } else if(is.null(args$xlab)) {
         if(type=="all") plot(x$lag, x1, type="h", ylim=c(0,1), xlab="", ...)
	 else if(type=="rho") plot(x$lag, x1, type="h", ylim=c(0,1), xlab="lag", ...)
         if(type=="all" || type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), xlab="lag", ...)
      } else {
         if(type=="all" || type=="rho") plot(x$lag, x1, type="h", ylim=c(0,1), ...)
         if(type=="all" || type=="rhobar") plot(x$lag, x2, type="h", ylim=c(-1,1), ...)
      }
   } # end of if else 'main' stmts.
   invisible()
} # end of 'plot.atdf' function.

taildep <- function(x, y, u, type=c("all", "chi", "chibar"), na.rm=FALSE) {
    type <- match.arg(type)
    if(missing(y)) {
	if((is.matrix(x) || is.data.frame(x) || is.matrix(x)) && dim(x)[2] == 2) y <- x[,2]
	else stop("taildep: must supply y or else x must be a two column matrix or data frame object.")
	x <- x[,1]
    }
    n <- length(x)
    if(n != length(y)) stop("taildep: y must have same length as x.")
    if(na.rm) {
	good <- !is.na(x) & !is.na(y)
	x <- x[good]
	y <- y[good]
	n <- length(x)
    }
    xun <- sort(x)[floor(n*u)]
    yun <- sort(y)[floor(n*u)]
    id <- (x > xun) & (y > yun)
    if(type=="all" || type=="chi") chi <- sum(id, na.rm=TRUE)/(n*(1-u))
    if(type=="all" || type=="chibar") chibar <- 2*log(1 - u)/log(mean(id)) - 1
    if(type=="all") {
	res <- c(chi, chibar)
	names(res) <- c("chi", "chibar")
    } else if(type=="chi") {
	res <- chi
	names(res) <- "chi"
    } else {
	res <- chibar
	names(res) <- "chibar"
    }
    return(res)
} # end of 'taildep' function.

taildep.test <- function(x, y, cthresh=-0.5, trans="relative.rank", na.action=na.fail, ...) {
    out <- list()
    out$call <- match.call()
    out$data.name <- c(deparse(substitute(x)),deparse(substitute(y)))
    out$method <- "Reiss-Thomas (13.35)"
    out$transformation <- trans
    PARAMETER <- c(cthresh, NA, NA, NA)
    names(PARAMETER) <- c("threshold", "m", "n", "exceedance rate (%)")
    if(missing(y)) {
        if((is.matrix(x) || is.data.frame(x)) && dim(x)[2] == 2) y <- x[,2]
        else stop("taildep.test: must supply y or else x must be a two column matrix or data frame object.")
        x <- x[,1]
    }
    tmp <- cbind(x,y)
    tmp <- na.action(tmp)
    n <- dim(tmp)[1]
    x <- tmp[,1]
    y <- tmp[,2]
    out$trans.args <- list(...)
    u <- do.call(trans, c(list(x=x), list(...)))
    v <- do.call(trans, c(list(x=y), list(...)))
    u <- u - 1
    v <- v - 1
    ci <- (u + v)
    out$c.full <- ci
    ci <- ci[ci > cthresh]
    m <- length(ci)
    PARAMETER[2:4] <- c(m, n, round(m/n, digits=3)*100)
    ci <- log(ci/cthresh)
    STATISTIC <- -(sum(ci,na.rm=TRUE) + m)/sqrt(m)
    names(STATISTIC) <- "statistic"
    out$statistic <- STATISTIC
    out$alternative <- "less"
    PVAL <- pnorm(STATISTIC, lower.tail=TRUE)
    out$parameter <- PARAMETER
    out$p.value <- PVAL
    class(out) <- "htest"
    return(out)
} # end of 'taildep.test' function.

relative.rank <- function(x, div="n", ...) {
    if(div=="n") n <- length(x)
    else if(is.element(div, c("n+1","n + 1", "n+ 1", "n +1"))) n <- length(x)+1
    else stop("relative.rank: div must be one of n or n+1.")
    x <- (rank(x, ...)/n)
    return(x)
} # end of 'relative.rank' funciton.

extremalindex <- function(x, threshold, method=c("intervals","runs"), run.length=1, na.action=na.fail, ...) {

    method <- tolower(method)
    method <- match.arg(method)

    xout <- x

    data.name <- deparse(substitute(x))
    dname <- as.character(substitute(x))
    if(length(dname) > 1) dname <- c(dname[2], dname[length(dname)])
    else dname <- data.name

    u <- threshold

    x <- na.action(x)
    n <- length(x)

    eid <- x > u
    N.u <- sum(eid)

    if(sum(eid)==0) {
	res <- numeric(0)
	attr(res, "theta.msg") <- "extremalindex: No values of x exceed the threshold."
	attr(res, "method") <- "none"
	class(res) <- "extremalindex"
	return(res)
    }

    if(method=="intervals") {

	T.u <- diff((1:n)[eid])

	if(!any(T.u > 2)) {
            hold <- colSums(cbind(T.u, T.u^2))
	    theta.type <- "theta.hat"
	} else {
	    hold <- colSums(cbind(T.u - 1, (T.u - 1) * (T.u - 2)))
	    theta.type <- "theta.tilde"
	}

	theta <- 2 * (hold[1])^2/((N.u - 1) * hold[2])

	if(theta > 1) {
	    theta <- 1
	    theta.msg <- "Extremal index estimate > 1.  Re-set to 1."
	} else theta.msg <- "Valid extremal index estimated"

	K <- ifelse(round(theta * N.u, digits=0) != theta * N.u, ceiling(theta * N.u), theta * N.u)

	o <- order(T.u, na.last=TRUE, decreasing=TRUE)
	oT.u <- T.u[o]
	r <- oT.u[K]

	##
	## Handle ties
	## 
	T.u2 <- c(diff(oT.u), 0)
	ind0 <- T.u2 == 0
	ind1 <- !ind0
	if((K > 1) && (T.u2[K - 1] == 0)) {
	    ind2.0 <- (1:(N.u - 1))[ind0]
	    ind2.1 <- (1:(N.u - 1))[ind1]
	    ind3 <- (ind2.0 < K) & (ind2.0 > max(ind2.1[ind2.1 < K]))
	    K <- min(ind2.0[ind3])
	    r <- oT.u[K]
	}

	res <- c(theta, K, r)
	attr(res, "theta.type") <- theta.type

    } else if(method == "runs") {

	tmp <- decluster(x=x, threshold=threshold, method="runs", r=run.length, ...)
	nc <- length(unique(attributes(tmp)$clusters))
	N.r <- sum(c(tmp) > threshold)
	res <- c(N.r/N.u, nc, run.length)
	attr(res, "cluster") <- attributes(tmp)$cluster

    } else stop("extremalindex: sorry, currently only methods runs and intervals available.")

    names(res) <- c("extremal.index", "number.of.clusters", "run.length")

    attr(res, "data") <- xout
    attr(res, "method") <- method
    attr(res, "data.name") <- dname
    attr(res, "data.call") <- data.name
    attr(res, "call") <- match.call()
    attr(res, "na.action") <- na.action
    attr(res, "threshold") <- u
    class(res) <- "extremalindex"

    return(res)

} # end of 'extremalindex' function.

ci.extremalindex <- function(x, alpha=0.05, R=502, return.samples=FALSE, ...) {

    tmp <- attributes(x)

    r <- x["run.length"]
    nc <- x["number.of.clusters"]

    y <- datagrabber(x)
    n <- length(y)
    u <- tmp$threshold
    eid <- y > u
    Nu <- sum(eid)
    s <- (1:n)[eid]
    Tu <- diff(s)

    if(tmp$method=="runs") cluster <- tmp$cluster
    else if(tmp$method=="intervals") {
	cluster <- rep(1, Nu)
        if(Nu > 1) cluster[2:Nu] <- 1 + cumsum(Tu > r)
    }

    # Set of intercluster times.
    # ict <- c(unlist(aggregate(c(Tu,NA), by=list(cluster), function(x) return(x[1]))$x))
    # ict <- ict[!is.na(ict)]

    # Data by cluster.
    csizes <- tabulate(cluster)
    csizes <- csizes[!is.na(csizes)]
    csizes <- csizes[csizes != 0]
    mcs <- max(csizes)
    Y <- Thold <- numeric(mcs * nc) + NA

    ind <- unlist(c(apply(cbind(csizes, 0:(nc - 1)), 1, function(x, sz) 1:x[1] + x[2] * sz, sz=mcs)))
    Y[ind] <- y[eid]
    Y <- matrix(Y, mcs, nc)
    Thold[ind] <- c(1, Tu)
    Thold <- matrix(Thold, mcs, nc)

    sam <- matrix(NA, R, 3)

    for(i in 1:R) {

	ic <- sample(1:nc, replace=TRUE)
	Tub <- Thold[,ic]
	if(nc > 1) Tub[1,2:nc] <- sample(Thold[1, 2:nc], replace=TRUE)
	# if(nc > 1) Tub[1,2:nc] <- sample(ict, replace=TRUE)

	Tub <- Tub[!is.na(Tub)]
	sb <- cumsum(Tub)

	Yb.star <- Y[,ic]
	Yb <- rep(u, max(sb))
	Yb[sb] <- c(Yb.star[!is.na(Yb.star)])

	sam[i,] <- extremalindex(x=Yb, threshold=u, method=tmp$method, run.length=r)

    } # end of for 'i' loop.

    if(is.null(colnames(sam))) colnames(sam) <- c("extremal.index", "number.of.clusters", "run.length")
    if(tmp$method=="runs") sam <- sam[,-3]

    if(return.samples) return(sam)

    bds <- t(apply(sam, 2, quantile, probs=c(alpha/2, 1 - alpha/2)))

    conf.level <- paste((1 - alpha) * 100, "%", sep="")
    if(tmp$method=="runs") out <- cbind(bds[,1], x[1:2], bds[,2])
    else out <- cbind(bds[,1], x, bds[,2])
    colnames(out) <- c(paste(conf.level, " lower CI", sep=""), "Estimate", paste(conf.level, " upper CI", sep=""))

    attr(out, "data.name") <- tmp$call
    attr(out, "method") <- "Ferro-Segers Bootstrap"
    attr(out, "ei.method") <- tmp$method
    if(tmp$method=="runs") {
	names(r) <- "run.length"
	attr(out, "parameter") <- r
    }
    attr(out, "conf.level") <- conf.level
    class(out) <- "ci"
    return(out)
} # end of 'ci.extremalindex' function.

print.extremalindex <- function(x, ...) {
    tmp <- attributes(x)
    if(tmp$method=="intervals") {
	cat("\n", "Intervals Method Estimator for the Extremal Index\n")
	print(tmp$theta.msg)
	if(tmp$theta.type=="theta.hat") cat("\n", "theta.hat used because there are no inter-exceedance times > 2.\n") 
	else if(tmp$theta.type=="theta.tilde") cat("\n", "theta.tilde used because there exist inter-exceedance times > 2.\n")
    } else if(tmp$method=="runs") {
	cat("\n", "Runs Estimator for the Extremal Index\n")
    } else if(tmp$method=="none") {
	print(tmp$theta.msg)
    } else stop("print: sorry, currently no print method for anything but intervals or runs estimation methods.")
    y <- x
    attributes(y) <- NULL
    names(y) <- tmp$names
    print(y)
    invisible()
} # end of 'print.extremalindex' function.

decluster <- function(x, threshold, ...) {
    UseMethod("decluster", x)
} # end of 'decluster' function.

# decluster.formula <- function(x, threshold, ..., method=c("runs","intervals"), clusterfun="max") {
#    method <- tolower(method)
#    method <- match.arg(method)
#    class(x) <- c(method, class(x))
#    UseMethod("decluster", x)
# }

decluster.default <- function(x, threshold, ..., method=c("runs","intervals"), clusterfun="max") {

    method <- tolower(method)
    method <- match.arg(method)
    class(x) <- c(method, class(x))
    UseMethod("decluster", x)

} # end of 'decluster.default' function.

decluster.data.frame <- function(x, threshold, ..., which.cols, method=c("runs","intervals"), clusterfun="max") {

    if(missing(which.cols)) stop("decluster: Must provide a which.cols argument for data frames.")
    if(!is.element(length(which.cols), 1:2)) stop("decluster: which.cols must have length one or two.")

    y <- x[,which.cols[1]]
    if(length(which.cols) == 2) groups <- x[,which.cols[2]]
    else groups <- NULL

    out <- decluster(x=y, threshold=threshold, method=method, clusterfun=clusterfun, groups=groups, ...)
    attr(out, "data.name") <- c(deparse(substitute(x)), which.cols)
    if(is.null(groups)) attr(out, "groups") <- "NULL"

    return(out)

} # end of 'decluster.data.frame' function.

decluster.runs <- function(x, threshold, ..., data, r=1, clusterfun="max", groups=NULL, replace.with, na.action=na.fail) {

    clusterfun <- match.fun(clusterfun)
    if(missing(replace.with)) replace.with <- threshold

    xout <- x

    dname <- deparse(substitute(x))
    thname <- deparse(substitute(threshold))

#     if(!missing(data)) {
# 	# if(missing(data)) stop("decluster: must supply data argument if x is a formula.")
# 	datc <- colnames(data)
# 	if(is.element(dname, datc)) {
# 	    id <- dname == datc
# 	    x.fun <- ifelse(dname == "substitute(x)", deparse(x), dname)
# 	    x <- c(data[,id])
#             x.fun <- formula(paste(x.fun, "~ 1"))
# 	} else if(is.formula(x)) {
# 	    x.fun <- x
# 	    x.str <- unlist(strsplit(as.character(x.fun), split="~"))
# 	    x.str <- x.str[x.str != ""]
# 	    ns <- length(x.str)
# 	    if(x.str[ns] != 1) stop("decluster: invalid x argument.")
# 	    x <- eval(parse(text=x.str[1]), envir=data)
#         } else x.fun <- NULL
# 
# 	if(is.element(thname, datc)) {
# 	    id <- thname == datc
# 	    th.fun <- ifelse(thname == "substitute(theshold)", deparse(threshold), thname)
# 	    threshold <- c(data[,id])
# 	    th.fun <- formula(paste(th.fun, "~ 1"))
# 	} else if(is.formula(threshold)) {
# 	    th.fun <- threshold
# 	    th.str <- unlist(strsplit(as.character(th.fun), split="~"))
# 	    th.str <- th.str[th.str != ""]
# 	    ns <- length(th.str)
# 	    if(th.str[ns] != 1) stop("decluster: invalid threshold argument.")
# 	    threshold <- eval(parse(text=th.str[1]), envir=data)
#         }
#     } else {
# 	if(is.formula(x) || is.formula(threshold)) stop("decluster: data argument missing.")
# 	x.fun <- th.fun <- NULL
#     }

    if(length(threshold==1)) x <- na.action(x)
    else {
	look <- na.action(cbind(x, threshold))
	x <- look[,1]
	threshold <- look[,2]
    }
    n <- length(x)

    eid <- x > threshold
    xu <- x[eid]

    Nu <- sum(eid)
    s <- (1:n)[eid]
    cluster <- rep(1, Nu)

    if(is.null(groups)) {

	res <- numeric(n) + replace.with
	res[!eid] <- x[!eid]

	T.u <- diff(s)
	if(Nu > 1) {

	    cluster[2:Nu] <- 1 + cumsum(T.u > r)

	    # For some unknown reason, the next command didn't work, but the subsequent one does.
	    # cres <- c(aggregate(xu, by=list(cluster), clusterfun, ...)$x)[cluster]
	    cres <- c(aggregate(xu, by=list(cluster), clusterfun)$x)[cluster]
	    cf.id <- xu == cres

	    posfun <- function(ind) {
	        m <- length(ind)
	        res <- logical(m)
	        if(any(ind)) {
	        	ind2 <- min((1:m)[ind])
	        	res[ind2] <- TRUE
	        } else res[1] <- TRUE
	        return(res)
	    } # end of internal 'posfun' function.
	    pos <- ((1:n)[eid])[unlist(c(t(aggregate(cf.id, by=list(cluster), posfun)$x)))]
	    res[pos] <- x[pos]
	} else if(any(eid)) res[eid] <- x[eid]
	
    } else {

	b <- groups
	if(length(threshold) == 1) {

	    res <- aggregate(x, by=list(b), FUN=decluster, threshold=threshold, method="runs", r=r, clusterfun=clusterfun, groups=NULL)
	    res <- unlist(c(t(res$x)))

	} else {

	    dcfun <- function(dat, r, clusterfun) {
		xdat <- dat[,1]
		th <- dat[,2]
		return(decluster(x=xdat, threshold=th, method="runs", r=r, clusterfun=clusterfun))
	    } # end of internal 'dcfun' function.

	    hold <- list()
	    ub <- unique(b)
	    nb <- length(ub)
	    for(i in 1:nb) hold[[i]] <- cbind(x[b==ub[i]], threshold[b==ub[i]])
	    res <- lapply(hold, FUN=dcfun, r=r, clusterfun=cluster)
	    res <- c(unlist(res))
	} # end of if else 'threshold' is constant or not stmts.

	T.u <- diff(s)
	k <- c(aggregate(groups, by=list(groups), length)$x)
	bl <- rep(1:length(unique(groups)), k)[eid] - 1
        if(Nu > 1) cluster[2:Nu] <- 1 + cumsum(T.u > r) + bl[2:Nu]

    } # end of if else 'groups' are NULL stmts.

    attr(res, "call") <- match.call()
    # attr(res, "x.fun") <- x.fun
    # attr(res, "th.fun") <- th.fun
    attr(res, "data") <- xout
    attr(res, "data.name") <- dname
    attr(res, "decluster.function") <- clusterfun
    attr(res, "method") <- "runs"
    attr(res, "threshold") <- threshold
    attr(res, "groups.name") <- deparse(substitute(groups))
    attr(res, "groups") <- groups
    attr(res, "run.length") <- r
    attr(res, "na.action") <- na.action
    attr(res, "clusters") <- cluster
    class(res) <- "declustered" 

    return(res)

} # end of 'decluster.runs' function.

decluster.intervals <- function(x, threshold, ..., clusterfun="max", groups=NULL, replace.with, na.action=na.fail) {

    xout <- x

#     if(is.formula(x)) {
#         if(missing(data)) stop("decluster: must supply data argument if x is a formula.")
#         x.fun <- x
# 	x.fun <- ifelse(deparse(substitute(x)) == "substitute(x)", deparse(x), deparse(substitute(x)))
#         # x.fun <- formula(paste(x.fun, "~ 1"))
#         x <- model.response(model.matrix(x.fun, data=data))
#     } else x.fun <- NULL
# 
#     if(is.formula(threshold)) {
#         if(missing(data)) stop("decluster: must supply data argument if threshold is a formula.")
#         th.fun <- threshold
# 	th.fun <- ifelse(deparse(substitute(threshold)) == "substitute(threshold)", deparse(threshold), deparse(substitute(threshold)))
#         # th.fun <- formula(paste(th.fun, "~ 1"))
#         threshold <- model.response(model.matrix(th.fun, data=data))
#     } else th.fun <- NULL

    th <- extremalindex(x=x, threshold=threshold, na.action=na.action)
    if(th[1] >= 1) r <- 0
    else r <- th[3]

    if(missing(replace.with)) tmp <- decluster(x=x, threshold=threshold, clusterfun=clusterfun, groups=groups, na.action=na.action)
    else tmp <- decluster(x=x, threshold=threshold, clusterfun=clusterfun, groups=groups, replace.with=replace.with, na.action=na.action)

    hold <- attributes(tmp)
    res <- c(tmp)

    attr(res, "call") <- match.call()
    # attr(res, "x.fun") <- x.fun
    # attr(res, "th.fun") <- th.fun
    attr(res, "data") <- xout
    attr(res, "data.name") <- deparse(substitute(x))
    attr(res, "decluster.function") <- clusterfun
    attr(res, "method") <- "intervals"
    attr(res, "threshold") <- threshold
    attr(res, "groups.name") <- deparse(substitute(groups))
    attr(res, "groups") <- groups
    attr(res, "run.length") <- r
    attr(res, "na.action") <- na.action
    attr(res, "clusters") <- hold$clusters
    class(res) <- "declustered"

    return(res)

} # end of 'decluster.intervals' function.


plot.declustered <- function(x, which.plot=c("scatter", "atdf"), qu=0.85, xlab=NULL, ylab=NULL, main=NULL, col="gray", ...) {

    which.plot <- tolower(which.plot)
    which.plot <- match.arg(which.plot)

    tmp <- attributes(x)
     if(length(tmp$data.name) == 1) {
 	# x0 <- try(get(tmp$data.name), silent=TRUE)
 	# if(class(x0) == "try-error") x0 <- NULL
	x0 <- tmp$data
     } else x0 <- datagrabber(x)

    if(!is.null(x0)) {
        if(is.data.frame(x0) || is.matrix(x0)) {
	    bl <- x0[,2]
	    x0 <- x0[,1]
        } else bl <- NULL
    } else bl <- NULL

    x <- c(x)
    # if(tmp$data.name == "x") warning("plot.declustered: attempt to plot original data failed because original data has same name (x) as declustered argument.")
    # x0 <- tmp$na.action(x0)
    if(!is.null(x0)) n <- length(x0)
    else n <- length(x)
    m <- length(x)
 
    if(is.null(xlab)) xlab <- "index"
    if(is.null(ylab)) {
	ylab <- tmp$data.name
	if(length(ylab) == 3) ylab <- paste(ylab[1], ": ", ylab[2], " grouped by ", ylab[3], sep="")
	else if(length(ylab) == 2) ylab <- paste(ylab[1], ": ", ylab[2], sep="")
    }
    if(is.null(main) && which.plot == "scatter") main <- deparse(tmp$call)
    else main <- ""

    # if(which.plot == "primary") par(mfrow=c(5,1), oma=c(0,0,2,0), mar=c(5.1, 4, 0.1, 2.1))
    # else
    if(which.plot == "atdf") {
	if(!is.null(x0)) par(mfrow=c(2,2), oma=c(0,0,2,0), mar=c(4.1, 4.1, 2.1, 2.1))
	else par(mfrow=c(2,1), oma=c(0,0,2,0), mar=c(4.1, 4.1, 2.1, 2.1))
    }

    if(is.element(which.plot, c("primary", "scatter"))) {
        if(!is.null(x0)) plot(x0, xlab=xlab, ylab=ylab, main=main, col=col, ...)
	else plot(x, xlab=xlab, ylab=ylab, main=main, col=col, ...)
        u <- tmp$threshold
        if(length(u) == 1 || all(u == u[1])) {
	    u <- u[1]
	    abline(h=u, lty=2)
        } else if(length(u) == n) lines(1:n, u, lty=2)
        else lines(1:length(u), u, lty=2)
        if(tmp$groups.name != "NULL") {
	    if(is.null(bl)) bl <- tmp$groups
	    if( !is(bl, "try-error" ) ) # warning("plot.declustered: unable to get groups for clusters.")
	    {
		m2 <- length(unique(bl))
	        m3 <- length(bl)
	        k <- c(aggregate(bl, by=list(bl), length)$x)
	        bl <- rep(1:m2, k)
	        bl2 <- (1:m3)[diff(bl) == 1]
	        abline(v=bl2, col="blue")
	    }
        }
        points(1:m, x)
    }

    if(is.element(which.plot, c("primary", "atdf"))) {
	if(!is.null(x0)) hold0 <- atdf(x0, u=qu, plot=FALSE, na.action=tmp$na.action)
	hold1 <- atdf(x, u=qu, plot=FALSE, na.action=tmp$na.action)
	if(which.plot=="atdf") {
	    if(!is.null(x0)) plot(hold0, type="rho", ylab="rho (original data)", xlab="", main=paste(ylab[1], " (original)\n", "auto-tail dependence function", sep=""))
	    plot(hold1, type="rho", ylab="rho (declustered data)", xlab="", main=paste(ylab[1], " (declustered)\n", "auto-tail dependence function", sep=""))
 	} else {
	    if(!is.null(x0)) plot(hold0, type="rho", ylab="rho (original data)", xlab="", main="")
            plot(hold1, type="rho", ylab="rho (declustered data)", xlab="", main="")
	}
	if(which.plot == "primary" && !is.null(x0)) plot(hold0, type="rhobar", ylab="rhobar (original data)", xlab="", main="")
	else if(!is.null(x0)) plot(hold0, type="rhobar", ylab="rhobar (original data)", main="")
	plot(hold1, type="rhobar", ylab="rhobar (declustered data)", main="")
    }

    if(is.element(which.plot, c("primary","atdf"))) mtext(tmp$call, line=0.5, outer=TRUE)
    invisible()
} # end of 'plot.declustered' function.

print.declustered <- function(x, ...) {
    tmp <- attributes(x)
    cat("\n", tmp$data.name, " declustered via", tmp$method, " declustering.\n")
    hold <- extremalindex(x, threshold=tmp$threshold)
    cat("\n", "Estimated extremal index (intervals estimate) = ", hold[1], "\n")
    cat("\n", "Number of clusters = ", length(unique(tmp$clusters)), "\n")
    cat("\n", "Run length = ", tmp$run.length, "\n")
    invisible()
} # end of 'print.declustered' function.