R/ews.R

In codyna: Complex Dynamic Systems

#' Detect Early Warning Signals in a Time Series
#'
#' @export
#' @param data \[`ts`, `numeric()`]\cr Univariate time series data.
#' @param method \[`character(1)`]\cr Name of the analysis method.
#'   Either `"rolling"` or `"expanding"` for rolling window and expanding
#'   window, respectively.
#' @param metrics \[`character(1)`]\cr Names of the EWS metrics to compute.
#'   The default is `"all"` computing all metrics. The available options are:
#'
#'   * `"ar1"`: The autoregressive coefficient of an AR1 model.
#'   * `"sd"`: Standard deviation.
#'   * `"skew"`: Skewness.
#'   * `"kurt"`: Kurtosis.
#'   * `"cv"`: Coefficient of variation.
#'   * `"rr"`: Return rate (`1 - ar1`).
#'   * `"all"`: All of the above.
#'
#' @param window \[`numeric(1)`]\cr Window size as a proportion of the total
#'   series length (default `0.5`).
#' @param burnin \[`numeric(1)`]\cr Burn-in period as a proportion of the total
#'   series length (default `0.1`).
#' @param demean \[`logical(1)`]\cr Should the time series be demeaned before
#'   analysis?  If `TRUE` (the default), the `"ar1"` metric will be based on an
#'   AR1 model where the mean of the observations is first subtracted.
#'   See [stats::ar.ols()] for details.
#' @param detrend \[`character(1)`]\cr Name of the detrending method to
#'   apply to the time series data before computing the metrics.
#'   The default is `"none"` for no detrending. The available options are:
#'
#'   * `"gaussian"`: Estimates a smooth curve via kernel-based regression
#'     using [stats::ksmooth()] with a Gaussian kernel which is then subtracted
#'     from the time series.
#'   * `"loess"`: Estimates a smooth curve via local polynomial regression
#'     using [stats::loess()] which is then subtracted from the time series.
#'   * `"linear"`: Fits a linear regression model via [stats::lm()] and uses
#'     the residuals for computing the metrics.
#'   * `"first-diff"`: Uses the differences between the time series and its
#'     first-order lagged values.
#'   * `"none"`: Use the original time series data.
#'
#' @param threshold \[`numeric(1)`]\cr The z-score threshold value for the
#'   expanding window method. The default is `2.0`.
#' @param consecutive \[`integer(1)`]\cr The number of times the `threshold`
#'   has to be crossed consecutively to be counted as a detection. The default
#'   is `2`.
#' @param bandwidth See [stats::ksmooth()].
#' @param span See [stats::loess()].
#' @param degree See [stats::loess()].
#' @return An object of class `ews` containing the EWS results as a `tibble`.
#' @examples
#' set.seed(123)
#' ts_data <- stats::arima.sim(list(order = c(1, 1, 0), ar = 0.6), n = 200)
#'
#' # Rolling window (default)
#' ews_roll <- detect_warnings(ts_data)
#'
#' # Expanding window
#' ews_exp <- detect_warnings(ts_data, method = "expanding")
#'
detect_warnings <- function(data, method = "rolling",
                            metrics = "all", window = 0.5, burnin = 0.1,
                            demean = TRUE, detrend = "none",
                            threshold = 2.0, consecutive = 2L,
                            bandwidth, span, degree) {
  data <- prepare_timeseries_data(data)
  values <- data$values
  time <- data$time
  method <- check_match(method, c("rolling", "expanding"))
  available_metrics <- c("ar1", "sd", "skew", "kurt", "cv", "rr")
  metrics <- metrics %m% available_metrics
  metrics <- check_match(
    metrics,
    c(available_metrics, "all"),
    several.ok = TRUE
  )
  metrics <- ifelse_("all" %in% metrics, available_metrics, metrics)
  check_range(window)
  check_range(burnin)
  check_values(consecutive, strict = TRUE)
  check_flag(demean)
  detrend <- check_match(
    detrend,
    c("none", "gaussian", "loess", "linear", "first-diff")
  )
  n <- length(values)
  window <- floor(window * n)
  burnin <- floor(burnin * n)
  bandwidth <- bandwidth %m% round(window / 2)
  span <- span %m% 0.25
  degree <- degree %m% 2
  values <- detrend_ts(values, time, detrend, window, bandwidth, span, degree)
  ifelse_(
    method == "rolling",
    rolling_ews(values, time, metrics, window, demean),
    expanding_ews(values, time, metrics, burnin, demean, threshold, consecutive)
  )
}

