Q: Is there a good work-around broken user-defined transform handling in fable?

Following Forecasting Principles and Practice, Ch. 13.3, Ensuring forecasts stay within limits on a dataset, I found strange behavior around user-defined transforms, back transforms, and forecasts. In particular, for the scaled logit transform. See below.

Is there a good hack/fix/work-around for this? I can always manually transform and backtransform the data myself, before passing it to model(). But it’d be nice if someone knew of way to avoid this.

library(fpp3)

y <- c(145.1, 27.3, 37.4, 3.2, 3.8, 106.9, 4.3, 43.2, 16.2, 251.9, 1.1, 20.8, 
       17.7, 47.1, 34.2, 85.3, 61.5, 34.8, 216.3, 28.9, 4.8, 5.6, 104.3, 10.7, 
       120.8, 37.2, 435.5, 22.3, 23.9, 3.6, 10.3, 13.4, 19.9, 4.1, 5.9, 4.1, 
       756.3, 68.5, 218, 32, 227.1, 79.9, 141.2, 24.5, 179.6, 39.3, 391.4, 12.2, 
       9.2, 13.1, 8, 15.9, 2.3, 24.1, 14.7, 1.4, 4.8, 6.6, 7.8, 4.4, 0.1, 0.1, 0.1, 
       0.1, 0.2, 0.5, 0.1, 0.1, 0.1, 0.1, 0.2, 0.1, 0.1, 0.2, 0.3, 0.1, 0.1, 9.4, 2.3, 
       30.5, 1.3, 0.2, 0.4, 0.1, 0.8, 2.2, 0.1, 2.7, 0.1, 5.4, 0.2, 0.1, 0.6, 0.1, 
       0.1, 27.9, 12, 21.8, 5.4, 18217.6, 22806.3, 46311.4, 12808, 10253.8, 
       11275.2, 8353.5, 11929.1, 12147.8, 1385.8, 17915.9, 343, 4152.8, 5520.9)

# upper and lower bounds of data
max(y) # 46311.4
min(y) # 0.1

# construct data
tibble(idx=1:113,y=y) %>% tsibble(index=idx) -> dat

# define scaled logit
# via Forecasting Priniples and Practice, Ch. 13.3, Ensuring forecasts stay
# within limits
# https://otexts.com/fpp3/limits.html
scaled_logit <- function(x, lower = 0, upper = 1) {
  log((x - lower) / (upper - x))
}
inv_scaled_logit <- function(x, lower = 0, upper = 1) {
  (upper - lower) * exp(x) / (1 + exp(x)) + lower
}
my_scaled_logit <- new_transformation(
  scaled_logit, inv_scaled_logit)

# forecasting, setting bounds so that all data between
dat %>%
  model(
    y.lgt.arima=ARIMA(
      my_scaled_logit(y,lower=0,upper=50000)~pdq(d=0)+PDQ(0,1,0,period=54)
      )
    ) ->
  fit

# fits stay within bounds
fit %>% fitted() %>% pull(.fitted) %>% max() # 29475.22
fit %>% fitted() %>% pull(.fitted) %>% min() # 0.1744393

# forecast goes outside of bounds
fit %>% forecast(h=100) -> fcs
fcs %>% pull(.mean) %>% max() # 82827.8
fcs %>% pull(.mean) %>% min() # -119000

fcs %>% autoplot(level=NULL)+autolayer(dat)