R/main.R

In spooky: Time Feature Extrapolation Using Spectral Analysis and Jack-Knife Resampling

#' Spectral forecasting with jackknife resampling
#'
#' @param df Numeric or categorical time-feature data frame.
#' @param seq_len Forecast horizon or range of horizons to search.
#' @param lno Jackknife leave-out value or search range.
#' @param n_samp Number of candidate configurations.
#' @param n_windows Number of validation windows.
#' @param ci Confidence level for prediction summaries.
#' @param smoother Apply a moving-average smoother to numeric data.
#' @param dates Optional Date vector matching `df`.
#' @param error_scale Scale used by numeric error metrics.
#' @param error_benchmark Benchmark used by relative error metrics.
#' @param seed Random seed.
#' @return An object of class `spooky_fit`.
#' @export
spooky <- function(df, seq_len = NULL, lno = NULL, n_samp = 30,
                   n_windows = 3, ci = 0.8, smoother = FALSE,
                   dates = NULL, error_scale = "naive",
                   error_benchmark = "naive", seed = 42) {
  if (!is.data.frame(df) || !nrow(df) || !ncol(df)) {
    stop("df must be a non-empty data frame")
  }
  numeric_columns <- vapply(df, is.numeric, logical(1))
  categorical_columns <- vapply(
    df,
    function(x) is.factor(x) || is.character(x),
    logical(1)
  )
  if (!(all(numeric_columns) || all(categorical_columns))) {
    stop("use either all numeric or all categorical columns")
  }
  if (!is.null(dates) &&
      (!inherits(dates, "Date") || length(dates) != nrow(df))) {
    stop("dates must match df and be Date")
  }
  if (ci <= 0 || ci >= 1 || n_samp < 1 || n_windows < 1) {
    stop("invalid search or confidence parameters")
  }
  if (is.null(seq_len)) {
    seq_len <- c(1, max(1, round(sqrt(nrow(df)))))
  }
  if (is.null(lno)) {
    lno <- c(1, max(1, round(sqrt(nrow(df)))))
  }

  set.seed(seed)
  binary <- all(categorical_columns)
  if (binary) {
    levels <- sort(unique(unlist(lapply(df, as.character))))
    df <- as.data.frame(lapply(df, function(x) {
      as.numeric(factor(as.character(x), levels = levels)) - 1
    }))
  } else {
    df <- as.data.frame(lapply(df, function(x) {
      x <- as.numeric(x)
      if (anyNA(x)) {
        x <- stats::approx(
          which(!is.na(x)), x[!is.na(x)], seq_along(x), rule = 2
        )$y
      }
      x
    }))
  }
  if (smoother && !binary) {
    df <- as.data.frame(lapply(df, function(x) {
      smoothed <- stats::filter(x, rep(1 / 3, 3), sides = 2)
      smoothed[is.na(smoothed)] <- x[is.na(smoothed)]
      smoothed
    }))
  }
  if (length(seq_len) == 1 && length(lno) == 1) n_samp <- 1
  sequence_values <- if (n_samp == 1) {
    seq_len[1]
  } else {
    round(stats::runif(n_samp, min(seq_len), max(seq_len)))
  }
  leave_out_values <- if (n_samp == 1) {
    lno[1]
  } else {
    round(stats::runif(n_samp, min(lno), max(lno)))
  }

  candidates <- lapply(seq_len(n_samp), function(i) {
    .candidate(
      df, sequence_values[i], leave_out_values[i], n_windows, ci,
      binary, error_scale, error_benchmark
    )
  })
  scores <- do.call(rbind, lapply(candidates, `[[`, "avg_errors"))
  best <- if (n_samp == 1) 1 else which.min(rowMeans(scale(scores)))
  chosen <- candidates[[best]]
  history <- data.frame(
    seq_len = sequence_values,
    lno = leave_out_values,
    scores,
    check.names = FALSE
  )
  errors <- do.call(rbind, chosen$errors)
  rownames(errors) <- colnames(df)
  result <- list(
    exploration = candidates,
    history = history,
    best_model = list(
      testing_errors = errors,
      preds = lapply(chosen$model, `[[`, "quantile_preds"),
      plots = lapply(chosen$model, `[[`, "plot")
    ),
    time_log = proc.time()[["elapsed"]]
  )
  class(result) <- c("spooky_fit", "list")
  result
}

