Time Series Forecast — See Where Your Numbers Are Heading

suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

raw_date <- df$date
  non_blank <- !(is.na(raw_date) | !nzchar(trimws(as.character(raw_date))))
  d <- parse_dates_robust(raw_date)
  n_nonblank <- sum(non_blank)
  if (n_nonblank == 0) {
    stop(sprintf("The &#x27;%s' column is empty — no dates to forecast from.", date_name))
  }
  n_unparsed <- sum(non_blank & is.na(d))
  if (n_unparsed / n_nonblank > 0.05) {
    stop(sprintf(
      "%d of %d values in &#x27;%s' (%.0f%%) could not be read as dates. Please use a recognizable date format (e.g. 2024-01-31 or 01/31/2024).",
      n_unparsed, n_nonblank, date_name, 100 * n_unparsed / n_nonblank))
  }

v <- df$value
  if (!is.numeric(v)) {
    conv <- suppressWarnings(as.numeric(as.character(v)))
    n_orig <- sum(!is.na(v) & nzchar(as.character(v)))
    if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
      v <- conv
    } else {
      stop(sprintf("The &#x27;%s' column is not numeric — a forecast needs a numeric metric.",
                   value_name))
    }
  }
  v <- as.numeric(v)

  keep <- !is.na(d) & !is.na(v)
  n_dropped_rows <- initial_rows - sum(keep)
  work <- data.frame(date = d[keep], value = v[keep])

gaps <- as.numeric(diff(agg$date))
  med_gap <- median(gaps)
  if (med_gap <= 1.5) {
    freq_label <- "daily";   step_days <- 1;  freq <- 7;  period_word <- "day";   cycle_word <- "week"
  } else if (med_gap >= 5.5 && med_gap <= 8.5) {
    freq_label <- "weekly";  step_days <- 7;  freq <- 52; period_word <- "week";  cycle_word <- "year"
  } else if (med_gap >= 26 && med_gap <= 35) {
    freq_label <- "monthly"; step_days <- NA; freq <- 12; period_word <- "month"; cycle_word <- "year"
  } else {
    freq_label <- "irregular"; step_days <- max(1, round(med_gap)); freq <- 1
    period_word <- sprintf("%d-day period", step_days); cycle_word <- "cycle"
  }

if (freq_label == "monthly") {
    snapped <- as.Date(format(agg$date, "%Y-%m-01"))
    agg <- aggregate(value ~ date, data = data.frame(date = snapped, value = agg$value),
                     FUN = sum)
    agg <- agg[order(agg$date), , drop = FALSE]
    grid <- seq(min(agg$date), max(agg$date), by = "month")
  } else {
    grid <- seq(min(agg$date), max(agg$date), by = step_days)
  }
  interp <- approx(x = as.numeric(agg$date), y = agg$value,
                   xout = as.numeric(grid), rule = 2)
  values <- interp$y
  n_interp <- sum(!(as.numeric(grid) %in% as.numeric(agg$date)))
  if (freq_label != "irregular" && n_interp / length(grid) > 0.5) {
    stop(sprintf(
      "More than half of the %s periods between the first and last &#x27;%s' are missing — too sparse to forecast reliably.",
      freq_label, date_name))
  }

  n_points <- length(values)
  final_rows <- n_points
  rows_removed <- max(0, initial_rows - sum(keep))

horizon <- max(2, min(floor(0.25 * n_points), 2 * max(freq, 3), 30))
  fitres <- fit_chain(values, freq, horizon)
  method <- fitres$method
  method_label <- fitres$method_label
  seasonal_fit <- isTRUE(fitres$seasonal)

  if (freq_label == "monthly") {
    future_dates <- seq(grid[length(grid)], by = "month", length.out = horizon + 1)[-1]
  } else {
    future_dates <- grid[length(grid)] + step_days * seq_len(horizon)
  }

mape <- NA_real_; mae <- NA_real_; holdout_n <- 0L
  backtest_note <- NULL
  holdout_target <- min(max(freq, 6), floor(n_points / 3))
  if (holdout_target >= 4 && (n_points - holdout_target) >= 10) {
    bt <- tryCatch({
      tr <- values[seq_len(n_points - holdout_target)]
      fit_chain(tr, freq, holdout_target)
    }, error = function(e) NULL)
    if (!is.null(bt)) {
      act <- values[(n_points - holdout_target + 1):n_points]
      pred <- bt$mean
      mae <- mean(abs(act - pred))
      nz <- abs(act) > 1e-9
      if (any(nz)) mape <- mean(abs((act[nz] - pred[nz]) / act[nz])) * 100
      holdout_n <- holdout_target
    } else {
      backtest_note <- "Backtest skipped — the model could not be refitted on the shortened training window."
    }
  } else {
    backtest_note <- "Backtest skipped — the series is too short to hold out a meaningful test window."
  }

last_actual <- values[n_points]
  fc_last <- fitres$mean[horizon]
  pct_change <- if (abs(last_actual) > 1e-9) {
    (fc_last - last_actual) / abs(last_actual) * 100
  } else NA_real_
  direction_word <- if (is.na(pct_change)) "move" else if (pct_change > 2) "rise" else if (pct_change < -2) "fall" else "hold roughly steady"

idx <- seq_len(n_points)
  trend_slope <- as.numeric(coef(lm(values ~ idx))[2])
  level_now <- last_actual
  seasonal_amplitude <- 0
  seasonal_coefs <- numeric(0)
  if (seasonal_fit && !is.null(fitres$fit)) {
    co <- coef(fitres$fit)
    if ("a" %in% names(co)) level_now <- as.numeric(co["a"])
    if ("b" %in% names(co)) trend_slope <- as.numeric(co["b"])
    seasonal_coefs <- as.numeric(co[grep("^s[0-9]+$", names(co))])
    if (length(seasonal_coefs) > 0) {
      seasonal_amplitude <- max(seasonal_coefs) - min(seasonal_coefs)
    }
  } else if (method == "hw_trend" && !is.null(fitres$fit)) {
    co <- coef(fitres$fit)
    if ("a" %in% names(co)) level_now <- as.numeric(co["a"])
    if ("b" %in% names(co)) trend_slope <- as.numeric(co["b"])
  }