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). ``` r 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) ``` ![](https://i.sstatic.net/EdPId.png) ``` r 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 ``` <sup>Created on 2020-10-24 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)</sup>