R/adjustSeasonalSeeds.R

In forecast: Forecasting Functions for Time Series and Linear Models

###############################################################################

# TBATS code
cutWTBATS <- function(use.beta, w.tilda.transpose, seasonal.periods, p=0, q=0) {
  mask.vector <- numeric(length(seasonal.periods))
  i <- length(seasonal.periods)
  while (i > 1) {
    for (j in 1:(i - 1)) {
      if ((seasonal.periods[i] %% seasonal.periods[j]) == 0) {
        mask.vector[j] <- 1
      }
    }
    i <- i - 1
  }

  w.pos.counter <- 1
  w.pos <- 1
  if (use.beta) {
    w.pos <- w.pos + 1
  }
  for (s in seasonal.periods) {
    if (mask.vector[w.pos.counter] == 1) {
      w.tilda.transpose <- w.tilda.transpose[, -((w.pos + 1):(w.pos + s))]
    } else if (mask.vector[w.pos.counter] < 0) {
      # Cut more than one off
      w.pos <- w.pos + s
      w.tilda.transpose <- w.tilda.transpose[, -c((w.pos + mask.vector[w.pos.counter] + 1):w.pos)]
      w.pos <- w.pos + mask.vector[w.pos.counter]
    } else {
      w.pos <- w.pos + s
      w.tilda.transpose <- w.tilda.transpose[, -w.pos]
      w.pos <- w.pos - 1
    }
    w.pos.counter <- w.pos.counter + 1
  }
  if ((p != 0) | (q != 0)) {
    end.cut <- ncol(w.tilda.transpose)
    start.cut <- end.cut - (p + q) + 1
    w.tilda.transpose <- w.tilda.transpose[, -c(start.cut:end.cut)]
  }

  return(list(matrix = w.tilda.transpose, mask.vector = mask.vector))
}

# BATS code below
#########
cutW <- function(use.beta, w.tilda.transpose, seasonal.periods, p=0, q=0) {
  mask.vector <- numeric(length(seasonal.periods))
  i <- length(seasonal.periods)
  while (i > 1) {
    for (j in 1:(i - 1)) {
      if ((seasonal.periods[i] %% seasonal.periods[j]) == 0) {
        mask.vector[j] <- 1
      }
    }
    i <- i - 1
  }
  if (length(seasonal.periods) > 1) {
    for (s in length(seasonal.periods):2) {
      for (j in (s - 1):1) {
        hcf <- findGCD(seasonal.periods[s], seasonal.periods[j])
        if (hcf != 1) {
          if ((mask.vector[s] != 1) && (mask.vector[j] != 1)) {
            mask.vector[s] <- hcf * -1
          }
        }
      }
    }
  }
  w.pos.counter <- 1
  w.pos <- 1
  if (use.beta) {
    w.pos <- w.pos + 1
  }
  for (s in seasonal.periods) {
    if (mask.vector[w.pos.counter] == 1) {
      w.tilda.transpose <- w.tilda.transpose[, -((w.pos + 1):(w.pos + s))]
    } else if (mask.vector[w.pos.counter] < 0) {
      # Cut more than one off
      w.pos <- w.pos + s
      w.tilda.transpose <- w.tilda.transpose[, -c((w.pos + mask.vector[w.pos.counter] + 1):w.pos)]
      w.pos <- w.pos + mask.vector[w.pos.counter]
    } else {
      w.pos <- w.pos + s
      w.tilda.transpose <- w.tilda.transpose[, -w.pos]
      w.pos <- w.pos - 1
    }
    w.pos.counter <- w.pos.counter + 1
  }
  if ((p != 0) | (q != 0)) {
    end.cut <- ncol(w.tilda.transpose)
    start.cut <- end.cut - (p + q) + 1
    w.tilda.transpose <- w.tilda.transpose[, -c(start.cut:end.cut)]
  }

  return(list(matrix = w.tilda.transpose, mask.vector = mask.vector))
}

calcSeasonalSeeds <- function(use.beta, coefs, seasonal.periods, mask.vector, p=0, q=0) {
  x.pos.counter <- 1
  sum.k <- 0
  if (use.beta) {
    x.pos <- 2
    new.x.nought <- matrix(coefs[1:2], nrow = 2, ncol = 1)
  } else {
    x.pos <- 1
    new.x.nought <- matrix(coefs[1], nrow = 1, ncol = 1)
  }
  x.pos.counter <- 1
  for (s in seasonal.periods) {
    if (mask.vector[x.pos.counter] == 1) {
      # Make a vector of zeros
      season <- matrix(0, nrow = s, ncol = 1)
      new.x.nought <- rbind(new.x.nought, season)
    } else if (mask.vector[x.pos.counter] < 0) {
      extract <- coefs[(x.pos + 1):(x.pos + s + mask.vector[x.pos.counter])]
      # print("extract:")
      # print(extract)
      # Find k
      k <- sum(extract)
      # update sum.k
      sum.k <- sum.k + k / s
      # create the current.periodicity vector
      current.periodicity <- extract - k / s
      current.periodicity <- matrix(current.periodicity, nrow = length(current.periodicity), ncol = 1)
      additional <- matrix(-k / s, nrow = (-1 * mask.vector[x.pos.counter]), ncol = 1)
      current.periodicity <- rbind(current.periodicity, additional)
      new.x.nought <- rbind(new.x.nought, current.periodicity)
      x.pos <- x.pos + s + mask.vector[x.pos.counter]
    } else {
      # Find k
      k <- sum(coefs[(x.pos + 1):(x.pos + s - 1)])
      # update sum.k
      sum.k <- sum.k + k / s
      # create the current.periodicity vector
      current.periodicity <- coefs[(x.pos + 1):(x.pos + s - 1)] - k / s
      current.periodicity <- c(current.periodicity, -k / s)
      current.periodicity <- matrix(current.periodicity, nrow = length(current.periodicity), ncol = 1)
      new.x.nought <- rbind(new.x.nought, current.periodicity)
      x.pos <- x.pos + s - 1
    }
    # Adjust L(t)
    x.pos.counter <- x.pos.counter + 1
  }
  # print(new.x.nought)
  # Lastly, get the arma error seed states, if they exist.
  if ((p != 0) | (q != 0)) {
    arma.seed.states <- numeric((p + q))
    arma.seed.states <- matrix(arma.seed.states, nrow = length(arma.seed.states), ncol = 1)
    # Final value of x.nought
    x.nought <- rbind(new.x.nought, arma.seed.states)
  } else {
    x.nought <- new.x.nought
  }
  return(x.nought)
}

findGCD <- function(larger, smaller) {
  remainder <- larger %% smaller
  if (remainder != 0) {
    return(findGCD(smaller, remainder))
  } else {
    return(smaller)
  }
}