Passing different forecasting method to hierarchical time series forecast in R?

Question

I have a hierarchical time series, the bottom level series of which all exhibit intermittent demand. It seems advantageous to use Hyndman's HTS package for optimal combination within the hierarchy. It also seems advantageous to use Kourentzes' MAPA package for multiple aggregation prediction of the intermittent demand. In essence, I want to do something like:

forecast(my_hts, method='comb', fmethod='MAPA')

However, it is unclear to me if / how I can combine the two, since forecast.gts() only accepts fmethod=c("ets", "arima", "rw").

Is there a clever way to pass different forecasting methods to forecast.gts() without having to tear up the code?

Example to clarify what I mean:

library(hts)
library(MAPA)
set.seed(1)

#note intermittent demand of bottom level time series
x <- ts(rpois(365, lambda=0.05), frequency=365, start=2014)
y <- ts(rpois(365, lambda=0.07), frequency=365, start=2014)

#it's easy to make a MAPA forecast for the top-level time series
#but this isn't an optimal hierarchical forecast
mapasimple(x+y)

#it's also easy to make this a HTS and make an optimal hierarchical forecast
#but now I cannot use MAPA
z <- hts(data.frame(x,y)))
z_arima <- forecast(z, fmethod="arima")
z_rw <- forecast(z, fmethod="rw")
z_ets <- forecast(z, fmethod="ets")

#z_MAPA <- ?

Thank you for adding a reproducible example. We can now try to migrate this to [SO] where programming questions belong & are readily answered. — gung - Reinstate Monica, Feb 23 '15 at 04:00
I obviously defer to your moderation, but as both authors of the mentioned packages are (to my knowledge) active in CrossValidated, I thought this was a better place to post than general SO. — user1569317, Feb 23 '15 at 04:08
That was a reasonable guess, @user1569317, & this is a tricky & somewhat contentious issue. There is in fact a large & active group of R users on SO (more than here, I believe). My criterion is what the OP needs to be explained; if it is a statistical concept, the Q belongs here, if it is how the code works, the Q belongs on SO. I read your Q as "Is there a clever way to pass different forecasting methods to forecast.gts() without having to tear up the code?", not as 'how does forecasting work'. Hence I think SO is more suitable. (I also believe you will get a better / faster A there.) — gung - Reinstate Monica, Feb 23 '15 at 04:30

score 5 · Accepted Answer · edited Feb 24 '15 at 11:29

I'm posting because after a closer look at the hts documentation (insert well-deserved RTFM here), I think I found a work-around using the combinef() function from hts, which can be used to optimally combine forecasts outside of the forecast.gts() environment. I'll leave it up for a while before accepting as an answer so others can tell me if I'm wrong.

fh <- 8

library(hts)
library(MAPA)
set.seed(1)

x <- ts(rpois(365, lambda=0.05), frequency=365, start=2014)
y <- ts(rpois(365, lambda=0.07), frequency=365, start=2014)

my_hts <- hts(data.frame(x,y))

ally <- aggts(my_hts)

allf <- matrix(NA, nrow = fh, ncol = ncol(ally))

for(i in 1:ncol(ally)){
    allf[,i] <- mapafor(ally[,i], 
                        mapaest(ally[,i], outplot=0), 
                        fh = fh, 
                        outplot=0)$outfor
}

allf <- ts(allf)

y.f <- combinef(allf, my_hts$nodes, weights=NULL, keep="bottom")

#here's what the non-reconciled, bottom-level MAPA forecasts look like
print(allf[1,1:2])

 Series 1   Series 2
1 0.1343304 0.06032574

#here's what the reconciled MAPA bottom-level forecasts look like
#notice that they're different
print(y.f[1,])

[1] 0.06030926 0.07402938

This looks completely correct to me. (I took the liberty of correcting a few typos.) It should work in your setup, with time-invariant forecasts. If your forecasts vary over time (e.g., with trends, seasonality or causal factors), then you could run into the problem I mentioned above: your "optimal" forecasts may turn negative. Otherwise fine, +1! — Stephan Kolassa, Feb 24 '15 at 11:30

score 3 · Answer 2 · answered Feb 23 '15 at 07:37

"Is there a clever way to pass different forecasting methods to forecast.gts() without having to tear up the code?" - Almost certainly not. forecast.gts() does not allow you to plug in forecast methods along the lines of the family parameter to glm() and similar.

It will likely be easiest to run the MAPA forecasts and then reimplement the optimal combination yourself. It's really not all that hard, I did it myself a few times before the hts package came along. Look at the references in the hts package; most of them are available as preprints on Rob Hyndman's website.

One key problem in combining the optimal combination approach with intermittent demand forecasts is that the optimal combination may well yield negative forecasts. These may be "optimal" in the GLS sense used here, but they make no sense for demands. My solution was to use pcls() in the mgcv package for the actual reconciliation, in order to constrain the solution (i.e., the bottom level reconciled forecasts) to be nonnegative.

Thanks for the answer. Mind taking a look at my proposed solution below? If it is not valid, I will accept your answer. — user1569317, Feb 23 '15 at 20:25

Nikos Kourentzes · Answer 3 · 2015-02-24T12:49:57.807

I have run successfully in the past small experiments using MAPA forecasts and combinef. Stephan's point about the negative forecasts remains an issue. You might be interested to look into iMAPA in the tsintermittent package that is specifically designed to produce forecasts for intermittent demand time series. In contrast to MAPA that one implements Croston's method, SBA and exponential smoothing to choose from in the various temporal aggregation levels.

I would be very interested in learning more about your results. Perhaps a relevant research was presented in last year's ISF by E. Spiliotis et al. Examining the effect of temporal aggregation on forecasting accuracy for hierarchical energy consumption time series.

Thanks for the heads-up - both iMAPA and the presentation look very promising. If I come up with something interesting, I'll email results. — user1569317, Feb 24 '15 at 14:56

Passing different forecasting method to hierarchical time series forecast in R?

3 Answers3

Linked