rolling_ews <- function(x, time, metrics, window, demean) {
  w <- window
  n <- length(x)
  m <- n - w + 1L
  idx <- seq_len(m)
  rolling_ar1 <- numeric(m)
  rolling_ar1_demean <- numeric(m)
  rolling_mean <- numeric(m)
  rolling_var <- numeric(m)
  rolling_skew <- numeric(m)
  rolling_kurt <- numeric(m)
  y <- x[1:w]
  s1 <- sum(y)
  s2 <- sum(y^2)
  s3 <- sum(y^3)
  s4 <- sum(y^4)
  s_cur <- s1 - y[1L]
  s_lag <- s1 - x[w]
  s_lag2 <- s2 - x[w]^2
  s_prod <- sum(y[-1L] * y[-w])
  mu <- s1 / w
  m2 <- (s2 - s1^2 / w) / w
  m3 <- (s3 - 3 * mu * s2 + 2 * w * mu^3) / w
  m4 <- (s4 - 4 * mu * s3 + 6 * mu^2 * s2 - 3 * w * mu^4) / w
  mu_sq <- (w - 1) * mu^2
  rolling_ar1[1L] <- ifelse_(
    demean,
    (s_prod - mu * (s_cur + s_lag) + mu_sq) / (s_lag2 - 2 * mu * s_lag + mu_sq),
    s_prod / s_lag2
  )
  rolling_mean[1L] <- mu
  rolling_var[1L] <- m2 * w / (w - 1)
  rolling_skew[1L] <- m3 / (m2^(3/2))
  rolling_kurt[1L] <- m4 / (m2^2)
  for (i in seq(2, m)) {
    x_new <- x[i + w - 1L]
    x_old <- x[i - 1L]
    s1 <- s1 + x_new - x_old
    s2 <- s2 + x_new^2 - x_old^2
    s3 <- s3 + x_new^3 - x_old^3
    s4 <- s4 + x_new^4 - x_old^4
    s_lag2 <- s_lag2 + x[i + w - 2L]^2 - x_old^2
    s_prod <- s_prod - x_old * x[i] + x_new * x[i + w - 2L]
    mu <- s1 / w
    m2 <- (s2 - s1^2 / w) / w
    m3 <- (s3 - 3 * mu * s2 + 2 * w * mu^3) / w
    m4 <- (s4 - 4 * mu * s3 + 6 * mu^2 * s2 - 3 * w * mu^4) / w
    if (demean) {
      s_lag <- s_lag + x[i + w - 2L] - x_old
      s_cur <- s_cur + x_new - x[i]
      mu_sq <- (w - 1) * mu^2
      rolling_ar1[i] <- (s_prod - mu * (s_cur + s_lag) + mu_sq) /
        (s_lag2 - 2 * mu * s_lag + mu_sq)
    } else {
      rolling_ar1[i] <- s_prod / s_lag2
    }
    rolling_mean[i] <- mu
    rolling_var[i] <- m2 * w / (w - 1)
    rolling_skew[i] <- m3 / (m2^(3/2))
    rolling_kurt[i] <- m4 / (m2^2)
  }
  rolling_sd <- sqrt(rolling_var)
  rolling_metrics <- data.frame(
    time = time[w:n],
    ar1 = rolling_ar1,
    mean = rolling_mean,
    sd = rolling_sd,
    skew = rolling_skew,
    kurt = rolling_kurt,
    cv = rolling_sd / rolling_mean,
    rr = 1 - rolling_ar1
  )
  rolling_metrics <- rolling_metrics[, c("time", metrics)]
  kendall_tau <- apply(rolling_metrics[, -1, drop = FALSE], 2, function(z) {
    stats::cor.test(
      x = idx,
      y = z,
      alternative = "two.sided",
      conf.level = 0.95,
      method = "kendall"
    )$estimate
  })
  long <- tidyr::pivot_longer(
    rolling_metrics,
    cols = !(!!rlang::sym("time")),
    names_to = "metric",
    values_to = "score"
  ) |>
    dplyr::group_by(!!rlang::sym("metric")) |>
    dplyr::mutate(std = as.numeric(scale(!!rlang::sym("score")))) |>
    dplyr::ungroup()
  structure(
    long,
    orig_values = x,
    orig_time = time,
    cor = kendall_tau,
    method = "rolling",
    class = c("ews", "tbl_df", "tbl", "data.frame")
  )
}

