R/heuristicSTR.R

In stR: Seasonal Trend Decomposition Using Regression

#' @rdname heuristicSTR
#' @name heuristicSTR
#' @title Automatic STR decomposition with heuristic search of the parameters
#' @description Automatically selects parameters (lambda coefficients) for an STR decomposition of time series data.
#' Heuristic approach can give a better estimate compare to a standard optmisaton methods used in \code{\link{STR}}.
#'
#' If a parallel backend is registered for use before \code{STR} call,
#' \code{heuristicSTR} will use it for n-fold cross validation computations.
#'
#' @seealso \code{\link{STR}} \code{\link{STRmodel}} \code{\link{AutoSTR}}
#' @inheritParams data
#' @inheritParams predictors
#' @inheritParams confidence
#' @inheritParams lambdas
#' @inheritParams pattern
#' @inheritParams nFold
#' @inheritParams reltol
#' @inheritParams gapCV
#' @inheritParams solver
#' @inheritParams trace
#' @param ratioGap Ratio to define hyperparameter bounds for one-dimensional search.
#' @param relCV Minimum improvement required after all predictors tried. It is used to exit heuristic serach of lambda parameters.
#' @templateVar class STR
#' @templateVar topLevel1 \item \strong{cvMSE} -- optional cross validated (leave one out) Mean Squared Error.
#' @templateVar topLevel2 \item \strong{optim.CV.MSE} or \strong{optim.CV.MAE} -- best cross validated Mean Squared Error or Mean Absolute Error (n-fold) achieved during minimisation procedure.
#' @templateVar topLevel3 \item \strong{nFold} -- the input \code{nFold} parameter.
#' @templateVar topLevel4 \item \strong{gapCV} -- the input \code{gapCV} parameter.
#' @templateVar topLevel5 \item \strong{method} -- contains strings \code{"STR"} or \code{"RSTR"} depending on used method.
#' @template returnValue
#' @examples
#' \donttest{
#' TrendSeasonalStructure <- list(
#'   segments = list(c(0, 1)),
#'   sKnots = list(c(1, 0))
#' )
#' WDSeasonalStructure <- list(
#'   segments = list(c(0, 48), c(100, 148)),
#'   sKnots = c(as.list(c(1:47, 101:147)), list(c(0, 48, 100, 148)))
#' )
#'
#' TrendSeasons <- rep(1, nrow(electricity))
#' WDSeasons <- as.vector(electricity[, "WorkingDaySeasonality"])
#'
#' Data <- as.vector(electricity[, "Consumption"])
#' Times <- as.vector(electricity[, "Time"])
#' TempM <- as.vector(electricity[, "Temperature"])
#' TempM2 <- TempM^2
#'
#' TrendTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 116)
#' SeasonTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 24)
#'
#' TrendData <- rep(1, length(Times))
#' SeasonData <- rep(1, length(Times))
#'
#' Trend <- list(
#'   name = "Trend",
#'   data = TrendData,
#'   times = Times,
#'   seasons = TrendSeasons,
#'   timeKnots = TrendTimeKnots,
#'   seasonalStructure = TrendSeasonalStructure,
#'   lambdas = c(1500, 0, 0)
#' )
#' WDSeason <- list(
#'   name = "Dayly seas",
#'   data = SeasonData,
#'   times = Times,
#'   seasons = WDSeasons,
#'   timeKnots = SeasonTimeKnots,
#'   seasonalStructure = WDSeasonalStructure,
#'   lambdas = c(0.003, 0, 240)
#' )
#' StaticTempM <- list(
#'   name = "Temp Mel",
#'   data = TempM,
#'   times = Times,
#'   seasons = NULL,
#'   timeKnots = NULL,
#'   seasonalStructure = NULL,
#'   lambdas = c(0, 0, 0)
#' )
#' StaticTempM2 <- list(
#'   name = "Temp Mel^2",
#'   data = TempM2,
#'   times = Times,
#'   seasons = NULL,
#'   timeKnots = NULL,
#'   seasonalStructure = NULL,
#'   lambdas = c(0, 0, 0)
#' )
#' Predictors <- list(Trend, WDSeason, StaticTempM, StaticTempM2)
#'
#' elec.fit <- heuristicSTR(
#'   data = Data,
#'   predictors = Predictors,
#'   gapCV = 48 * 7
#' )
#'
#' plot(elec.fit,
#'   xTime = as.Date("2000-01-11") + ((Times - 1) / 48 - 10),
#'   forecastPanels = NULL
#' )
#'
#' ########################################
#'
#' TrendSeasonalStructure <- list(
#'   segments = list(c(0, 1)),
#'   sKnots = list(c(1, 0))
#' )
#' DailySeasonalStructure <- list(
#'   segments = list(c(0, 48)),
#'   sKnots = c(as.list(1:47), list(c(48, 0)))
#' )
#' WeeklySeasonalStructure <- list(
#'   segments = list(c(0, 336)),
#'   sKnots = c(as.list(seq(4, 332, 4)), list(c(336, 0)))
#' )
#' WDSeasonalStructure <- list(
#'   segments = list(c(0, 48), c(100, 148)),
#'   sKnots = c(as.list(c(1:47, 101:147)), list(c(0, 48, 100, 148)))
#' )
#'
#' TrendSeasons <- rep(1, nrow(electricity))
#' DailySeasons <- as.vector(electricity[, "DailySeasonality"])
#' WeeklySeasons <- as.vector(electricity[, "WeeklySeasonality"])
#' WDSeasons <- as.vector(electricity[, "WorkingDaySeasonality"])
#'
#' Data <- as.vector(electricity[, "Consumption"])
#' Times <- as.vector(electricity[, "Time"])
#' TempM <- as.vector(electricity[, "Temperature"])
#' TempM2 <- TempM^2
#'
#' TrendTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 116)
#' SeasonTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 24)
#' SeasonTimeKnots2 <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 12)
#'
#' TrendData <- rep(1, length(Times))
#' SeasonData <- rep(1, length(Times))
#'
#' Trend <- list(
#'   name = "Trend",
#'   data = TrendData,
#'   times = Times,
#'   seasons = TrendSeasons,
#'   timeKnots = TrendTimeKnots,
#'   seasonalStructure = TrendSeasonalStructure,
#'   lambdas = c(1500, 0, 0)
#' )
#' WSeason <- list(
#'   name = "Weekly seas",
#'   data = SeasonData,
#'   times = Times,
#'   seasons = WeeklySeasons,
#'   timeKnots = SeasonTimeKnots2,
#'   seasonalStructure = WeeklySeasonalStructure,
#'   lambdas = c(0.8, 0.6, 100)
#' )
#' WDSeason <- list(
#'   name = "Dayly seas",
#'   data = SeasonData,
#'   times = Times,
#'   seasons = WDSeasons,
#'   timeKnots = SeasonTimeKnots,
#'   seasonalStructure = WDSeasonalStructure,
#'   lambdas = c(0.003, 0, 240)
#' )
#' TrendTempM <- list(
#'   name = "Trend temp Mel",
#'   data = TempM,
#'   times = Times,
#'   seasons = TrendSeasons,
#'   timeKnots = TrendTimeKnots,
#'   seasonalStructure = TrendSeasonalStructure,
#'   lambdas = c(1e7, 0, 0)
#' )
#' TrendTempM2 <- list(
#'   name = "Trend temp Mel^2",
#'   data = TempM2,
#'   times = Times,
#'   seasons = TrendSeasons,
#'   timeKnots = TrendTimeKnots,
#'   seasonalStructure = TrendSeasonalStructure,
#'   lambdas = c(0.01, 0, 0)
#' ) # Starting parameter is too far from the optimal value
#' Predictors <- list(Trend, WSeason, WDSeason, TrendTempM, TrendTempM2)
#'
#' elec.fit <- heuristicSTR(
#'   data = Data,
#'   predictors = Predictors,
#'   gapCV = 48 * 7
#' )
#'
#' plot(elec.fit,
#'   xTime = as.Date("2000-01-11") + ((Times - 1) / 48 - 10),
#'   forecastPanels = NULL
#' )
#'
#' plotBeta(elec.fit, predictorN = 4)
#' plotBeta(elec.fit, predictorN = 5)
#' }
#' @author Alexander Dokumentov
#' @references Dokumentov, A., and Hyndman, R.J. (2022)
#' STR: Seasonal-Trend decomposition using Regression,
#' \emph{INFORMS Journal on Data Science}, 1(1), 50-62.
#' \url{https://robjhyndman.com/publications/str/}
#' @export

