Forecasting Hierarchical Time Series

Question

I am trying to perform 6 months forecasting over production data for three power plants, I built my data as an hts object that has 3 levels. However, when I am performing the forecast function and then try to see the accuracy using test data I get the following error: "Error in x - fcasts: non-conformable arrays"

Furthermore, when I try to apply the "arima" as a forecasting method on the hts object I get the following (the warning message is repeated 9 times, as I have 9 time series in the hts object):

forecasts <- forecast(data,h = 6 , method = "bu" , fmethod = "arima")

I used the following instructions to get the hts object:

and the data has the following structure: I am not sure where I am going wrong. Anyone can help with some thoughts??

Thank you!

The data:

structure(list(LarGroup1 = c(188.3, 187.2, 94.7, 109.2, 202.7, 
146.6, 121.9, 151.3, 111.1, 103.4, 188.1, 168.1, 233.9, 230.7, 
187.1, 0, 98.9, 173.5, 149.4, 168.6, 4.7, 14.8, 91.8, 166.5, 
170.5, 123.6, 85.2, 64.4), LarGroup2 = c(159.1, 127.7, 210.3, 
199.8, 113, 143.4, 144.5, 83.8, 41.6, 35.1, 95.2, 178.2, 241.1, 
236.4, 181.9, 194.3, 196.1, 92.4, 154.6, 78.9, 35.7, 0, 74.5, 
75.1, 140, 142.5, 3.8, 17.5), RibGroup1 = c(49.4, 102.4, 50.8, 
118.8, 108.4, 139.5, 121.7, 69.6, 53.4, 28, 113.3, 96.3, 70.8, 
124.4, 54.4, 128.7, 63.3, 2.1, 41.3, 0.4, 0.6, 0, 5.4, 57.9, 
9.9, 30, 221, 167.2), RibGroup2 = c(32.7, 32, 98.1, 6.3, 85.5, 
96.6, 41.1, 44.9, 50.4, 27.3, 0, 45.4, 199.1, 179.2, 86.1, 0, 
58.4, 43.3, 41.8, 42.1, 22.1, 11.8, 71.8, 112, 204.1, 40.9, 24.5, 
210.9), RibGroup3 = c(90.8, 15.4, 10.5, 124.4, 33.9, 8.4, 38.3, 
56.9, 13.5, 0, 32.6, 132.8, 160.7, 168.7, 60.7, 131.9, 110.8, 
29.2, 131.3, 62.1, 6.1, 0, 0, 3.4, 23.9, 192.7, 165.5, 0), SinGroup1 = c(235.2, 
225.4, 226.1, 234.4, 222.1, 232.3, 233.4, 201.9, 195.3, 209.4, 
233.6, 223.6, 222.2, 232, 224, 149.8, 201.6, 220.2, 203.1, 212.1, 
71.9, 82.3, 183.2, 210.6, 198.6, 230.8, 218, 163.2), SinGroup2 = c(233.4, 
225.6, 227, 51.6, 76, 230.7, 233.1, 202.7, 200.2, 207.2, 228.4, 
226.2, 183.9, 230.4, 222.3, 227.7, 177.9, 152, 218.6, 210.6, 
80.9, 63.2, 188.1, 209.5, 233.2, 210.1, 226.5, 200.5), SinGroup3 = c(233.2, 
188.5, 226.9, 234.7, 222.8, 234.6, 220.6, 156.4, 209.2, 218.7, 
232.9, 226.1, 215.4, 231, 222.7, 222.7, 183.7, 203.8, 216.8, 
112, 0, 39.6, 180.8, 203.6, 221.1, 228.9, 202.8, 186.7), SinGroup4 = c(218, 
215.5, 226.8, 235.6, 223.6, 234.8, 234.9, 69.3, 192, 207.8, 235.2, 
217.2, 235.1, 231.8, 223.5, 230.5, 225.6, 220.1, 220, 211.9, 
114.8, 44.5, 158.5, 206.3, 231.8, 179, 225.3, 198.6)), class = "data.frame", row.names = c(NA, 
-28L))

I don't know much about hierarchical time-series. Try reading otexts.org/fpp2/hierarchical.html about them. This online book is written by Rob Hyndman and explains step-by-step how to run a hierarchical forecast using hts. — hmhensen
When I create a data variable using your pasted data.frame, then call forecasts <- forecast(data, h=6, method="bu"), I get a different error: "Error in ets(object, lambda = lambda, biasadj = biasadj, allow.multiplicative.trend = allow.multiplicative.trend, : y should be a univariate time series". Is there more detail we're missing? — Ken Williams
@hmhensen I am following the chapters of this book, and I restructured my data based on the data set used in the example but still, I am not getting the results that I wish for. — salam abdeen
My main problem is applying the accuracy function, since what I need to do is comparing several method and choose the most fit one to the data — salam abdeen

