Get Prediction Intervals from hts package R - Hierarchical and Grouped Time Series in R

Question

I'm still pretty new to R and want to calculate the Prediction Intervals of both of my Time Series. The two datasets are already prepared as below-mentioned. I'm not sure how to get the specific Prediction Interval for my two ARIMA models. The most important to me is to get the Intervals from the top levels, but it would also be nice to know how to calculate them for any other level of my grouped Time Series.

Thank you in advance for the help!

This is my dataset so far. I already prepared it. (The second dataset has the exact same structure, but different values)

> tsib_data
# A tsibble: 448 x 5 [1Y]
# Key:       char, sex, age [32]
    year value char   sex   age  
   <int> <dbl> <chr> <chr> <chr>
 1  2006    62 M     m      5 
 2  2007    57 M     m      5 
 3  2008    80 M     m      5 
 4  2009    88 M     m      5 
 5  2010    90 M     m      5 
 6  2011    86 M     m      5 
 7  2012    87 M     m      5 
 8  2013    83 M     m      5 
 9  2014    99 M     m      5 
10  2015   103 M     m      5 
# ℹ 438 more rows
# ℹ Use `print(n = ...)` to see more rows

With the hts package from R I already found out that the following ARIMA models performed best. For them I want to calculate the Prediction Intervals.

#My best performing ARIMA models

#Model 1
forecast1 <- hts::forecast.gts(train, h = 40, method = "comb", weights = "mint", fmethod = "arima", nonnegative = TRUE, lambda = 0)


#Model 2
forecast2 <- forecast.gts(train, h = 40, method = "comb", weights = "ols", fmethod= "arima", nonnegative= TRUE, lambda=NULL)

Now I've already seen the following approach to calculate a Prediction Interval with ARIMA models from the hts/fable package:

#The approach I've seen already

fcsts <- data %>% 
  aggregate_key(sex * age, value = sum(value)) %>%
  model(arima = ARIMA(value)) %>%
  mutate(mint = min_trace(arima)) %>%
  forecast(h = 32)

fcsts %>% 
  hilo(level=95) %>%
  filter(is_aggregated(sex) & is_aggregated(age))

In my case I would then need to add "char":

# Adding column "char" to get the top level prediction interval

fcsts <- data %>% 
  aggregate_key(sex * age * char, value = sum(value)) %>%
  model(arima = ARIMA(value)) %>%
  mutate(mint = min_trace(arima)) %>%
  forecast(h = 40)

fcsts %>% 
  hilo(level=95) %>%
  filter(is_aggregated(sex) & is_aggregated(age) & is_aggregated(char))

I dont know if this approach is giving me the specific Prediction Intervals for my two ARIMA models. How do I need to change the code to get the Intervals I am looking for?

Can anyone help me with this?