R/dcor_serial_test.R

In dcorBSS: Distance-Correlation Based Methods for Blind Source Separation and Dependence Analysis

#' Distance Dependogram
#'
#' Computes lagwise distance covariance or distance correlation values for a
#' univariate or multivariate time series and returns an object that can be
#' plotted as a dependogram using [plot()].
#'
#' @param x A numeric vector or numeric matrix containing a time series. If
#'   \code{x} is a matrix, rows are interpreted as time points and columns as
#'   components.
#' @param transform Character string indicating whether the data should be
#'   transformed before computing the lagwise dependence measure. One of
#'   \code{"none"}, \code{"bowl"}, or \code{"biloop"}.
#' @param measure Character string specifying the lagwise dependence measure.
#'   Either \code{"dcov"} for distance covariance or \code{"dcor"} for distance
#'   correlation. The default is \code{"dcov"}.
#' @param lags Integer vector of lags to compute. If \code{NULL}, the default
#'   is \code{1:floor(n^(2/5))}.
#' @param block_size Either \code{NULL} or one positive integer. If
#'   \code{NULL}, the standard full computation is used. If positive, the
#'   blockwise computation is used with that block size.
#' @param ... Additional arguments passed to [bowl_transform()] or
#'   [biloop_transform()] when \code{transform} is not \code{"none"}.
#'
#' @return
#' An object of class \code{"sdt"} containing the lagwise dependence values in
#' \code{estimate}. It can be plotted with [plot()].
#'
#' @details
#' The function computes either \eqn{\mathrm{dCov}(X_t, X_{t+k})} or
#' \eqn{\mathrm{dCor}(X_t, X_{t+k})} for each lag \eqn{k}. The choice is
#' controlled by \code{measure}.
#' 
#' If \code{block_size = NULL}, the standard full pairwise-distance computation
#' is used at each lag. If \code{block_size} is positive, a blockwise
#' computation is used. The block size does not need to divide \code{n - k};
#' the final block may be smaller.
#'
#' If \code{transform = "bowl"} or \code{transform = "biloop"}, the
#' transformation is applied before computing the lagwise distance
#' autocorrelations. For multivariate input, \code{"biloop"} is allowed but is
#' generally not recommended.
#'
#' The returned object is compatible with [plot.sdt()], so the lagwise values
#' can be shown as a barplot or line plot in the same style as the dependograms
#' produced by [dcor_serial_test()] and [hsic_serial_test()].
#'
#' @examples
#' x <- as.numeric(arima.sim(list(ar = 0.6), n = 200))
#' fit <- dacf_curve(x, lags = 1:10)
#' plot(fit)
#'
#' y <- as.numeric(arima.sim(list(ar = 0.6), n = 200))
#' fit2 <- dacf_curve(y, lags = 1:10, block_size = 64L)
#' plot(fit2, type = "line")
#'
#' # Multivariate example
#' set.seed(1)
#' X <- cbind(
#'   as.numeric(arima.sim(list(ar = 0.5), n = 200)),
#'   as.numeric(arima.sim(list(ar = -0.3), n = 200))
#' )
#' fit3 <- dacf_curve(X, lags = 1:8)
#' plot(fit3)
#'
#' @export
dacf_curve <- function(x,
                       transform = c("none", "bowl", "biloop"),
                       measure = c("dcov", "dcor"),
                       lags = NULL,
                       block_size = NULL,
                       ...) {
  transform <- match.arg(transform)
  measure <- match.arg(measure)
  data.name <- deparse(substitute(x))

  x <- .as_matrix_input(x)

  if (!is.numeric(x)) {
    stop("'x' must be numeric.", call. = FALSE)
  }

  if (anyNA(x)) {
    stop("'x' must not contain missing values.", call. = FALSE)
  }

  if (nrow(x) < 5L) {
    stop("'x' must contain at least 5 non-missing observations.", call. = FALSE)
  }
  
  .check_block_size_or_null_dcov(block_size)
  
  obs <- .serial_dcor_statistic(
    x,
    transform = transform,
    measure = measure,
    type = "BP",
    lags = lags,
    block_size = block_size,
    ...
  )
  
  measure_name <- if (measure == "dcov") "dCov" else "dCor"
  
  out <- list(
    estimate = stats::setNames(obs$dependence, paste0("lag_", obs$lags)),
    lags = obs$lags,
    measure_name = measure_name,
    measure = measure,
    transform = transform,
    method = paste0(
      measure_name,
      " dependogram",
      if (!is.null(block_size)) " (blockwise computation)" else "",
      if (transform != "none") paste0(" with ", transform, " transformation") else ""
    ),
    data.name = data.name
  )

  class(out) <- "sdt"
  out
}


