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

forecasting, hierarchical, r, time-series

Solution

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

Problem

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

Original source