Statsmodels ARIMA: how to get confidence/prediction interval?

In general, the forecast and predict methods only produce point predictions, while the get_forecast and get_prediction methods produce full results including prediction intervals.

In your example, you can do:

forecast = model.get_forecast(123)
yhat = forecast.predicted_mean
yhat_conf_int = forecast.conf_int(alpha=0.05)

If your data is a Pandas Series, then yhat_conf_int will be a DataFrame with two columns, lower <name> and upper <name>, where <name> is the name of the Pandas Series.

If your data is a numpy array (or Python list), then yhat_conf_int will be an (n_forecasts, 2) array, where the first column is the lower part of the interval and the second column is the upper part.

To generate prediction intervals as opposed to confidence intervals (which you have neatly made the distinction between, and is also presented in Hyndman's blog post on the difference between prediction intervals and confidence intervals), then you can follow the guidance available in this answer.

You could also try to compute bootstrapped prediction intervals, which is laid out in this answer.

Below, is my attempt at implementing this (I'll update it when I get the chance to check it in more detail):

def bootstrap_prediction_interval(y_train: Union[list, pd.Series],
                                  y_fit: Union[list, pd.Series],
                                  y_pred_value: float,
                                  alpha: float = 0.05,
                                  nbootstrap: int = None,
                                  seed: int = None):
    """
    Bootstraps a prediction interval around an ARIMA model's predictions.
    Method presented clearly here:
        - https://stats.stackexchange.com/a/254321
    Also found through here, though less clearly:
        - https://otexts.com/fpp3/prediction-intervals.html
    Can consider this to be a time-series version of the following generalisation:
        - https://saattrupdan.github.io/2020-03-01-bootstrap-prediction/

    :param y_train: List or Series of training univariate time-series data.
    :param y_fit: List or Series of model fitted univariate time-series data.
    :param y_pred_value: Float of the model predicted univariate time-series you want to compute P.I. for.
    :param alpha: float = 0.05, the prediction uncertainty.
    :param nbootstrap: integer = 1000, the number of bootstrap sampling of the residual forecast error.
        Rules of thumb provided here:
            - https://stats.stackexchange.com/questions/86040/rule-of-thumb-for-number-of-bootstrap-samples
    :param seed: Integer to specify if you want deterministic sampling.

    :return: A list [`lower`, `pred`, `upper`] with `pred` being the prediction
    of the model and `lower` and `upper` constituting the lower- and upper
    bounds for the prediction interval around `pred`, respectively.
    """

    # get number of samples
    n = len(y_train)

    # compute the forecast errors/resid
    fe = y_train - y_fit

    # get percentile bounds
    percentile_lower = (alpha * 100) / 2
    percentile_higher = 100 - percentile_lower

    if nbootstrap is None:
        nbootstrap = np.sqrt(n).astype(int)
    if seed is None:
        rng = np.random.default_rng()
    else:
        rng = np.random.default_rng(seed)

    # bootstrap sample from errors
    error_bootstrap = []
    for _ in range(nbootstrap):
        idx = rng.integers(low=n)
        error_bootstrap.append(fe[idx])

    # get lower and higher percentiles of sampled forecast errors
    fe_lower = np.percentile(a=error_bootstrap, q=percentile_lower)
    fe_higher = np.percentile(a=error_bootstrap, q=percentile_higher)

    # compute P.I.
    pi = [y_pred_value + fe_lower, y_pred_value, y_pred_value + fe_higher]

    return pi