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 <- ?

Answer 1

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

Answer 2

"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.