R/filters_ma.R

In trendseries: Extract Trends from Time Series

#' Moving Average Filtering Methods
#'
#' @description Internal functions for various moving average trend extraction methods.
#' These methods use C-optimized implementations (RcppRoll) for performance where possible,
#' with custom R implementations for exponential moving averages.
#'
#' @details
#' All moving average functions use C-optimized implementations (via RcppRoll) for speed.
#' NAs are preserved at the beginning of the series as expected for moving averages.
#'
#' Parameter notes:
#' - **SMA**: window parameter specifies the number of periods
#' - **WMA**: weighted MA with custom or linear weights
#' - **EWMA**: alpha parameter (0 < alpha < 1) controls smoothing strength
#' - **Triangular**: double-smoothed MA with triangular weights
#'
#' @name ma-filters
#' @keywords internal

#' Extract simple moving average trend
#' @noRd
.extract_ma_trend <- function(ts_data, window, align, .quiet) {
  # Validate window parameter
  n <- length(ts_data)
  if (window < 2) {
    cli::cli_abort("Moving average window must be at least 2, got {window}")
  }
  if (window > n) {
    cli::cli_abort(
      "Moving average window ({window}) cannot exceed series length ({n})"
    )
  }

  # Validate align parameter
  if (!align %in% c("left", "center", "right")) {
    cli::cli_abort(
      "Moving average align must be 'left', 'center', or 'right', got {.val {align}}"
    )
  }

  # Check if we need 2xN MA (even window + centered alignment)
  use_2x <- .use_2xn(window, align)

  if (!.quiet) {
    if (use_2x) {
      cli::cli_inform(
        "Computing 2x{window}-period MA (auto-adjusted for even-window centering)"
      )
    } else {
      cli::cli_inform("Computing {window}-period MA with {align} alignment")
    }
  }

  # Use appropriate implementation
  if (use_2x) {
    return(.ma_2x(ts_data, window))
  } else {
    return(.sma(ts_data, window, align))
  }
}

#' Simple Moving Average with alignment options
#' @noRd
.sma <- function(ts_data, window = 10, align = "center") {
  # Use RcppRoll for C++ optimized rolling mean
  ma_result <- RcppRoll::roll_mean(
    as.numeric(ts_data),
    n = window,
    align = align,
    fill = NA,
    na.rm = FALSE
  )

  # Convert back to ts object
  trend_ts <- stats::ts(
    ma_result,
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )
  return(trend_ts)
}

#' 2xN moving average weights for even-window centered MAs
#'
#' @description The symmetric filter behind an even centered window, with
#' half weight on the two endpoints: `c(1/2, 1, ..., 1, 1/2) / N`. The weight
#' vector has odd length `N + 1`, so `stats::filter(sides = 2)` places it on
#' the anchor exactly and pads `N / 2` positions at each end.
#'
#' Composing two `RcppRoll` passes produces the same weights but anchors them
#' one period early, because `align = "center"` puts `n / 2` observations after
#' the anchor and only `n / 2 - 1` before it.
#'
#' Shared with [roll_series()], so `stats = "mean"` and the `ma` trend method
#' agree on even centered windows.
#' @noRd
.ma_2xn <- function(v, window, na_rm = FALSE) {
  v <- as.numeric(v)
  if (length(v) < window + 1) {
    return(rep(NA_real_, length(v)))
  }
  weights <- c(0.5, rep(1, window - 1), 0.5) / window
  if (na_rm) {
    observed <- !is.na(v)
    v[!observed] <- 0
    mass <- as.numeric(stats::filter(as.numeric(observed), weights, sides = 2))
    out <- as.numeric(stats::filter(v, weights, sides = 2)) / mass
    out[is.na(mass) | mass == 0] <- NA_real_
  } else {
    out <- stats::filter(v, weights, sides = 2)
  }

  return(as.numeric(out))
}

#' 2xN Moving Average for even-window centered MAs
#' @description Implements econometrically correct centered MA for even windows.
#' This is the standard approach used in X-13ARIMA-SEATS for seasonal adjustment.
#' @noRd
.ma_2x <- function(ts_data, window) {
  trend_ts <- stats::ts(
    .ma_2xn(ts_data, window),
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )
  return(trend_ts)
}

