THIS SEEMS TO BE AN OPEN PROBLEM. Let's look at some possible solutions and their drawbacks.
First, you propose this Yhat +- 1.96 * std(residuals). Let's put that in more standard statistics notation and then break down what the pieces mean.
$$
\hat y_i \pm 1.96\hat\sigma_{\hat\varepsilon}
$$
$\hat y_i$ is the point prediction.
$1.96$ comes from the $95\%$ bounds of a Gaussian distribution.
$\pm$ means that you go $1.96\hat\sigma_{\hat\varepsilon}$ above and $1.96\hat\sigma_{\hat\varepsilon}$ below the point prediction.
$\hat\sigma_{\hat\varepsilon}$ is the (constant) standard deviation of an implicit error term, estimated by the standard deviation of the residuals.
Thus, this formulation implies that the way to get the $95\%$ prediction interval is to use the point estimate $\hat y_i$ as the mean of a conditional Gaussian and the residual variance as the variance of that conditional Gaussian. We then take the $0.025$ and $0.975$ quantiles of this $N\left(\hat y_i, \hat\sigma^2_{\hat\varepsilon}\right)$ conditional distribution.
Under the assumption that there is a Gaussian error term with constant variance, this makes some sense. However, that is extremely restrictive. What about if the error term is not Gaussian? What about if the error term lacks a constant variance? What if the non-Gaussian error term lacks constant variance? Under this quite plausible deviations from the $N\left(\hat y_i, \hat\sigma^2_{\hat\varepsilon}\right)$ ideal, the prediction interval could wind up being terrible.
Second, you propose Yhat +- 1.96 * std(yhat). I interpret this to mean $\hat y_i\pm 1.96\sigma_{\hat y_i}$.
$\hat y_i$ is the point prediction (again).
$1.96$ comes from the $95\%$ bounds of a Gaussian distribution (again).
$\pm$ means that you go $1.96\hat\sigma_{\hat y_i}$ above and $1.96\hat\sigma_{\hat y_i}$ below the point prediction.
This is somewhat less restrictive than the earlier formulation, as it allows for each $i$ to have a different variance. However, the $1.96$ still implies a Gaussian conditional distribution that might not be reasonable. Additionally, estimating the variance at each $i$ is not entirely straightforward.
Third, you propose Yhat +- 1.96 * se(Yhat). I think this is the same as the previous formulation but with some confusion about the difference between standard deviation and standard error. The standard error should shrink toward zero as the sample size increases. However, we do not want such behavior. That's why prediction intervals, such as the first proposed formulation, do not shrink to have zero width.
Fourth, you propose to fit the Sklearn model with a quantile parameter set to .05 to generate upper and lower bounds. My suspicion is that this implements the first of the proposed formulations, complete with the caveats about non-Gaussian conditional distributions and non-constant error-term variance.
Another idea is to use explicit quantile modeling, that is, quantile regression. This has appeal, since quantile modeling is not bound to a particular distribution of the error term the way the above techniques are bound to a Gaussian error term. Sure, you could replace that $1.96$ with, say, the quantiles of a t-distribution if you are comfortable assuming t-distributed errors with a particular number of degrees of freedom, but even that is quite restrictive unless you really believe that is conditional distribution.
By using explicit quantile modeling, you let the data tell you what the conditional quantiles are, rather than making assumptions about particular distribution families that might be wildly inaccurate.
However, quantile modeling is not a silver bullet! Flexibility in the model is needed to allow for prediction intervals from the quantile models to vary in width. You might want wide prediction intervals in certain regions and narrowed intervals in others. For instance, in the (granted, contrived) example below, there would have to be flexibility in the modeling to capture the oscillating error variance.
set.seed(2024)
N <- 1000
x <- runif(N, -pi, pi)
Ey <- 1 - x
e <- rnorm(N, 0, sin(x)^2)
y <- Ey + e
plot(x, y)
"Fine, I'll be throwing a bunch of flexibility at the problem by using deep learning, anyway," you say? That flexibility can result in what I remember an arXiv paper referring to as the "embarrassing" problem of quantile crossover, where high predicted quantiles like $0.975$ are lower than low predicted quantiles like $0.025$. That is, you predict the lower limit of the prediction range to be $7$ while the upper limit of the prediction range to be $1$. Maybe you just flip them in that situation to say the prediction interval is $(1, 7)$, but then you are not using the explicit quantiles, and the prediction interval endpoints lack their usual interpretation and meaning. They might not wind up being very good prediction interval endpoints!
That arXiv paper made attempts to keep quantile crossover from happening in their quantile modeling. I recall these attempts not being entirely successful.
Thus, I call this an open problem that lacks a definite solution.
