How to create a prediction interval with the fact that the residuals follow a specific distribution (in python)

Question

I am looking at a software development pipeline where I am predicting the lead time of different products flowing through the pipeline.

After applying a boxcox transformation on the lead time (target variable) and creating a XGBoost regressor model I can see that the residuals follow a t-locationScale distribution.

So now I looked at this guide which describes a method to create a prediction interval for any regression model assuming that the residuals are normally distributed. https://qucit.com/a-simple-technique-to-estimate-prediction-intervals-for-any-regression-model_en/

But I tried to tweak it to my distribution.

So a t-locationScale distribution has a $\sigma$, $\mu$ and $\nu$ parameter. The variance is only defined for $\nu>2$. My specific distribution has $\nu = 2.56$ and $\mu = 0.04$, $\sigma = 0.97$ So I could take the 95% interval of this distribution and say that for any $\hat{y}$, the prediction interval is the 95% interval of the residual distribution.

But I want to take into consideration that the prediction interval should change with different inputs. I created a regressor model, which I trained and then made predictions using the validation set. I then took the square of the error and trained an additional error model on this data. Such that the error model could predict the variance of the residuals distribution.

  xgb = XGBoostRegressor()
  xgb.fit(X_train,y_train)
  y_hat = xgb.predict(X_val)
  val_error = (y_hat-y_val)**2

  xgb_error = XGBoostRegressor()
  xgb_error.fit(X_val, val_error)

  variance_hat_residuals = xgb_error.predict(X_test)

The relationship between variance and $\sigma$ and $\nu$ for a t-locationScale distribution is

var = $\sigma^2 *\frac{\nu}{\nu-2}$

Now here is where I make an assumption which I am not sure makes sense.

I assume that the degrees of freedom $\nu$ is the same as for all residuals, $\nu = 2.56$ and then I solve for $\sigma$ through the following.

$\hat{\sigma} = \sqrt{\frac{\hat{var}*(\nu-2)}{\nu}}$

And estimate the lower and upper quantiles from this distribution.

 residual_distribution_lower_quantile = scipy.stats.t.ppf(q = 0.025, df = 2.56, scale = sigma)
 residual_distribution_upper_quantile = scipy.stats.t.ppf(q = 0.0975, df = 2.56, scale = sigma)

I then predict the lead time $\hat{y}$ and say that the mean of the distribution is $\hat{y}$

   pred = xgb.prediction(X_test)
   lower_interval = pred + residual_distribution_lower_quantile
   upper_interval = pred + residual_distribution_upper_quantile

Does it make sense to make the claim of $\nu$ is static? My score for the prediction interval is now $81\%$ since I am clearly simplifying the problem.

Any suggestions for improving my method?