#' Extract EWMA trend
#' @noRd
.extract_ewma_trend <- function(ts_data, window = NULL, alpha = NULL, .quiet) {
  # Validate parameters - exactly one should be provided
  if (!is.null(window) && !is.null(alpha)) {
    cli::cli_abort("Provide either 'window' or 'alpha' for EWMA, not both")
  }

  # Validate window if provided
  if (!is.null(window)) {
    n <- length(ts_data)
    if (window < 2) {
      cli::cli_abort("EWMA window must be at least 2, got {window}")
    }
    if (window > n) {
      cli::cli_abort("EWMA window ({window}) cannot exceed series length ({n})")
    }
  }

  # Validate alpha if provided
  if (!is.null(alpha)) {
    if (alpha <= 0 || alpha >= 1) {
      cli::cli_abort(
        "EWMA alpha must be between 0 and 1 (exclusive), got {alpha}"
      )
    }
  }

  if (!.quiet) {
    if (!is.null(window)) {
      cli::cli_inform("Computing EWMA with window = {window}")
    } else {
      display_alpha <- alpha %||% 0.1
      cli::cli_inform("Computing EWMA with alpha = {display_alpha}")
    }
  }

  return(.ewma(ts_data, window = window, alpha = alpha))
}

#' Exponentially Weighted Moving Average
#' @noRd
.ewma <- function(ts_data, window = NULL, alpha = NULL) {
  # A window maps to alpha as in TTR::EMA()
  if (!is.null(window)) {
    alpha <- 2 / (window + 1)
  }
  alpha <- alpha %||% 0.1

  # S_t = alpha * y_t + (1 - alpha) * S_{t-1}, starting from S_1 = y_1
  y <- as.numeric(ts_data)
  ema_result <- stats::filter(
    alpha * y,
    1 - alpha,
    method = "recursive",
    init = y[1]
  )

  # Convert back to ts object
  trend_ts <- stats::ts(
    ema_result,
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )
  return(trend_ts)
}


#' Extract WMA trend
#' @noRd
.extract_wma_trend <- function(ts_data, window, weights, align, .quiet) {
  # Validate window parameter
  n <- length(ts_data)
  if (window < 2) {
    cli::cli_abort("WMA window must be at least 2, got {window}")
  }
  if (window > n) {
    cli::cli_abort(
      "WMA window ({window}) cannot exceed series length ({n})"
    )
  }

  # Validate weights if provided
  if (!is.null(weights)) {
    if (length(weights) != window) {
      cli::cli_abort(
        "WMA weights length ({length(weights)}) must match window ({window})"
      )
    }
    if (!is.numeric(weights) || any(weights < 0)) {
      cli::cli_abort("WMA weights must be non-negative numeric values")
    }
  }

  # Validate align parameter
  if (!align %in% c("left", "center", "right")) {
    cli::cli_abort(
      "WMA align must be 'left', 'center', or 'right', got {.val {align}}"
    )
  }

  if (!.quiet) {
    weight_msg <- if (is.null(weights)) "linear weights" else "custom weights"
    cli::cli_inform(
      "Computing {window}-period weighted MA with {weight_msg}, {align} alignment"
    )
  }

  return(.wma(ts_data, window, weights, align))
}

#' Weighted Moving Average (WMA)
#' @noRd
.wma <- function(ts_data, window = 10, weights = NULL, align = "center") {
  # Default to linear weights if not provided
  if (is.null(weights)) {
    weights <- 1:window
  }

  # Normalize weights
  weights <- weights / sum(weights)

  # Use RcppRoll for C++ optimized rolling mean with weights
  wma_result <- RcppRoll::roll_mean(
    as.numeric(ts_data),
    n = window,
    weights = weights,
    align = align,
    fill = NA,
    na.rm = FALSE
  )

  # Convert back to ts object
  trend_ts <- stats::ts(
    wma_result,
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )
  return(trend_ts)
}

