R/filters_signal.R

In trendseries: Extract Trends from Time Series

#' Signal Processing Methods
#'
#' @description Internal functions for signal processing-based trend extraction
#' methods including kernel smoothing and Kalman smoothing.
#'
#' @name signal-filters
#' @keywords internal

#' Extract kernel smoother trend
#' @noRd
.extract_kernel_trend <- function(ts_data, bandwidth, kernel_type, .quiet) {
  if (!.quiet) {
    bandwidth_msg <- if (is.null(bandwidth)) "auto" else "{bandwidth}"
    cli::cli_inform(
      "Computing kernel smoother with bandwidth = {bandwidth_msg}, kernel = {kernel_type}"
    )
  }

  return(.kernel_smooth(ts_data, bandwidth, kernel_type))
}

#' Kernel smoothing implementation
#' @noRd
.kernel_smooth <- function(ts_data, bandwidth = NULL, kernel = "normal") {
  time_index <- as.numeric(stats::time(ts_data))
  values <- as.numeric(ts_data)

  # Calculate bandwidth using theoretically sound approach
  if (is.null(bandwidth)) {
    # Use Silverman's rule of thumb (optimal bandwidth)
    bandwidth <- stats::bw.nrd0(time_index)
  } else {
    # Interpret bandwidth as a multiplier of the optimal bandwidth
    # This makes smoothing scale-invariant and frequency-appropriate
    auto_bandwidth <- stats::bw.nrd0(time_index)
    bandwidth <- bandwidth * auto_bandwidth
  }

  # Use stats::ksmooth for kernel regression
  # This is equivalent to Nadaraya-Watson estimator
  smooth_result <- stats::ksmooth(
    x = time_index,
    y = values,
    kernel = kernel,
    bandwidth = bandwidth,
    x.points = time_index
  )

  trend_ts <- stats::ts(
    smooth_result$y,
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )
  return(trend_ts)
}

#' Extract Kalman smoother trend
#' @noRd
.extract_kalman_trend <- function(
  ts_data,
  measurement_noise,
  process_noise,
  .quiet,
  smoothing = NULL
) {
  if (!.quiet) {
    noise_msg <- if (is.null(measurement_noise)) {
      "auto"
    } else {
      "{measurement_noise}"
    }
    cli::cli_inform(
      "Computing Kalman smoother with measurement noise = {noise_msg}"
    )
  }

  return(.kalman_smooth(
    ts_data,
    measurement_noise,
    process_noise,
    .quiet,
    smoothing
  ))
}

#' Kalman smoothing implementation
#' @noRd
.kalman_smooth <- function(
  ts_data,
  measurement_noise = NULL,
  process_noise = NULL,
  .quiet = FALSE,
  smoothing = NULL
) {
  # Use dlm package's optimized Kalman filtering
  y <- as.numeric(ts_data)

  if (!is.null(smoothing)) {
    if (
      !is.numeric(smoothing) ||
        length(smoothing) != 1 ||
        is.na(smoothing) ||
        !is.finite(smoothing) ||
        smoothing <= 0
    ) {
      cli::cli_abort(
        "Kalman {.arg smoothing} must be one finite, positive noise ratio"
      )
    }
  }
  noises <- list(
    measurement_noise = measurement_noise,
    process_noise = process_noise
  )
  for (name in names(noises)) {
    noise <- noises[[name]]
    if (
      !is.null(noise) &&
        (!is.numeric(noise) ||
          length(noise) != 1 ||
          is.na(noise) ||
          !is.finite(noise) ||
          noise < 0)
    ) {
      cli::cli_abort(
        "Kalman {.val {name}} must be one finite, non-negative variance"
      )
    }
  }

  y_var <- stats::var(y, na.rm = TRUE)
  if (!is.null(smoothing)) {
    # Explicit variances take precedence; the ratio fills unspecified ones.
    if (is.null(process_noise)) {
      process_noise <- if (is.null(measurement_noise)) {
        y_var * 0.01
      } else {
        measurement_noise / smoothing
      }
    }
    measurement_noise <- measurement_noise %||% (smoothing * process_noise)
  } else {
    measurement_noise <- measurement_noise %||% (y_var * 0.1)
    process_noise <- process_noise %||% (y_var * 0.01)
  }

  # Build local level model (random walk + noise)
  mod <- dlm::dlmModPoly(order = 1, dV = measurement_noise, dW = process_noise)

  # Apply Kalman filter and smoother
  filtered <- dlm::dlmFilter(y, mod)
  smoothed <- dlm::dlmSmooth(filtered)

  if (is.null(smoothed$s)) {
    cli::cli_abort("Kalman smoother returned NULL states")
  }
  trend_values <- smoothed$s[-1]

  trend_ts <- stats::ts(
    trend_values,
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )
  return(trend_ts)
}