expanding_ews <- function(x, time, metrics, burnin, demean,
                          threshold, consecutive) {
  w <- burnin + 1L
  n <- length(x)
  m <- n - w + 1L
  idx <- seq_len(m)
  expanding_ar1 <- numeric(m)
  expanding_mean <- numeric(m)
  expanding_var <- numeric(m)
  expanding_skew <- numeric(m)
  expanding_kurt <- numeric(m)
  y <- x[1:w]
  s1 <- sum(y)
  s2 <- sum(y^2)
  s3 <- sum(y^3)
  s4 <- sum(y^4)
  s_cur <- s1 - y[1L]
  s_lag <- s1 - x[w]
  s_lag2 <- s2 - x[w]^2
  s_prod <- sum(y[-1L] * y[-w])
  mu <- s1 / w
  m2  <- (s2 - s1^2 / w) / w
  m3  <- (s3 - 3 * mu * s2 + 2 * w * mu^3) / w
  m4  <- (s4 - 4 * mu * s3 + 6 * mu^2 * s2 - 3 * w * mu^4) / w
  mu_sq <- (w - 1) * mu^2
  expanding_ar1[1] <- ifelse_(
    demean,
    (s_prod - mu * (s_cur + s_lag) + mu_sq) / (s_lag2 - 2 * mu * s_lag + mu_sq),
    s_prod / s_lag2
  )
  expanding_mean[1] <- mu
  expanding_var[1] <- m2 * w / (w - 1)
  expanding_skew[1] <- m3 / (m2^(3/2))
  expanding_kurt[1] <- m4 / (m2^2)
  for (i in seq(2, m)) {
    w <- w + 1L
    x_new <- x[w]
    s1 <- s1 + x_new
    s2 <- s2 + x_new^2
    s3 <- s3 + x_new^3
    s4 <- s4 + x_new^4
    s_lag2 <- s_lag2 + x[w - 1L]^2
    s_prod <- s_prod + x_new * x[w - 1L]
    mu <- s1 / w
    m2  <- (s2 - s1^2 / w) / w
    m3  <- (s3 - 3 * mu * s2 + 2 * w * mu^3) / w
    m4  <- (s4 - 4 * mu * s3 + 6 * mu^2 * s2 - 3 * w * mu^4) / w
    if (demean) {
      s_lag <- s_lag + x[w - 1L]
      s_cur <- s_cur + x_new
      mu_sq <- (w - 1) * mu^2
      expanding_ar1[i] <- (s_prod - mu * (s_cur + s_lag) + mu_sq) /
        (s_lag2 - 2 * mu * s_lag + mu_sq)
    } else {
      expanding_ar1[i] <- s_prod / s_lag2
    }
    expanding_mean[i] <- mu
    expanding_var[i] <- m2 * w / (w - 1)
    expanding_skew[i] <- m3 / (m2^(3/2))
    expanding_kurt[i] <- m4 / (m2^2)
  }
  expanding_sd <- sqrt(expanding_var)
  expanding_metrics <- data.frame(
    time = time[(burnin + 1):n],
    ar1 = expanding_ar1,
    mean = expanding_mean,
    sd = expanding_sd,
    skew = expanding_skew,
    kurt = expanding_kurt,
    cv = expanding_sd / expanding_mean,
    rr = 1 - expanding_ar1
  )
  expanding_metrics <- expanding_metrics[, c("time", metrics)]
  signs <- rep(1, ncol(expanding_metrics) - 1)
  names(signs) <- names(expanding_metrics[-1])
  if ("rr" %in% names(signs)) {
    signs["rr"] <- -1
  }
  long <- tidyr::pivot_longer(
    expanding_metrics,
    cols = !(!!rlang::sym("time")),
    names_to = "metric",
    values_to = "score"
  ) |>
    dplyr::group_by(!!rlang::sym("metric")) |>
    dplyr::mutate(
      z_score = expanding_z(!!rlang::sym("score")),
      detected = as.integer(
        check_lags(
          !!rlang::sym("z_score") * signs[!!rlang::sym("metric")],
          threshold,
          consecutive
        )
      )
    ) |>
    dplyr::ungroup()
  structure(
    long,
    orig_values = x,
    orig_time = time,
    threshold = threshold,
    classification = classify_ews(long),
    method = "expanding",
    class = c("ews", "tbl_df", "tbl", "data.frame")
  )
}

