R/extract_trends.R

In trendseries: Extract Trends from Time Series

#' Extract trends from time series objects
#'
#' @description
#' Extract trend components from time series objects using various econometric methods.
#' Designed for monthly and quarterly economic data analysis. Returns trend components
#' as time series objects or a list of time series.
#'
#' @param ts_data A time series object (`ts`, `xts`, or `zoo`) or any object
#'   convertible via tsbox. Missing values inside the observed span are
#'   rejected: methods disagree on what to do with them, and several fail
#'   silently. Impute them first. Leading and trailing missing values are
#'   excluded from estimation and returned as `NA`.
#' @param methods Character vector of trend methods.
#'   Options: `"hp"`, `"bk"`, `"cf"`, `"ma"`, `"stl"`, `"loess"`, `"spline"`, `"poly"`,
#'   `"bn"`, `"ucm"`, `"hamilton"`, `"spencer"`, `"henderson"`, `"ewma"`, `"wma"`,
#'   `"triangular"`, `"kernel"`, `"kalman"`, `"median"`,
#'   `"gaussian"`.
#'   Default is `"stl"`.
#' @param window Unified window/period parameter for moving
#'   average methods (ma, wma, triangular, stl, ewma, median, gaussian,
#'   henderson). Must be positive.
#'   If NULL, uses frequency-appropriate defaults. For EWMA, the window is
#'   converted to the smoothing factor via `alpha = 2 / (window + 1)`. Cannot be
#'   used simultaneously with `smoothing` for EWMA method.
#'   For `ma`, `median`, and `henderson` methods, a numeric vector is accepted
#'   (e.g., `c(9, 13, 23)`), which runs the method once per window value and returns
#'   a named list with keys like `henderson_9`, `henderson_13`, `henderson_23`.
#'   Other methods ignore extra values (with a warning).
#' @param smoothing Unified smoothing parameter for smoothing
#'   methods (hp, loess, spline, ewma, kernel, kalman).
#'   For hp: a value above 1 is lambda itself; a value of 1 or less is a
#'   fraction of the default lambda for the frequency, `1600 * (frequency / 4)^4`.
#'   For EWMA: specifies the alpha parameter (0-1) for traditional exponential smoothing.
#'   Cannot be used simultaneously with `window` for EWMA method.
#'   For kernel: multiplier of optimal bandwidth (1.0 = optimal, <1 = less smooth, >1 = more smooth).
#'   For kalman: a finite, positive ratio of measurement to process noise
#'   (higher = more smoothing). An explicit noise variance in `params` determines
#'   the other variance from this ratio. If both variances are supplied, they
#'   take precedence over `smoothing`. Without a ratio, unspecified measurement
#'   and process variances default to 0.1 and 0.01 times the series variance.
#'   For others: typically 0-1 range.
#' @param band Unified band parameter for bandpass filters
#'   (bk, cf). Provide as `c(low, high)`: the shortest and longest cycle to
#'   remove, in periods of the series (months for monthly data). Both values
#'   must be positive. Defaults to cycles of 1.5 to 8 years: `c(6, 32)` for
#'   quarterly data, `c(18, 96)` for monthly, and `c(2, 8)` for annual.
#' @param align Unified alignment parameter for moving average
#'   methods (ma, wma, triangular, gaussian). Valid values: `"center"` (default, uses
#'   surrounding values), `"right"` (causal, uses past values only), `"left"` (anti-causal,
#'   uses future values only). Note: triangular only supports `"center"` and `"right"`.
#'   If NULL, uses `"center"` as default.
#' @param params Optional list of method-specific parameters for fine control:
#'   - **HP Filter**: `hp_onesided` (logical, default FALSE) - Use one-sided (real-time) filter instead of two-sided
#'   - **STL**: `stl_s_window` or `s.window` (numeric/"periodic", default "periodic") - Seasonal window,
#'     `stl_t_window` or `t.window` (numeric/NULL, default NULL) - Trend window,
#'     `stl_robust` or `robust` (logical, default FALSE) - Use robust fitting.
#'     Note: Both dot notation (`s.window`) and underscore notation (`stl_s_window`) are accepted.
#'   - **Spline**: `spline_cv` (logical/NULL) - Cross-validation method: NULL (none), TRUE (leave-one-out), FALSE (GCV)
#'   - **Polynomial**: `poly_degree` (integer, default 1), `poly_raw` (logical, default FALSE for orthogonal polynomials)
#'   - **UCM**: `ucm_type` (character) - Model type: "level", "trend", or "BSM".
#'     Defaults to "BSM" for frequencies 2 to 12 and "level" otherwise.
#'     Explicit "BSM" requests require frequency at most 12.
#'   - **Others**: `bn_ar_order`, `hamilton_h`, `hamilton_p`,
#'     `kernel_type`, `kalman_measurement_noise`, `kalman_process_noise`,
#'     `median_endrule`, `gaussian_sigma`, `wma_weights`.
#'   - **Note**: Alignment parameters (`ma_align`, `wma_align`, `triangular_align`, `gaussian_align`)
#'     can still be passed via `params` but it's recommended to use the unified `align` parameter instead.
#' @param .quiet If `TRUE`, suppress informational messages.
#'
#' @return If single method, returns a `ts` object. If multiple methods, returns
#'   a named list of `ts` objects.
#'
#' @importFrom cli cli_abort cli_inform cli_warn
#' @importFrom stats is.ts frequency start time ts fitted lm poly loess smooth.spline stl filter var AIC residuals
#' @importFrom hpfilter hp1 hp2
#' @importFrom tsbox ts_ts ts_df
#' @importFrom lubridate year month quarter
#' @importFrom RcppRoll roll_mean
#'
#' @details
#' This function focuses on monthly (frequency = 12) and quarterly (frequency = 4)
#' economic data. It uses established econometric methods with appropriate defaults:
#'
#' - **HP Filter**: lambda = 1600 (quarterly), 129600 (monthly), 6.25 (annual), following Ravn and Uhlig (2002). Supports both two-sided and one-sided (real-time) variants
#' - **Baxter-King**: Bandpass filter for business cycles (1.5 to 8 years by default)
#' - **Christiano-Fitzgerald**: Asymmetric bandpass filter
#' - **Moving Average**: Centered, frequency-appropriate windows
#' - **STL**: Seasonal-trend decomposition
#' - **Loess**: Local polynomial regression
#' - **Spline**: Smoothing splines
#' - **Polynomial**: Linear/polynomial trends
#' - **Beveridge-Nelson**: Permanent/transitory decomposition
#' - **UCM**: Unobserved Components Model (basic structural model up to monthly data, local level otherwise)
#' - **Hamilton**: Regression-based alternative to HP filter
#' - **Advanced MA**: EWMA with various implementations
#' - **Kernel Smoother**: Non-parametric regression with various kernel functions
#' - **Kalman Smoother**: Adaptive filtering for noisy time series
#' - **Median Filter**: Robust filtering using running medians to remove outliers
#' - **Gaussian Filter**: Weighted average with Gaussian (normal) density weights
#'
#' **Parameter Usage Notes**:
#' - **HP Filter**: Use `hp_onesided=TRUE` for real-time analysis or when future data should not
#'   influence current estimates. One-sided filter is appropriate for nowcasting, policy analysis,
#'   and avoiding look-ahead bias. Default two-sided filter is optimal for historical analysis.
#' - **EWMA**: Use either `window` (converted to `alpha = 2 / (window + 1)`) OR `smoothing` (alpha parameter), not both
#' - **Kalman**: Use `smoothing` parameter or `params` list for fine control of noise parameters
#' - **Spline**: Use `spline_cv` to control cross-validation (NULL=none, TRUE=LOO-CV, FALSE=GCV)
#' - **Polynomial**: Use `poly_raw=FALSE` for orthogonal polynomials (more stable for degree > 2)
#'   or `poly_raw=TRUE` for raw polynomials. Warning issued for degree > 3 (overfitting risk).
#' - **UCM**: Choose model type - "level" (simplest), "trend" (time-varying slope), or
#'   "BSM" (with seasonal component, requires seasonal data). Variances are
#'   estimated by maximum likelihood, so `smoothing` does not apply. The trend
#'   is the smoothed level. On seasonal data, "level" and "trend" can absorb
#'   the seasonality into the level and return the series itself, which is
#'   why the default is "BSM" up to monthly data. "BSM" carries one state per
#'   season and becomes very slow on weekly or daily data.
#'
#' @examples
#' # Single method
#' hp_trend <- extract_trends(AirPassengers, methods = "hp")
#'
#' # Multiple methods with unified smoothing
#' smooth_trends <- extract_trends(
#'   AirPassengers,
#'   methods = c("hp", "loess", "ewma"),
#'   smoothing = 0.3
#' )
#'
#' # EWMA with window (alpha derived from window size)
#' ewma_window <- extract_trends(AirPassengers, methods = "ewma", window = 12)
#'
#' # EWMA with alpha (traditional formula)
#' ewma_alpha <- extract_trends(AirPassengers, methods = "ewma", smoothing = 0.2)
#'
#' # Moving averages with unified window
#' ma_trends <- extract_trends(
#'   AirPassengers,
#'   methods = c("ma", "wma", "triangular"),
#'   window = 8
#' )
#'
#' # Bandpass filters with unified band
#' bp_trends <- extract_trends(
#'   AirPassengers,
#'   methods = c("bk", "cf"),
#'   band = c(18, 96)
#' )
#'
#' # Moving average with right alignment (causal filter)
#' ma_causal <- extract_trends(
#'   AirPassengers,
#'   methods = "ma",
#'   window = 12,
#'   align = "right"
#' )
#'
#' # Signal processing methods with specific parameters
#' finance_trends <- extract_trends(
#'   AirPassengers,
#'   methods = c("kalman", "gaussian"),
#'   window = 9,  # For Gaussian filter
#'   params = list(kalman_measurement_noise = 0.1)  # Kalman-specific parameter
#' )
#'
#' # Spline with cross-validation options
#' spline_trends <- extract_trends(
#'   AirPassengers,
#'   methods = "spline",
#'   params = list(spline_cv = FALSE)  # Use GCV instead of default
#' )
#'
#' # Polynomial with orthogonal vs raw polynomials
#' poly_trends <- extract_trends(
#'   AirPassengers,
#'   methods = "poly",
#'   params = list(poly_degree = 2, poly_raw = FALSE)  # Orthogonal (default)
#' )
#'
#' # UCM with different model types
#' ucm_trends <- extract_trends(
#'   AirPassengers,
#'   methods = "ucm",
#'   params = list(ucm_type = "BSM")  # Basic Structural Model with seasonality
#' )
#'
#' # HP Filter: One-sided (real-time) vs Two-sided (historical)
#' hp_realtime <- extract_trends(
#'   AirPassengers,
#'   methods = "hp",
#'   params = list(hp_onesided = TRUE)  # For nowcasting and real-time analysis
#' )
#'
#' # STL with custom parameters via params (both notations work)
#' stl_custom1 <- extract_trends(
#'   AirPassengers,
#'   methods = "stl",
#'   params = list(s.window = 21, robust = TRUE)  # Dot notation
#' )
#'
#' stl_custom2 <- extract_trends(
#'   AirPassengers,
#'   methods = "stl",
#'   params = list(stl_s_window = 21, stl_robust = TRUE)  # Underscore notation
#' )
#'
#' # Advanced: fine-tune specific methods
#' custom_trends <- extract_trends(
#'   AirPassengers,
#'   methods = c("median", "kalman"),
#'   window = 7,
#'   params = list(median_endrule = "constant")
#' )
#'
#' @export
extract_trends <- function(
  ts_data,
  methods = "stl",
  window = NULL,
  smoothing = NULL,
  band = NULL,
  align = NULL,
  params = list(),
  .quiet = FALSE
) {
  # Input validation
  if (is.null(ts_data)) {
    cli::cli_abort("{.arg ts_data} cannot be NULL")
  }

  # Validate methods
  .validate_methods(methods)

  # Convert to ts object using tsbox if needed
  if (!stats::is.ts(ts_data)) {
    tryCatch(
      {
        ts_data <- tsbox::ts_ts(ts_data)
      },
      error = function(e) {
        cli::cli_abort(c(
          "Failed to convert input to time series object.",
          "i" = "Input must be convertible to ts via tsbox package.",
          "x" = "Error: {e$message}"
        ))
      }
    )
  }

  # Validate unified parameters
  .validate_unified_params(window, smoothing, band, align, params)

  # Validate frequency
  freq <- stats::frequency(ts_data)
  if (freq < 1 || freq > 365) {
    cli::cli_abort(
      "Frequency must be between 1 (annual) and 365 (daily), got {freq}"
    )
  }

  # Reject interior gaps, and set leading/trailing missing values aside so the
  # methods see a complete series. The data.frame path never reaches here with
  # an NA: it drops those rows and validates the period grid beforehand.
  trimmed <- .trim_to_observed(ts_data, .observed_span(ts_data))
  na_template <- if (is.null(trimmed)) NULL else ts_data
  if (!is.null(trimmed)) {
    ts_data <- trimmed
  }

  # Warn for frequency-sensitive methods with non-standard frequencies
  if (!freq %in% c(1, 4, 12) && any(methods %in% .FREQ_SENSITIVE_METHODS)) {
    freq_sensitive <- intersect(methods, .FREQ_SENSITIVE_METHODS)
    if (!.quiet) {
      cli::cli_warn(
        "Methods {.val {freq_sensitive}} are optimized for standard economic frequencies.
         Using frequency = {freq} with default parameters may produce suboptimal results.
         Consider specifying parameters explicitly via the {.arg params} argument."
      )
    }
  }

  # The default lambda was calibrated on annual to monthly data. Warn
  # regardless of .quiet: on weekly or daily data it rarely gives the trend the
  # user expects, and the result would otherwise look plausible.
  hp_lambda_set <- !is.null(smoothing) || !is.null(params$hp_lambda)
  if ("hp" %in% methods && freq > 12 && !hp_lambda_set) {
    default_lambda <- .default_hp_lambda(freq)
    cli::cli_warn(c(
      "The HP filter is rarely used on data with frequency {freq}.",
      "i" = "Using the default {.code lambda = {format(default_lambda, scientific = FALSE)}}.",
      "i" = "Set {.arg smoothing} to choose lambda explicitly."
    ))
  }

  # Check minimum observations
  min_obs <- 3 * freq
  if (length(ts_data) < min_obs && !.quiet) {
    cli::cli_warn(
      "Series has {length(ts_data)} observations.
       Minimum {min_obs} recommended for reliable trend extraction."
    )
  }

  # Handle vector window: expand ma/median methods into one call per window value
  if (!is.null(window) && length(window) > 1) {
    window_methods <- intersect(methods, .WINDOW_VECTOR_METHODS)
    other_methods <- setdiff(methods, .WINDOW_VECTOR_METHODS)

    first_window_methods <- intersect(other_methods, .WINDOW_METHODS)
    if (length(window_methods) == 0) {
      first_window_methods <- other_methods
    }
    if (length(first_window_methods) > 0) {
      cli::cli_warn(c(
        "Multiple {.arg window} values are only supported for {.val ma}, {.val median}, and {.val henderson} methods.",
        "i" = "Using first value ({window[1]}) for method(s) {.val {first_window_methods}}."
      ))
    }
    if (length(window_methods) == 0) {
      window <- window[1]
    } else {
      results <- list()

      for (method in other_methods) {
        results[[method]] <- .extract_trends_impl(
          ts_data = ts_data,
          methods = method,
          freq = freq,
          na_template = na_template,
          window = window[1],
          smoothing = smoothing,
          band = band,
          align = align,
          params = params,
          .quiet = .quiet
        )
      }

      for (w in window) {
        for (method in window_methods) {
          results[[paste0(method, "_", w)]] <- .extract_trends_impl(
            ts_data = ts_data,
            methods = method,
            freq = freq,
            na_template = na_template,
            window = w,
            smoothing = smoothing,
            band = band,
            align = align,
            params = params,
            .quiet = .quiet
          )
        }
      }

      if (length(results) == 1) {
        return(.restore_time_base(results[[1]], na_template))
      }
      return(.restore_time_base(results, na_template))
    }
  }

  return(.extract_trends_impl(
    ts_data = ts_data,
    methods = methods,
    freq = freq,
    na_template = na_template,
    window = window,
    smoothing = smoothing,
    band = band,
    align = align,
    params = params,
    .quiet = .quiet
  ))
}