#' Test Serial Dependence Using Distance Covariance or Distance Correlation
#'
#' Test whether a univariate time series exhibits serial dependence using a
#' portmanteau-type test based on lagwise distance covariance or distance
#' correlation statistics.
#'
#' @param x A numeric vector, or a one-column numeric matrix, containing a
#'   univariate time series.
#' @param transform Character string indicating whether the data should be
#'   transformed before computing the lagwise dependence measure. One of
#'   \code{"none"}, \code{"bowl"}, or \code{"biloop"}.
#' @param measure Character string specifying the lagwise dependence measure.
#'   Either \code{"dcov"} for distance covariance or \code{"dcor"} for distance
#'   correlation. The default is \code{"dcov"}.
#' @param type Character string specifying the portmanteau statistic. Either
#'   \code{"BP"} for a Box--Pierce type statistic or \code{"FP"} for the
#'   kernel-weighted Fokianos-Pitsillou statistic.
#' @param kernel Character string specifying the kernel used when
#'   \code{type = "FP"}. Currently only \code{"bartlett"} is implemented.
#' @param B Number of permutation resamples used to approximate the null
#'   distribution.
#' @param lags Integer vector of lags to include in the test statistic. If
#'   \code{NULL}, the default is \code{1:floor(n^(2/5))}. In case of subject
#'    knowledge it is recommended not to use the default but a vector of the 
#'    form \code{1:m}, where \code{m} is the largest meaningful lag of interest 
#'    strictly smaller than the length of the time series.
#' @param seed Optional random seed.
#' @param block_size Either \code{NULL} or one positive integer. If
#'   \code{NULL}, the standard full computation is used for the lagwise
#'   distance covariance or distance correlation values. If positive, the
#'   blockwise computation is used with that block size.
#' @param ... Additional arguments passed to [bowl_transform()] or
#'   [biloop_transform()] when \code{transform} is not \code{"none"}.
#'
#' @return An object of classes \code{"sdt"} and \code{"htest"}.
#'
#' @details
#' For each lag \eqn{k}, the procedure computes a dependence measure
#' \eqn{D_k} between \eqn{X_t} and \eqn{X_{t+k}}. The measure is either
#' distance covariance or distance correlation.
#'
#' If \code{type = "BP"}, the test statistic is
#' \deqn{T_{BP} = n \sum_{k \in \mathcal{L}} D_k^2.}
#'
#' If \code{type = "FP"}, the test statistic is
#' \deqn{T_{FP} = \sum_{k \in \mathcal{L}} (n-k) K^2(k/m) D_k^2,}
#' where \eqn{K} is a kernel function and \eqn{m} is the number of selected
#' lags. Currently, the Bartlett kernel is used:
#' \deqn{
#' K(z) =
#' \begin{cases}
#' 1 - |z|, & |z| \le 1 \\
#' 0, & \text{otherwise}
#' \end{cases}
#' }
#' Hence, for \code{type = "FP"}, lags larger than \eqn{m} and the largest lag
#' of \code{1:m} get weight zero. A warning is given if selected lags other than
#' the largest lag of \code{1:m} do not contribute to the test statistic.
#'
#' The null distribution is approximated by random permutations of the observed
#' series. This targets the null hypothesis of serial independence. The
#' procedure is simple and broadly applicable, but it is computationally more
#' intensive than classical autocorrelation-based tests.
#'
#' The procedure assumes that the time series is completely observed. Missing
#' observations should therefore be handled before applying the test.
#'
#' In practical applications, a larger number of permutation resamples is
#' recommended to obtain stable and reliable p-values. The default
#' \code{B = 499} is intended for exploratory use, whereas values such as
#' \code{B = 2000} or larger are recommended for reporting final results.
#'
#' If \code{transform = "bowl"} or \code{transform = "biloop"}, the series is
#' first transformed and the distance correlation is then computed on the
#' transformed data. By default, scaling is enabled for these transformations,
#' unless the user explicitly supplies \code{do_scale = FALSE}.
#'
#' Robust scaling can be supplied through the \code{location} and \code{scale}
#' arguments. For example, for univariate data one may use the sample median
#' together with \code{\link[robustbase:Qn]{robustbase::Qn}}.
#'
#' @references
#' Fokianos, K. and Pitsillou, M (2017).
#' Consistent testing for pairwise dependence in time series.
#' \emph{Technometrics}, \strong{59}(2), 262--270
#' <doi:10.1080/00401706.2016.1156024>
#' 
#' Székely, G. J., Rizzo, M. L., and Bakirov, N. K. (2007).
#' Measuring and testing dependence by correlation of distances.
#' \emph{Annals of Statistics}, \strong{35}(6), 2769--2794.
#' <doi:10.1214/009053607000000505>.
#'
#' Székely, G. J. and Rizzo, M. L. (2009).
#' Brownian distance covariance.
#' \emph{Annals of Applied Statistics}, \strong{3}(4), 1236--1265.
#' <doi:10.1214/09-AOAS312>.
#'
#' @examples
#' x <- rnorm(200)
#' dcor_serial_test(x, B = 99)
#'
#' y <- as.numeric(arima.sim(list(ar = 0.6), n = 200))
#'
#' dcor_serial_test(y, B = 99)
#' dcor_serial_test(y, B = 99, measure = "dcor")
#' dcor_serial_test(y, B = 99, type = "FP")
#' dcor_serial_test(y, B = 99, transform = "biloop")
#'
#' \donttest{
#' # robust scaling example
#' if (requireNamespace("robustbase", quietly = TRUE)) {
#'   ymat <- matrix(y, ncol = 1)
#'   loc <- apply(ymat, 2, stats::median)
#'   sc <- apply(ymat, 2, robustbase::Qn)
#'
#'   dcor_serial_test(
#'     y,
#'     transform = "bowl",
#'     B = 99,
#'     location = loc,
#'     scale = sc
#'   )
#' }
#'
#' dcor_serial_test(y, B = 99, block_size = 64L)
#' }
#'
#' @export
dcor_serial_test <- function(x,
                             transform = c("none", "bowl", "biloop"),
                             measure = c("dcov", "dcor"),
                             type = c("FP", "BP"),
                             kernel = "bartlett",
                             B = 499L,
                             lags = NULL,
                             seed = NULL,
                             block_size = NULL,
                             ...) {
  transform <- match.arg(transform)
  measure <- match.arg(measure)
  type <- match.arg(type)
  
  if (kernel != "bartlett") {
    stop("Currently only kernel = 'bartlett' is implemented.", call. = FALSE)
  }
  .check_B(B)
  data.name = deparse(substitute(x))

  x <- .as_matrix_input(x)

  if (ncol(x) != 1L) {
    stop("'x' must be a univariate time series.", call. = FALSE)
  }
  
  if (!is.numeric(x)) {
    stop("'x' must be numeric.", call. = FALSE)
  }
  
  if (anyNA(x)) {
    stop("'x' must not contain missing values.", call. = FALSE)
  }
  
  n <- nrow(x)
  
  if (n < 5L) {
    stop("'x' must contain at least 5 non-missing observations.", call. = FALSE)
  }
  
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  obs <- .serial_dcor_statistic(
    x,
    transform = transform,
    measure = measure,
    type = type,
    kernel = kernel,
    lags = lags,
    block_size = block_size,
    ...
  )

  if (type == "FP") {
    .warn_kernel_lags(obs$lags)
  }

  perm_stats <- vapply(
    seq_len(B),
    function(b) {
      .serial_dcor_statistic(
        sample(as.numeric(x), size = n, replace = FALSE),
        transform = transform,
        measure = measure,
        type = type,
        kernel = kernel,
        lags = obs$lags,
        block_size = block_size,
        ...
      )$statistic
    },
    numeric(1L)
  )
  
  p_value <- (1 + sum(perm_stats >= obs$statistic)) / (B + 1)
  
  measure_name <- if (measure == "dcov") "dCov" else "dCor"
  
  out <- list(
    statistic = c(T = obs$statistic),
    p.value = p_value,
    method = paste0(
      "Permutation test of serial dependence based on ",
      measure_name,
      " using a ",
      type,
      " portmanteau statistic",
      if (!is.null(block_size)) " (blockwise computation)" else "",
      if (transform != "none") paste0(" with ", transform, " transformation") else ""
    ),
    data.name = data.name,
    estimate = stats::setNames(obs$dependence, paste0("lag_", obs$lags)),
    parameter = c(B = B, lags = length(obs$lags)),
    null.value = NULL,
    alternative = "serial dependence",
    lags = obs$lags,
    measure_name = measure_name,
    measure = measure,
    type = type,
    kernel = kernel,
    transform = transform
  )
  
  class(out) <- c("sdt", "htest")
  out
}


