This question already has an answer here:

I am doing forecasts for a lot of time series. These consist of forecasts for specific items in specific regions, where all regions belong to a country. So this is a hierarchial series.

  • Product 1 - Country 1
    1. product 1 - country 1 - region 1-1
    2. product 1 - country 1 - region 1-2
  • Product 1 - Country 2
    1. Product 1 - Country 2 - region 2-1
    2. ........
  • ........

And then there is of course also the "Total" of all.

My matrix cmat contains all bottom level data, i.e. for all products and regions. Every column is a time series. The name of a column is for example: 000214ABCab. The numbers are the product ID, the next 3 letters are the country, and the next two letters are the region. I want forecasts on this level, but also on the product-country level.

First I was struggling with the performance, as I have a lot of data. In the example below I only use a subset of 1000 time series.

subs = cmat[,1:1000]
hh = hts(subs, characters = c(maxlength+3,2))
ally = aggts(hh)

#make list of data tables with time series 
subs2 = as.data.table(t(ally))
subs2$CODE = unlist(hh$labels)
datlist = subs2[, list(list(.SD)),by=CODE]$V1
	setattr(datlist, 'names', subs2$CODE) #code is the identifier of products/country/region

s1 = Sys.time()
fcast = NULL
for(i in 1:nrow(subs2)){
    fc = data.table(pmax(forecast(stlf(ts(matrix(datlist[[i]]),frequency=12)),h=12)$mean,0) )
    names(fc) = names(datlist)[i]

    fcast = rbind(fcast, t(fc))
s2 = Sys.time()

fcasts = t(fcast)
y.f = combinef(fcasts, get_nodes(hh),keep="all")

s3 = Sys.time()

I am also surprised by how fast combinef is. However, I only want positive or zero forecasts. Preferably also integer forecasts. I already do pmax(forecast, 0) in the loop. But by using combinef, I lose this again.

If I do now again pmax then I will also lose the functionality of the combinef, right? So how can I make the forecasts positive? Also, for me it is not necessary to have a "total forecast" at all, but I don't know if I can omit this?


marked as duplicate by Stephan Kolassa, mdewey, kjetil b halvorsen, jbowman, Peter Flom Mar 2 '18 at 13:33

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.


You can't ensure positivity and sum consistency of hierarchical forecasts if you use forecast::combinef().

What I personally found useful was to set up the summation matrix (see the original publication by Hyndman et al., 2011) and then solve the relevant (weighted) least squares problem with additional nonnegativity constraints. This will give you nonnegative and sum consistent forecasts. I have found repeatedly that this approach still results in better forecasts on all levels in the hierarchy.

One possible tool that solves (weighted) least squares with linear constraints is the pcls() function in the mgcv package. However, this is of course also not optimized to leverage the specific structure of potential forecast hierarchy matrices, so your performance may be significantly worse than if you use combinef().

  • $\begingroup$ Thank you for your answer! You have experience with using this method? $\endgroup$ – pk_22 Feb 26 '18 at 11:44
  • $\begingroup$ Yes, as I wrote. $\endgroup$ – Stephan Kolassa Feb 26 '18 at 12:17
  • $\begingroup$ Thank you. My data is on spare parts mainly, do you consider stlf() as a good method for this? $\endgroup$ – pk_22 Feb 26 '18 at 12:22
  • $\begingroup$ Spare parts demand is often intermittent-time-series, where seasonality is not overly evident. You may want to look at the MAPA algorithm by Kourentzes et al. $\endgroup$ – Stephan Kolassa Feb 26 '18 at 14:30
  • $\begingroup$ Thank you. I did, but I'm not satisfied with it, see my results in the update on my post. Or do I miss anything when using mapa? $\endgroup$ – pk_22 Feb 27 '18 at 8:24

Not the answer you're looking for? Browse other questions tagged or ask your own question.