Rob Hyndman Rob Hyndman · Accepted Answer · 2018-12-20T23:00:14

In the accuracy function, you need to include test data, not training data. You ask for 6 steps ahead, but your test data only consists of 4 time periods.
The seasonal differencing error suggests you are using an old version of the forecast package. Please update your packages.

The following code works using current CRAN packages (forecast v8.4, hts v

library(hts)
Production_data <- data.frame(
  LarGroup1 = c(
    188.3, 187.2, 94.7, 109.2, 202.7,
    146.6, 121.9, 151.3, 111.1, 103.4, 188.1, 168.1, 233.9, 230.7,
    187.1, 0, 98.9, 173.5, 149.4, 168.6, 4.7, 14.8, 91.8, 166.5,
    170.5, 123.6, 85.2, 64.4
  ), LarGroup2 = c(
    159.1, 127.7, 210.3,
    199.8, 113, 143.4, 144.5, 83.8, 41.6, 35.1, 95.2, 178.2, 241.1,
    236.4, 181.9, 194.3, 196.1, 92.4, 154.6, 78.9, 35.7, 0, 74.5,
    75.1, 140, 142.5, 3.8, 17.5
  ), RibGroup1 = c(
    49.4, 102.4, 50.8,
    118.8, 108.4, 139.5, 121.7, 69.6, 53.4, 28, 113.3, 96.3, 70.8,
    124.4, 54.4, 128.7, 63.3, 2.1, 41.3, 0.4, 0.6, 0, 5.4, 57.9,
    9.9, 30, 221, 167.2
  ), RibGroup2 = c(
    32.7, 32, 98.1, 6.3, 85.5,
    96.6, 41.1, 44.9, 50.4, 27.3, 0, 45.4, 199.1, 179.2, 86.1, 0,
    58.4, 43.3, 41.8, 42.1, 22.1, 11.8, 71.8, 112, 204.1, 40.9, 24.5,
    210.9
  ), RibGroup3 = c(
    90.8, 15.4, 10.5, 124.4, 33.9, 8.4, 38.3,
    56.9, 13.5, 0, 32.6, 132.8, 160.7, 168.7, 60.7, 131.9, 110.8,
    29.2, 131.3, 62.1, 6.1, 0, 0, 3.4, 23.9, 192.7, 165.5, 0
  ), SinGroup1 = c(
    235.2,
    225.4, 226.1, 234.4, 222.1, 232.3, 233.4, 201.9, 195.3, 209.4,
    233.6, 223.6, 222.2, 232, 224, 149.8, 201.6, 220.2, 203.1, 212.1,
    71.9, 82.3, 183.2, 210.6, 198.6, 230.8, 218, 163.2
  ), SinGroup2 = c(
    233.4,
    225.6, 227, 51.6, 76, 230.7, 233.1, 202.7, 200.2, 207.2, 228.4,
    226.2, 183.9, 230.4, 222.3, 227.7, 177.9, 152, 218.6, 210.6,
    80.9, 63.2, 188.1, 209.5, 233.2, 210.1, 226.5, 200.5
  ), SinGroup3 = c(
    233.2,
    188.5, 226.9, 234.7, 222.8, 234.6, 220.6, 156.4, 209.2, 218.7,
    232.9, 226.1, 215.4, 231, 222.7, 222.7, 183.7, 203.8, 216.8,
    112, 0, 39.6, 180.8, 203.6, 221.1, 228.9, 202.8, 186.7
  ), SinGroup4 = c(
    218,
    215.5, 226.8, 235.6, 223.6, 234.8, 234.9, 69.3, 192, 207.8, 235.2,
    217.2, 235.1, 231.8, 223.5, 230.5, 225.6, 220.1, 220, 211.9,
    114.8, 44.5, 158.5, 206.3, 231.8, 179, 225.3, 198.6
  )
)
Production_data_ts <- ts(Production_data, frequency = 12, start = c(2016, 7))
Production_data_hts <- hts(Production_data_ts, characters = c(3, 6))
data <- window(Production_data_hts, start = c(2016, 7), end = c(2018, 6))
test <- window(Production_data_hts, start = c(2018, 7), end = c(2018, 10))
forecasts <- forecast(data, h = 4, method = "bu")
accuracy(forecasts, test)

Forecasting Hierarchical Time Series

1 Answers