.forecast <- function(x, horizon, leave_out) {
  if (length(x) < 3) {
    return(matrix(rep(utils::tail(x, 1), horizon), nrow = 1))
  }
  differenced <- diff(x)
  integrated <- TRUE
  if (!is.finite(stats::sd(differenced)) || stats::sd(differenced) == 0) {
    differenced <- x
    integrated <- FALSE
  }
  spectrum <- stats::fft(differenced)
  keep <- seq_len(max(1, length(spectrum) -
    min(leave_out - 1, length(spectrum) - 1)))
  padded <- c(spectrum[keep], complex(length.out = horizon))
  forecast <- Re(stats::fft(padded, inverse = TRUE) / length(padded))
  matrix(
    utils::tail(forecast, horizon) +
      if (integrated) utils::tail(x, 1) else 0,
    nrow = 1
  )
}

.metrics <- function(actual, prediction, history, binary, scale, benchmark) {
  if (binary) {
    actual_binary <- round(actual)
    prediction_binary <- round(prediction)
    return(c(
      dice = 2 * sum(actual_binary * prediction_binary) /
        max(sum(actual_binary) + sum(prediction_binary), 1),
      mae = mean(abs(actual - prediction))
    ))
  }
  error <- actual - prediction
  metric_scale <- if (scale == "deviation") {
    stats::sd(history)
  } else {
    mean(abs(diff(history)))
  }
  baseline <- if (benchmark == "average") mean(history) else utils::tail(history, 1)
  c(
    me = mean(error),
    mae = mean(abs(error)),
    mse = mean(error^2),
    rmsse = sqrt(mean(error^2)) / max(metric_scale, 1e-12),
    mape = mean(abs(error / actual), na.rm = TRUE),
    rmae = mean(abs(error)) /
      max(mean(abs(actual - baseline)), 1e-12)
  )
}

.summary <- function(predictions, ci) {
  probabilities <- sort(unique(c(
    (1 - ci) / 2, 0.25, 0.5, 0.75, (1 + ci) / 2
  )))
  summary_values <- t(apply(predictions, 2, function(x) {
    c(
      min(x),
      stats::quantile(x, probabilities),
      max(x),
      mean(x),
      stats::sd(x)
    )
  }))
  colnames(summary_values) <- c(
    "min", paste0(probabilities * 100, "%"), "max", "mean", "sd"
  )
  as.data.frame(round(summary_values, 5))
}

.candidate <- function(df, horizon, leave_out, n_windows, ci, binary,
                       scale, benchmark) {
  if (floor(nrow(df) / (n_windows + 1)) <= horizon) {
    stop("not enough data for validation windows")
  }
  errors <- vector("list", ncol(df))
  models <- vector("list", ncol(df))
  cuts <- floor(seq(1, nrow(df), length.out = n_windows + 1))
  for (feature in seq_len(ncol(df))) {
    window_errors <- list()
    window_predictions <- list()
    for (window in seq_len(n_windows)) {
      training <- df[seq_len(cuts[window]), feature]
      truth <- df[(cuts[window] + 1):cuts[window + 1], feature]
      predictions <- .forecast(training, length(truth), leave_out)
      window_errors[[window]] <- .metrics(
        truth, colMeans(predictions), training, binary, scale, benchmark
      )
      window_predictions[[window]] <- predictions
    }
    common_horizon <- min(vapply(window_predictions, ncol, integer(1)))
    predictions <- do.call(rbind, lapply(window_predictions, function(x) {
      x[, seq_len(common_horizon), drop = FALSE]
    }))
    errors[[feature]] <- colMeans(do.call(rbind, window_errors))
    models[[feature]] <- list(
      quantile_preds = .summary(predictions, ci),
      plot = function(...) invisible(NULL)
    )
  }
  list(
    errors = errors,
    avg_errors = colMeans(do.call(rbind, errors)),
    model = models
  )
}

#' @export
print.spooky_fit <- function(x, ...) {
  cat(
    "Spooky 2.0 fit\nFeatures:", length(x$best_model$preds),
    "\nCandidates:", nrow(x$history), "\n"
  )
}

#' @export
plot.spooky_fit <- function(x, feature = 1, ...) {
  invisible(NULL)
}