How can I get the confidence interval of hierarchical forecasts reconciled with hts::combinef?

Question

I have hierarchical data that I want to forecast, reconcile, and then plot each model's forecast.

UPDATE: In my attempt to make a minimum reproducible example I made it too simple and didnt really ask the question I'm trying to ask. I changed the modeling function from auto.arima to a custom function that goes through a number of models. The reason for this is I dont just want to model / plot the best fit model but a number of them. Not shown here is that auto.arima is actually used in the modeling function to get the few best arima models. So in summary I want to model and plot a number of hierarchical series using a number of models and I want to be able to visualize them and their prediction intervals individually after reconciliation.

Using the example in the function vignette I can forecast each series in the hierarchy and plot it like this:

library(forecast)
library(hts)
h <- 10
ally <- aggts(htseg1)
## changed this from allf to three different tables showing I need to 
## be able to check different models, not just the best fit.
allf_single <- matrix(NA, nrow=h, ncol=ncol(ally)) 
allf_dub_add <- matrix(NA, nrow=h, ncol=ncol(ally))
allf_dub_add_damp <- matrix(NA, nrow=h, ncol=ncol(ally))
all_mod <- list()

## Adding custom function to show I need to look at multiple models per
## series. 
mod_func <- function(time_series, h) {
  single_es <- ses(time_series, h=h)
  double_es_add <- holt(time_series, damped=FALSE, exponential=FALSE, h=h)
  double_es_add_damp <- holt(time_series, damped=TRUE, exponential=FALSE, h=h)
  
  model_list <- list(single=single_es, 
                          double_add=double_es_add, 
                          double_add_damp=double_es_add_damp)
  return(model_list)
}

#forecast each series
for(i in 1:ncol(ally)) {
  series <- colnames(ally)[i]
  ## changed this from auto.arima to custom function
  all_mod[[series]] <- mod_func(ally[ , i], h=h) 
  ## capturing predictions for each model and each series
  ## (only single shown here)
  allf_single[ , i] <- all_mod[[series]][["single"]]$mean
  allf_dub_add[ , i] <- all_mod[[series]][["double_add"]]$mean
  allf_dub_add_damp[ , i] <- all_mod[[series]][["double_add_damp"]]$mean
}

#plot multiple models for single series
autoplot(all_mod[["Total"]][["single"]])
autoplot(all_mod[["Total"]][["double_add"]])
autoplot(all_mod[["Total"]][["double_add_damp"]])

I can then reconcile the hierarchy, substitute in the new forecast values, and plot an individual series like this:

#reconcile the hierarchy
allf_dub_add_damp <- ts(allf_dub_add_damp, start=51)
y.f <- combinef(allf_dub_add_damp, get_nodes(htseg1), keep="all", algorithms="lu")

#replace forecast values with reconciled values
for(i in 1:length(all_mod[["Total"]][["double_add_damp"]][["mean"]])) {
  all_mod[["Total"]][["double_add_damp"]][["mean"]][i] <- y.f[ , 1][i]
}

#plot single reconciled series
autoplot(all_mod[["Total"]][["double_add_damp"]])

Where I'm stuck is how to adjust the confidence interval of the model post reconciliation. I know if I can calculate the confidence interval I can substitute it into the forecast model object in the upper and lower slots like the forecast values, but I'm not sure how to calculate the confidence interval after reconciliation. The only thing I have found that seems relevant to this question is this post, but there is no answer to the question.

Any help would be greatly appreciated!

Answer 1

The fable package replaces the hts package and produces prediction intervals. It is also much easier to handle the aggregation structure. Here is some code using the same example as in your question (updated to include multiple models).

library(tsibble)
library(feasts)
library(fable)
library(dplyr)

df <- as_tsibble(hts::htseg1$bts) %>%
  mutate(
    level1 = substr(key, 1, 1),
    level2 = substr(key, 2, 2)
  ) %>%
  as_tsibble(index=index, key=c(level1,level2)) %>%
  select(-key) %>%
  aggregate_key(level1/level2, value = sum(value))

fc <- df %>%
  model(
    ses = ETS(value ~ trend("N") + season("N")),
    holt = ETS(value ~ trend("A") + season("N")),
    arima = ARIMA(value)
  ) %>%
  reconcile(
    ses = min_trace(ses),
    holt = min_trace(holt),
    arima = min_trace(arima)
  ) %>%
  forecast(h=10)
fc %>%
  filter(is_aggregated(level1)) %>%
  autoplot(df)

fc %>%
  mutate(PI = hilo(value, level=95))
#> # A fable: 240 x 7 [1Y]
#> # Key:     level1, level2, .model [24]
#>    level1 level2 .model index         value .mean                      PI
#>    <chr*> <chr*> <chr>  <dbl>        <dist> <dbl>                  <hilo>
#>  1 A      A      ses     2002 N(9.2, 0.031)  9.15 [8.807747,  9.497520]95
#>  2 A      A      ses     2003 N(9.2, 0.062)  9.15 [8.664913,  9.640354]95
#>  3 A      A      ses     2004 N(9.2, 0.093)  9.15 [8.555308,  9.749959]95
#>  4 A      A      ses     2005  N(9.2, 0.12)  9.15 [8.462905,  9.842362]95
#>  5 A      A      ses     2006  N(9.2, 0.15)  9.15 [8.381495,  9.923773]95
#>  6 A      A      ses     2007  N(9.2, 0.19)  9.15 [8.307893,  9.997374]95
#>  7 A      A      ses     2008  N(9.2, 0.22)  9.15 [8.240208, 10.065059]95
#>  8 A      A      ses     2009  N(9.2, 0.25)  9.15 [8.177208, 10.128059]95
#>  9 A      A      ses     2010  N(9.2, 0.28)  9.15 [8.118037, 10.187231]95
#> 10 A      A      ses     2011  N(9.2, 0.31)  9.15 [8.062070, 10.243198]95
#> # … with 230 more rows

^{Created on 2020-10-24 by the reprex package (v0.3.0)}