R/ch_ews.R

In earlywarnings: Early Warning Signals for Critical Transitions in Time Series

#' Conditional Heteroskedasticity
#'
#' \code{\link{ch_ews}} is used to estimate changes in conditional heteroskedasticity within rolling windows along a timeseries
#'
#' @param timeseries a numeric vector of the observed timeseries values or a numeric matrix where the first column represents the time index and the second the observed timeseries values. Use vectors/matrices with headings.
#' @param winsize is length of the rolling window expressed as percentage of the timeseries length (must be numeric between 0 and 100). Default is 10\%.
#' @param alpha is the significance threshold (must be numeric). Default is 0.1.
#' @param optim logical. If TRUE an autoregressive model is fit to the data within the rolling window using AIC optimization. Otherwise an autoregressive model of specific order \code{lags} is selected.
#' @param lags is a parameter that determines the specific order of an autoregressive model to fit the data. Default is 4.
#' @param logtransform logical. If TRUE data are logtransformed prior to analysis as log(X+1). Default is FALSE.
#' @param interpolate logical. If TRUE linear interpolation is applied to produce a timeseries of equal length as the original. Default is FALSE (assumes there are no gaps in the timeseries). 
#' 
#' @return \code{\link{ch_ews}} returns a matrix that contains:
#'   time the time index.
#'   r.squared the R2 values of the regressed residuals.
#'   critical.value the chi-square critical value based on the desired \code{alpha} level for 1 degree of freedom divided by the number of residuals used in the regression.
#'   test.result logical. It indicates whether conditional heteroskedasticity was significant.
#'   ar.fit.order the order of the specified autoregressive model- only informative if \code{optim} FALSE was selected.
#'
#' In addition, \code{\link{ch_ews}} plots the original timeseries and the R2 where the level of significance is also indicated.
#'  
#' @export
#' 
#' @author T. Cline, modified by V. Dakos
#' @references Seekell, D. A., et al (2011). 'Conditional heteroscedasticity as a leading indicator of ecological regime shifts.' \emph{American Naturalist} 178(4): 442-451
#' 
#' Dakos, V., et al (2012).'Methods for Detecting Early Warnings of Critical Transitions in Time Series Illustrated Using Simulated Ecological Data.' \emph{PLoS ONE} 7(7): e41010. doi:10.1371/journal.pone.0041010 
#' @seealso 
#' \code{\link{generic_ews}}; \code{\link{ddjnonparam_ews}}; \code{\link{bdstest_ews}}; \code{\link{sensitivity_ews}}; \code{\link{surrogates_ews}}; \code{\link{ch_ews}}; \code{movpotential_ews}; \code{livpotential_ews} 
# ; \code{\link{timeVAR_ews}}; \code{\link{thresholdAR_ews}}
#' @examples 
#' data(foldbif)
#' out=ch_ews(foldbif, winsize=50, alpha=0.05, optim=TRUE, lags)
#' @keywords early-warning
# Author: Timothy Cline, October 25, 2011.  Modified by: Vasilis Dakos, January
# 3, 2012.

ch_ews <- function(timeseries, winsize = 10, alpha = 0.1, optim = TRUE, lags = 4, 
    logtransform = FALSE, interpolate = FALSE) {
    
    timeseries <- ts(timeseries)  #strict data-types the input data as tseries object for use in later steps
    if (dim(timeseries)[2] == 1) {
        ts.in = timeseries
        timeindex = 1:dim(timeseries)[1]
    } else if (dim(timeseries)[2] == 2) {
        ts.in = timeseries[, 2]
        timeindex = timeseries[, 1]
    } else {
        warning("not right format of timeseries input")
    }
    
    # Interpolation
    if (interpolate) {
        YY <- approx(timeindex, ts.in, n = length(ts.in), method = "linear")
        ts.in <- YY$y
    }
    
    # Log-transformation
    if (logtransform) {
        ts.in <- log(ts.in + 1)
    }
    
    winSize = round(winsize * length(ts.in)/100)
    sto <- matrix(nrow = (length(ts.in) - (winSize - 1)), ncol = 5)  # creates a matrix to store output
    
    count <- 1  #place holder for writing to the matrix
    for (i2 in winSize:length(ts.in)) {
        # loop to iterate through the model values by window lengths of the input value
        
        # the next line applys the autoregressive model optimized using AIC then we omit
        # the first data point(s) which come back as NA and square the residuals
        if (optim == TRUE) {
            arm <- ar(ts.in[(i2 - (winSize - 1)):i2], method = "ols")
        } else {
            arm <- ar(ts.in[(i2 - (winSize - 1)):i2], aic = FALSE, order.max = lags, 
                method = "ols")
        }
        resid1 <- na.omit(arm$resid)^2
        
        l1 <- length(resid1)  # stores the number of residuals for many uses later
        lm1 <- lm(resid1[2:l1] ~ resid1[1:(l1 - 1)])  #calculates simple OLS model of describing t+1 by t
        
        # calculates the critical value: Chi-squared critical value using desired alpha
        # level and 1 degree of freedom / number of residuals used in regression
        critical <- qchisq((1 - alpha), df = 1)/(length(resid1) - 1)
        
        sto[count, 1] <- timeindex[i2]  # stores a time component
        sto[count, 2] <- summary(lm1)$r.squared  # stores the r.squared for this window
        sto[count, 3] <- critical  # stores the critical value
        
        # the next flow control group stores a simple 1 for significant test or 0 for
        # non-significant test
        if (summary(lm1)$r.squared > critical) {
            sto[count, 4] <- 1
        } else {
            sto[count, 4] <- 0
        }
        sto[count, 5] <- arm$order
        count <- count + 1  # increment the place holder
    }
    
    sto <- data.frame(sto)  # data types the matrix as a data frame
    colnames(sto) <- c("time", "r.squared", "critical.value", "test.result", "ar.fit.order")  # applies column names to the data frame
    
    # This next series of flow control statements will adjust the max and minimum
    # values to yield prettier plots In some circumstances it is possible for all
    # values to be far greater or far less than the critical value; in all cases we
    # want the critical line ploted on the figure
    if (max(sto$r.squared) < critical) {
        maxY <- critical + 0.02
        minY <- critical - 0.02
    } else if (min(sto$r.squared >= critical)) {
        minY <- critical - 0.02
        maxY <- critical + 0.02
    } else {
        maxY <- max(sto$r.squared)
        minY <- min(sto$r.squared)
    }
    
    # this creates a very simple plot that is well fitted to the data.  it also plots
    # the critical value line
    
    par(mar = (c(0, 4, 0, 1) + 0), oma = c(5, 1, 2, 1), mfrow = c(2, 1))
    plot(timeindex, ts.in, type = "l", ylab = "data", xlab = "", cex.axis = 0.8, 
        cex.lab = 0.8, xaxt = "n", las = 1, xlim = c(timeindex[1], timeindex[length(timeindex)]))
    plot(timeindex[winSize:length(timeindex)], sto$r.squared, ylab = expression(paste("R^2")), 
        xlab = "time", type = "b", cex = 0.5, cex.lab = 0.8, cex.axis = 0.8, las = 1, 
        ylim = c(min(sto$r.squared), max(sto$r.squared)), xlim = c(timeindex[1], 
            timeindex[length(timeindex)]))
    legend("topleft", "conditional heteroskedasticity", bty = "n")
    abline(h = sto$critical, lwd = 0.5, lty = 2, col = 2)
    mtext("time", side = 1, line = 2, cex = 0.8)  #outer=TRUE print on the outer margin
    
    return(sto)
}