NGBoost multi-step ahead probabilistic forecasting for 96 time steps

-1 I am building a probabilistic energy forecasting model using NGBoost (NGBRegressor) to predict household residual power (residual_power_w) at 15-minute resolution. My goal is to forecast the full next 24 hours (96 steps) with uncertainty estimates (prediction intervals, CRPS) at each step. The problem: NGBoost only supports a scalar target — fit() accepts a 1D array. Unlike XGBoost which has multi_strategy='multi_output_tree' for native multi-output, NGBoost has no equivalent. What I have tried: The NGBoost maintainer on GitHub Discussion #243 confirmed the only current approach is one model per horizon step: "At the moment there isn't a simple way to predict all of [yt+1 ... yt+n] in one shot, so you have to use one model for t+1, another model for t+2, and so on up to t+n" So training 96 separate NGBoost models is the documented workaround. However this is very expensive since NGBoost is already slow compared to XGBoost. The MultivariateNormal(k=96) distribution exists in NGBoost but requires k*(k+3)/2 = 4,752 parameters per row — completely impractical for k=96. My question: Is there a better approach to get probabilistic 96-step ahead forecasts from NGBoost without training 96 separate models? Specifically I need: Prediction intervals at each of the 96 horizon steps CRPS evaluation per step Any alternative that keeps the probabilistic output (not just point forecasts) is welcome. Environment: Python 3, ngboost 0.4+, sklearn