.transform_for_serial <- function(x,
                                  transform = c("none", "bowl", "biloop"),
                                  ...) {
  transform <- match.arg(transform)
  xmat <- .as_matrix_input(x)
  
  if (!is.numeric(xmat)) {
    stop("'x' must be numeric.", call. = FALSE)
  }
  
  if (transform == "none") {
    return(xmat)
  }
  
  dots <- list(...)
  
  if (transform == "bowl") {
    if (is.null(dots$do_scale)) {
      dots$do_scale <- TRUE
    }
    
    allowed <- c("alpha", "do_scale", "location", "scale")
    dots <- dots[intersect(names(dots), allowed)]
    
    args <- c(list(X = xmat), dots)
    names(args)[1L] <- "X"
    
    return(do.call(bowl_transform, args = args, quote = FALSE))
  }
  
  if (transform == "biloop") {
    if (ncol(xmat) > 1L) {
      warning(
        "'transform = \"biloop\"' is generally not recommended for multivariate input.",
        call. = FALSE
      )
    }
    
    if (is.null(dots$do_scale)) {
      dots$do_scale <- TRUE
    }
    
    allowed <- c("c1", "c2", "do_scale", "location", "scale")
    dots <- dots[intersect(names(dots), allowed)]
    
    args <- c(list(X = xmat), dots)
    names(args)[1L] <- "X"
    
    return(do.call(biloop_transform, args = args, quote = FALSE))
  }
  
  stop("Unknown transformation.", call. = FALSE)
}