heuristicSTR <- function(data, predictors,
                         confidence = NULL, lambdas = NULL,
                         pattern = extractPattern(predictors), nFold = 5, reltol = 0.005, gapCV = 1,
                         solver = c("Matrix", "cholesky"),
                         trace = FALSE,
                         ratioGap = 1e12, # Ratio to define bounds for one-dimensional search.
                         relCV = 0.01 # Minimum improvement required after all predictors tried. It is used to exit.
) {
  lambdas <- upLambdas(
    data = data, predictors = predictors,
    confidence = confidence, lambdas = lambdas,
    pattern = pattern, nFold = nFold,
    reltol = reltol, gapCV = gapCV,
    solver = solver, trace = trace, ratioGap = ratioGap
  )
  result <- STR(
    data = data, predictors = predictors,
    confidence = confidence, lambdas = lambdas,
    pattern = pattern, nFold = nFold,
    reltol = reltol, gapCV = gapCV,
    solver = solver, trace = trace
  )
  return(result)
}

# Tries to find optimal lambda parameters optimising them in groups (one group per predictor), over
# the residuals and allowing only to increase L1 norm of the group of the parameters
upLambdas <- function(data, predictors,
                      confidence = NULL, lambdas = NULL,
                      pattern = extractPattern(predictors), nFold = 5, reltol = 0.005, gapCV = 1,
                      solver = c("Matrix", "cholesky"),
                      trace = FALSE,
                      ratioGap = 1e6, # Ratio to define bounds for one-dimensional search
                      nPasses = 2,
                      maxLambda = 1e6) {
  if (is.null(lambdas)) {
    lambdas <- lapply(predictors, FUN = function(p) list(lambdas = p$lambdas))
  } else {
    lambdas <- lapply(lambdas, FUN = function(p) list(lambdas = p$lambdas))
  }

  lData <- length(data)
  if (nFold * gapCV > lData) {
    stop(paste0(
      "nFold*gapCV should be less or equal to the data length.\nnFold = ",
      nFold, "\ngapCV = ", gapCV, "\nlength of the data = ", lData
    ))
  }
  subInds <- lapply(1:nFold, FUN = function(i) sort(unlist(lapply(1:gapCV, FUN = function(j) seq(from = (i - 1) * gapCV + j, to = lData, by = nFold * gapCV)))))
  complInds <- lapply(subInds, FUN = function(s) setdiff(1:lData, s))

  strDesign <- STRDesign(predictors)
  C <- strDesign$cm$matrix
  fcastC <- lapply(subInds, FUN = function(si) C[si, ])
  trainC <- lapply(complInds, FUN = function(ci) C[ci, ])
  fcastData <- lapply(subInds, FUN = function(si) data[si])
  trainData <- lapply(complInds, FUN = function(ci) data[ci])
  rm <- strDesign$rm
  regMatrix <- rm$matrix
  regSeats <- rm$seats

  newLambdas <- lambdas
  if (trace) {
    dummy <- lapply(newLambdas, FUN = function(p) {
      cat(format(p$lambdas, digits = 4, scientific = T, width = 9))
      cat("  ")
    })
  }
  oldLambdas <- oldFit <- fit <- NULL
  for (j in seq_len(nPasses)) {
    for (i in seq_along(predictors)) {
      if (sum(abs(newLambdas[[i]]$lambdas)) == 0) next
      if (trace) {
        cat("\nPredictor: ")
        cat(i)
        cat("\n")
      }
      if (is.null(fit)) {
        fit <- try(STRmodel(data, strDesign = strDesign, lambdas = newLambdas, confidence = NULL, trace = F), silent = !trace)
        if ("try-error" %in% class(fit)) {
          if (trace) cat("\nError in STRmodel. Reverting lambda coefficients.")
          fit <- oldFit
          newLambdas <- oldLambdas
        }
      }
      newData <- fit$output$predictors[[i]]$data + fit$output$random$data
      newPredictors <- fit$input$predictors[i]
      startLambdas <- newLambdas[i]
      fit_ <- STR_(
        data = newData,
        predictors = newPredictors,
        confidence = confidence,
        lambdas = startLambdas,
        pattern = extractPattern(newPredictors),
        nFold = nFold,
        reltol = reltol,
        gapCV = gapCV,
        solver = solver,
        trace = F,
        ratioGap = ratioGap
      )
      oldNorm <- sum(abs(newLambdas[[i]]$lambdas))
      newNorm <- sum(abs(fit_$input$lambdas[[1]]$lambdas))
      if (newNorm < oldNorm) {
        next
      } else {
        oldLambdas <- newLambdas
        oldFit <- fit
        newLambdas[[i]]$lambdas <- pmin(fit_$input$lambdas[[1]]$lambdas, maxLambda)
        fit <- NULL
        if (trace) {
          dummy <- lapply(newLambdas, FUN = function(p) {
            cat(format(p$lambdas, digits = 4, scientific = T, width = 9))
            cat("  ")
          })
        }
      }
    }
  }
  return(newLambdas)
}