ForecastingMethodsEvaluation/SSAForecast.R at master · Remit/ForecastingMethodsEvaluation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
require(forecast)
require(Rssa)

# Function to derive the forecasts using SSA model
ssa.forecast <- function(train.timeseries, pred.steps) {
  model.fitting.start <- Sys.time()
  ssa.model = ssa(train.timeseries$series)
  # Heuristic to determine the groups using correlations
  wcor.matrix <- wcor(ssa.model, groups = 1:20)
  unnaceptable.wcor.pairs <- (wcor.matrix < 0.95) & (wcor.matrix > 0.01) # To folter approximate zeros
  # We will be interested in top-left quadrant of pairs with FALSE values (~ correlations close to 1) using symmetry
  unnaceptable.wcor.pairs.cols <- split(unnaceptable.wcor.pairs,
                                        rep(1:ncol(unnaceptable.wcor.pairs),
                                            each = nrow(unnaceptable.wcor.pairs)))

  contains.true <- function(vector.data) {
    index <- which(vector.data == TRUE)
    if(length(index) != 0) {
      index <- min(index)
    } else {
      index <- -1
    }

    return(index)
  }

  indices <- lapply(unnaceptable.wcor.pairs.cols, contains.true)

  indices.pairs.df <- data.frame(col.index = c(1:length(indices)),
                                 row.index = unlist(indices))
  indices.pairs.df.filtered <- indices.pairs.df[indices.pairs.df$row.index > 0, ]
  col.index.max <- min(indices.pairs.df.filtered[indices.pairs.df.filtered$row.index > indices.pairs.df.filtered$col.index, "row.index"]) - 1
  row.index.max <- min(indices.pairs.df.filtered[indices.pairs.df.filtered$row.index < indices.pairs.df.filtered$col.index, "col.index"]) - 1

  group.max <- 10 #default value
  if(col.index.max == row.index.max) {
    group.max <- col.index.max
  }

  ssa.forecast <- forecast(ssa.model,
                           groups = list(1:group.max),
                           method = "recurrent",
                           bootstrap = TRUE,
                           len = pred.steps,
                           R = 1000)

  mean.fc <- ssa.forecast$mean
  sd.fc <- sd(mean.fc)
  lb.80 <- mean.fc - 1.28 * sd.fc
  lb.95 <- mean.fc - 1.96 * sd.fc
  lower <- data.frame(lb.80 = lb.80,
                      lb.95 = lb.95)

  ub.80 <- mean.fc + 1.28 * sd.fc
  ub.95 <- mean.fc + 1.96 * sd.fc
  upper <- data.frame(ub.80 = ub.80,
                      ub.95 = ub.95)

  ssa.forecast$mean <- mean.fc
  ssa.forecast$lower <- lower
  ssa.forecast$upper <- upper
  model.fitting.end <- Sys.time()
  ssa.forecast$model.fitting.duration <- difftime(model.fitting.end, model.fitting.start, units = "secs")
  ssa.forecast$parameters.selection.duration = 0

  return(ssa.forecast)
}