.check_block_size_or_null_dcov <- function(block_size) {
  if (is.null(block_size)) {
    return(invisible(TRUE))
  }
  
  .check_block_size_dcov(block_size)
}


.serial_dcor_statistic <- function(x,
                                   transform = c("none", "bowl", "biloop"),
                                   measure = c("dcov", "dcor"),
                                   type = c("BP", "H96", "FP"),
                                   kernel = "bartlett",
                                   lags = NULL,
                                   block_size = NULL,
                                   ...) {
  transform <- match.arg(transform)
  measure <- match.arg(measure)
  type <- match.arg(type)
  
  if (kernel != "bartlett") {
    stop("Currently only kernel = 'bartlett' is implemented.", call. = FALSE)
  }
  
  x <- .as_matrix_input(x)
  
  if (!is.numeric(x)) {
    stop("'x' must be numeric.", call. = FALSE)
  }
  
  if (anyNA(x)) {
    stop("'x' must not contain missing values.", call. = FALSE)
  }
  
  n <- nrow(x)
  
  if (is.null(lags)) {
    lags <- seq_len(max(1L, floor(n^(2 / 5))))
  }
  
  lags <- as.integer(lags)
  lags <- lags[lags >= 1L & lags < n]
  
  if (length(lags) == 0L) {
    stop("No valid lags supplied.", call. = FALSE)
  }
  
  .check_block_size_or_null_dcov(block_size)
  
  x_trans <- .transform_for_serial(x, transform = transform, ...)
  
  dependence_by_lag <- vapply(
    lags,
    function(k) {
      x1 <- x_trans[seq_len(n - k), , drop = FALSE]
      x2 <- x_trans[(k + 1L):n, , drop = FALSE]
      
      if (measure == "dcov") {
        dcov_large(x1, x2, block_size = block_size)
      } else {
        dcor_large(x1, x2, block_size = block_size)
      }
    },
    numeric(1L)
  )
  
  if (type == "BP") {
    weights <- rep(1, length(lags)) * n
  } else if (type == "H96") {
    p <- length(lags)
    weights <- .serial_dcor_kernel(lags / p, kernel = kernel)^2 * n 
  } else if (type == "FP"){
    p <- length(lags)
    weights <- .serial_dcor_kernel(lags / p, kernel = kernel)^2 * (n - lags)
  } else {
    stop("Wrong type of test statistic.")
  }
  
  stat <- sum(weights * dependence_by_lag^2)
  
  list(
    statistic = unname(stat),
    dependence = dependence_by_lag,
    lags = lags,
    measure = measure,
    type = type,
    kernel = kernel,
    transform = transform,
    weights = weights
  )
}


.serial_dcor_kernel <- function(x, kernel = "bartlett") {
  if (kernel != "bartlett") {
    stop("Currently only kernel = 'bartlett' is implemented.", call. = FALSE)
  }
  
  pmax(1 - abs(x), 0)
}