expanding_z <- function(x) {
  n <- length(x)
  var <- numeric(n)
  cent <- numeric(n)
  y <- x[1:2]
  s1 <- sum(y)
  s2 <- sum(y^2)
  mu <- s1 / 2
  m2 <- (s2 - s1^2 / 2) / 2
  cent[2] <- (x[2] - mu)
  var[2] <- m2 * 2
  for (i in 3:n) {
    x_new <- x[i]
    s1 <- s1 + x_new
    s2 <- s2 + x_new^2
    mu <- s1 / i
    m2  <- (s2 - s1^2 / i) / i
    cent[i] <- (x[i] - mu)
    var[i] <- m2 * i / (i - 1)
  }
  c(0, cent[-1] / sqrt(var[-1]))
}

detrend_ts <- function(values, time, method, window, bandwidth, span, degree) {
  switch(method,
    `gaussian` = {
      smoothed <- stats::ksmooth(
        x = time,
        y = values,
        kernel = "normal",
        bandwidth = bandwidth,
        x.points = time
      )$y
      values - smoothed
    },
    `loess` = {
      fit <- stats::loess(
        values ~ time,
        span = span,
        degree = degree,
        normalize = FALSE
      )
      stats::residuals(fit)
    },
    `linear` = {
      fit <- stats::lm(values ~ time)
      stats::residuals(fit)
    },
    `first-diff` = {
      c(0, diff(values))
    },
    `none` = {
      values
    }
  )
}

check_lags <- function(x, threshold, consecutive) {
  n <- length(x)
  lagged <- vapply(
    0:(consecutive - 1L),
    function(k) {
      c(rep(NA, k), x[1:(n - k)])
    },
    numeric(n)
  )
  apply(lagged, 1L, function(row) all(row > threshold, na.rm = FALSE))
}

#' Analyzes the expanding window EWS results to classify the system
#' state at each time point based on the number of metrics showing warnings.
#' This provides a qualitative assessment of system stability.
#'
#' @param x A `ews` object from the expanding window method.
#' @return A `data.frame` with time, warning count, and state classifications.
#' @noRd
classify_ews <- function(x) {
  n_metrics <- length(unique(x$metric))
  d <- x |>
    dplyr::filter(!is.na(!!rlang::sym("z_score"))) |>
    dplyr::group_by(!!rlang::sym("time")) |>
    dplyr::summarize(
      time = (!!rlang::sym("time"))[1],
      count = sum(!!rlang::sym("detected"))
    )
  if (n_metrics == 1L) {
   d$state <- cut(
      d$count,
      breaks = c(-Inf, 0, 1, Inf),
      labels = c("Stable", "Warning", "Failing"),
      right = TRUE
    )
  } else if (n_metrics <= 3) {
    d$state <- cut(
      d$count,
      breaks = c(-Inf, 0, 1, n_metrics, Inf),
      labels = c("Stable", "Vulnerable", "Warning", "Failing"),
      right = TRUE
    )
  } else {
    d$state <- cut(
      d$count,
      breaks = c(-Inf, 0, 1, 2, 0.5 * n_metrics, Inf),
      labels = c(
        "Stable", "Vulnerable", "Warning", "Critical", "Failing"
      ),
      right = TRUE
    )
  }
  d
}