.extract_trends_impl <- function(
  ts_data,
  methods,
  freq,
  na_template,
  window,
  smoothing,
  band,
  align,
  params,
  .quiet
) {
  # Process unified parameters to get method-specific parameters
  unified_params <- .process_unified_params(
    methods = methods,
    window = window,
    smoothing = smoothing,
    band = band,
    align = align,
    params = params,
    frequency = freq,
    .quiet = .quiet
  )

  # Extract parameters from unified system with defaults
  .get_param <- function(name, default) unified_params[[name]] %||% default

  # Method-specific parameters
  hp_lambda <- .get_param("hp_lambda", .default_hp_lambda(freq))
  hp_onesided <- .get_param("hp_onesided", FALSE)
  ma_window <- .get_param("ma_window", freq)
  ma_align <- .get_param("ma_align", "center")
  stl_s_window <- .get_param("stl_s_window", "periodic")
  stl_t_window <- .get_param("stl_t_window", NULL)
  stl_robust <- .get_param("stl_robust", FALSE)
  loess_span <- .get_param("loess_span", 0.75)
  spline_spar <- .get_param("spline_spar", NULL)
  spline_cv <- .get_param("spline_cv", NULL)
  poly_degree <- .get_param("poly_degree", 1)
  poly_raw <- .get_param("poly_raw", FALSE)
  bn_ar_order <- .get_param("bn_ar_order", NULL)
  # A level model on seasonal data absorbs the seasonality into the level.
  # BSM carries one state per season, so it is impractical above monthly.
  ucm_default <- if (freq > 1 && freq <= 12) "BSM" else "level"
  ucm_type <- .get_param("ucm_type", ucm_default)
  # Frequency-aware defaults (Hamilton 2018): h = 8, p = 4 for quarterly,
  # h = 24, p = 12 for monthly, etc.
  hamilton_defaults <- .get_hamilton_params(freq)
  hamilton_h <- .get_param("hamilton_h", hamilton_defaults$h)
  hamilton_p <- .get_param("hamilton_p", hamilton_defaults$p)
  ewma_window <- .get_param("ewma_window", NULL)
  ewma_alpha <- .get_param("ewma_alpha", NULL)
  wma_window <- .get_param("wma_window", freq)
  wma_weights <- .get_param("wma_weights", NULL)
  wma_align <- .get_param("wma_align", "center")
  triangular_window <- .get_param("triangular_window", freq)
  triangular_align <- .get_param("triangular_align", "center")
  default_band <- .default_band(freq)
  bk_low <- .get_param("bk_low", default_band[1])
  bk_high <- .get_param("bk_high", default_band[2])
  cf_low <- .get_param("cf_low", default_band[1])
  cf_high <- .get_param("cf_high", default_band[2])
  kernel_bandwidth <- .get_param("kernel_bandwidth", NULL)
  kernel_type <- .get_param("kernel_type", "normal")
  kalman_smoothing <- .get_param("kalman_smoothing", NULL)
  kalman_measurement_noise <- .get_param("kalman_measurement_noise", NULL)
  kalman_process_noise <- .get_param("kalman_process_noise", NULL)
  median_window <- .get_param("median_window", 5)
  median_endrule <- .get_param("median_endrule", "median")
  gaussian_window <- .get_param("gaussian_window", 7)
  gaussian_sigma <- .get_param("gaussian_sigma", NULL)
  gaussian_align <- .get_param("gaussian_align", "center")
  henderson_window <- .get_param("henderson_window", NULL)

  # Extract trends
  trends <- list()

  for (method in methods) {
    trend <- switch(
      method,
      "hp" = .extract_hp_trend(ts_data, hp_lambda, hp_onesided, .quiet),
      "bk" = .extract_bk_trend(ts_data, bk_low, bk_high, .quiet),
      "cf" = .extract_cf_trend(ts_data, cf_low, cf_high, .quiet),
      "ma" = .extract_ma_trend(ts_data, ma_window, ma_align, .quiet),
      "stl" = .extract_stl_trend(
        ts_data,
        stl_s_window,
        stl_t_window,
        stl_robust,
        .quiet
      ),
      "loess" = .extract_loess_trend(ts_data, loess_span, .quiet),
      "spline" = .extract_spline_trend(ts_data, spline_spar, spline_cv, .quiet),
      "poly" = .extract_poly_trend(ts_data, poly_degree, poly_raw, .quiet),
      "bn" = .extract_bn_trend(ts_data, bn_ar_order, .quiet),
      "ucm" = .extract_ucm_trend(ts_data, ucm_type, .quiet),
      "hamilton" = .extract_hamilton_trend(
        ts_data,
        hamilton_h,
        hamilton_p,
        .quiet
      ),
      "spencer" = .extract_spencer_trend(ts_data, .quiet),
      "ewma" = .extract_ewma_trend(ts_data, ewma_window, ewma_alpha, .quiet),
      "wma" = .extract_wma_trend(
        ts_data,
        wma_window,
        wma_weights,
        wma_align,
        .quiet
      ),
      "triangular" = .extract_triangular_trend(
        ts_data,
        triangular_window,
        triangular_align,
        .quiet
      ),
      "kernel" = .extract_kernel_trend(
        ts_data,
        kernel_bandwidth,
        kernel_type,
        .quiet
      ),
      "kalman" = .extract_kalman_trend(
        ts_data,
        kalman_measurement_noise,
        kalman_process_noise,
        .quiet,
        smoothing = kalman_smoothing
      ),
      "median" = .extract_median_trend(
        ts_data,
        median_window,
        median_endrule,
        .quiet
      ),
      "gaussian" = .extract_gaussian_trend(
        ts_data,
        gaussian_window,
        gaussian_sigma,
        gaussian_align,
        .quiet
      ),
      "henderson" = .extract_henderson_trend(ts_data, henderson_window, .quiet)
    )

    trends[[method]] <- trend
  }

  # Return single ts if one method, list if multiple
  if (length(methods) == 1) {
    return(.restore_time_base(trends[[1]], na_template))
  } else {
    return(.restore_time_base(trends, na_template))
  }
}