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 search of lambda parameters.
#' @return A structure containing input and output data.
#' It is an \strong{S3} class \code{STR}, which is a list with the following components:
#' \itemize{
#' \item \strong{output} -- contains decomposed data. It is a list of three components:
#' \itemize{
#' \item \strong{predictors} -- a list of components where each component
#' corresponds to the input predictor. Every such component is a list containing the following:
#' \itemize{
#' \item \strong{data} -- fit/forecast for the corresponding predictor (trend, seasonal component, flexible or seasonal predictor).
#' \item \strong{beta} -- beta coefficients of the fit of the coresponding predictor.
#' \item \strong{lower} -- optional (if requested) matrix of lower bounds of confidence intervals.
#' \item \strong{upper} -- optional (if requested) matrix of upper bounds of confidence intervals.
#' }
#' \item \strong{random} -- a list with one component \strong{data}, which contains residuals of the model fit.
#' \item \strong{forecast} -- a list with two components:
#' \itemize{
#' \item \strong{data} -- fit/forecast for the model.
#' \item \strong{beta} -- beta coefficients of the fit.
#' \item \strong{lower} -- optional (if requested) matrix of lower bounds of confidence intervals.
#' \item \strong{upper} -- optional (if requested) matrix of upper bounds of confidence intervals.
#' }
#' }
#' \item \strong{input} -- input parameters and lambdas used for final calculations.
#' \itemize{
#' \item \strong{data} -- input data.
#' \item \strong{predictors} - input predictors.
#' \item \strong{lambdas} -- smoothing parameters used for final calculations (same as input lambdas for STR method).
#' }
#' \item \strong{cvMSE} -- optional cross validated (leave one out) Mean Squared Error.
#' \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.
#' \item \strong{nFold} -- the input \code{nFold} parameter.
#' \item \strong{gapCV} -- the input \code{gapCV} parameter.
#' \item \strong{method} -- contains strings \code{"STR"} or \code{"RSTR"} depending on used method.
#' }
#' @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)
}