applying a quantile loss function to optimize an exponential smoothing model for probabilistic forecasting

Question

I am trying to build a model that generate prediction intervals from a traditional winters-holt time serie. There is a built-in function that does it in statmodels but there are a few missing features preventing me from fully depend on it. I am looking for a critic eye and don't know anywhere else to post:

What I did is use a quantile loss pinball function to optimize the parameters (level, trend and seasonality) of the model and generate different forecast corresponding to different quantile values that I input in the quantile loss function.

first of all: is it good practice to optimize the parameters for each quantile

second, is the described process viable to get prediction intervals?

Now, in the results dataset, I observe none difference between predictions made with a quantile of 0.05 and a 0.95 or worst, there are even timeseries where the results give lower values for the 0.95 quantile value.

I am confused and wondered what could be wrong in the method knowing that the model works just fine with finding regular point series using mean square error loss function

1. here is my quantile loss function:

def quantile_loss(q,y_p, y):
    a = np.where((y > y_p), q *(y-y_p), (y_p - y)*(1-q))

    return a

and here is the code for the parameter optimization:

def HoltWinterLowHightimeseriesCVscore(params,quantile_values, data, loss_function=quantile_loss,slen=12):
    """
        Returns error on CV

        params - vector of parameters for optimization
        series - dataset with timeseries
        slen - season length for Holt-Winters model
    """
    # errors array
    errors = []

    values = data
    alpha, beta, gamma = params

    # set the number of folds for cross-validation
    tscv = TimeSeriesSplit(n_splits=3)

    # iterating over folds, train model on each, forecast and calculate error
    for train, test in tscv.split(values):
        model = HoltWintersLowHigh(series=values, slen=slen,
                            alpha=alpha, beta=beta, gamma=gamma, n_preds=12)
        model.triple_exponential_smoothing()
        predictions = model.result[-len(test):]
        actual = values[test]
        error = loss_function(quantile, predictions, actual)
        errors.append(error)

    return np.mean(np.array(errors))

and finally here is the final part where the functions are called to make the predictions:

forecast = {}

for i in seasonal_profile_df.index:
    quantile_values = [0.92]

    if seasonal_profile_df['trend'].loc[i] == 'trending' and seasonal_profile_df['seasonality'].loc[i] == 'seasonal' and seasonal_profile_df['demand_level'].loc[i] == 'low' or seasonal_profile_df['variability'].loc[i] == 'high':

        index = pd.DatetimeIndex(new_df.index)
        series = pd.Series(data=new_df.iloc[:, i], index=index)
        print(len(series))
        data = series[:-10]  #leave some data for testing
        x = [0, 0, 0]
        #for i in quantile:

        for j in quantile_values:
            quantile_values = j
    # Minimizing the loss function
            opt = minimize(HoltWinterLowHightimeseriesCVscore, x0=x,
                       args=( quantile_values,data, quantile_loss,),
                       method="TNC", bounds=((0, 1), (0, 1), (0, 1))
                       )
            alpha_final, beta_final, gamma_final = opt.x
            print(opt.x)
            model = HoltWintersLowHigh(series, slen=12,
                                alpha=alpha_final,
                                beta=beta_final,
                                gamma=gamma_final,
                                n_preds=12, scaling_factor=1.96)
            model.triple_exponential_smoothing()
            plotHoltWintersLowHigh(series, quantile_values)

            result= {"Id": seasonal_profile_df['Id'].loc[i]}
            result['results'] = model.result[-12:]

I am really hoping to get a fresh or more experienced eye on this or what is going wrong