#' Henderson filter weights
#' @description Computes symmetric Henderson filter weights for a filter of
#' length n. The weights minimise the sum of squared third differences of
#' the trend, while exactly reproducing polynomials up to degree 3.
#'
#' The formula is: w_j ∝ [(m+1)²−j²][(m+2)²−j²][(m+3)²−j²]·(η−j²)
#' where m=(n−1)/2 and η is chosen so that Σj²wⱼ = 0.
#'
#' @references
#' Henderson, R. (1916). Note on graduation by adjusted average.
#' Transactions of the Actuarial Society of America, 17, 43–48.
#'
#' ABS (2003). A Guide to Interpreting Time Series. Australian Bureau of Statistics.
#' @noRd
.henderson_weights <- function(n) {
  if (n %% 2 == 0) {
    cli::cli_abort("Henderson filter length must be odd, got {n}")
  }
  if (n < 5) {
    cli::cli_abort("Henderson filter length must be at least 5, got {n}")
  }

  m <- (n - 1L) / 2L
  j <- seq_len(n) - (m + 1L) # -m, ..., 0, ..., m

  # Product of three quadratic factors (base kernel)
  P <- ((m + 1L)^2L - j^2L) *
    ((m + 2L)^2L - j^2L) *
    ((m + 3L)^2L - j^2L)

  # Compute eta so that Σj²wⱼ = 0 (cubic polynomial reproduction)
  sum_j2P <- sum(j^2L * P)
  sum_j4P <- sum(j^4L * P)
  eta <- sum_j4P / sum_j2P

  # Full weights
  w <- P * (eta - j^2L)
  return(w / sum(w))
}

#' Extract Henderson trend
#' @noRd
.extract_henderson_trend <- function(ts_data, window, .quiet) {
  n <- length(ts_data)

  # Default window based on frequency
  if (is.null(window)) {
    freq <- stats::frequency(ts_data)
    window <- if (freq == 4L) 9L else 13L
  }

  window <- as.integer(window)

  # Enforce odd window
  if (window %% 2L == 0L) {
    window <- window + 1L
    if (!.quiet) {
      cli::cli_inform(
        "Henderson filter requires an odd window. Adjusted to {window}."
      )
    }
  }

  if (window < 5L) {
    cli::cli_abort("Henderson filter window must be at least 5, got {window}")
  }

  if (window > n) {
    cli::cli_abort(
      "Henderson filter window ({window}) cannot exceed series length ({n})"
    )
  }

  if (!.quiet) {
    cli::cli_inform("Computing {window}-term Henderson moving average")
  }

  weights <- .henderson_weights(window)
  result <- stats::filter(as.numeric(ts_data), filter = weights, sides = 2L)

  trend_ts <- stats::ts(
    as.numeric(result),
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )
  return(trend_ts)
}

#' Extract Triangular MA trend
#' @noRd
.extract_triangular_trend <- function(ts_data, window, align, .quiet) {
  # Validate window parameter
  n <- length(ts_data)
  if (window < 3) {
    cli::cli_abort("Triangular MA window must be at least 3, got {window}")
  }
  if (window > n) {
    cli::cli_abort(
      "Triangular MA window ({window}) cannot exceed series length ({n})"
    )
  }

  # Validate align parameter
  if (!align %in% c("center", "right")) {
    cli::cli_abort(
      "Triangular MA align must be 'center' or 'right', got {.val {align}}"
    )
  }

  if (!.quiet) {
    cli::cli_inform(
      "Computing {window}-period triangular MA with {align} alignment"
    )
  }

  return(.triangular(ts_data, window, align))
}

#' Triangular Moving Average (custom implementation using stats::filter)
#' @noRd
.triangular <- function(ts_data, window = 10, align = "center") {
  # Create triangular weights
  if (window %% 2 == 1) {
    # Odd window: symmetric triangle with peak at center
    mid <- (window + 1) / 2
    weights <- c(1:mid, (mid - 1):1)
  } else {
    # Even window: two middle values are equal (flat peak)
    mid <- window / 2
    weights <- c(1:mid, mid:1)
  }

  # For right alignment, reverse weights so recent observations get higher weights
  if (align == "right") {
    weights <- rev(weights)
  }

  # Normalize weights to sum to 1
  weights <- weights / sum(weights)

  # Set sides parameter based on alignment
  sides <- if (align == "center") 2L else 1L

  # Use stats::filter for efficient convolution
  result <- stats::filter(
    as.numeric(ts_data),
    filter = weights,
    method = "convolution",
    sides = sides
  )

  # Convert back to ts object
  trend_ts <- stats::ts(
    as.numeric(result),
    start = stats::start(ts_data),
    frequency = stats::frequency(ts_data)
  )

  return(trend_ts)
}