sklearn GP return std dev is zero for predictions where it must be large

Question

I am trying regression using Gaussian processes sklearn package. The standard deviation on predictions are zero, where it must be larger.

kernel = ConstantKernel() + 1.0 * DotProduct() ** 0.3 + 1.0 * WhiteKernel()

gpr = GaussianProcessRegressor(
    kernel=kernel, 
    alpha=0.3, 
    normalize_y=True, 
    random_state=123, 
    n_restarts_optimizer=0
)

gpr.fit(X_train, y_train)

Here I have shown the samples from posterior after training the model. It clearly shows the standard deviation increases along with x-axis.

This is the output I got. As the value increases along x-axis the stddev must increase, where as it is showing zero stddev.

Acutal results should look something like this.

Is it a bug ?

Full Code to reproduce the issue.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, WhiteKernel, DotProduct

df = pd.read_csv('train.csv')
X_train = df[:,0].to_numpy().reshape(-1,1)
y_train = df[:,1].to_numpy()

X_pred = np.linspace(0.01, 8.5, 1000).reshape(-1,1)

# Instantiate a Gaussian Process model
kernel = ConstantKernel() + 1.0 * DotProduct() ** 0.3 + 1.0 * WhiteKernel()

gpr = GaussianProcessRegressor(
    kernel=kernel, 
    alpha=0.3, 
    normalize_y=True, 
    random_state=123, 
    n_restarts_optimizer=0
)

gpr.fit(X_train, y_train)

print(
    f"Kernel parameters before fit:\n{kernel} \n"
    f"Kernel parameters after fit: \n{gpr.kernel_} \n"
    f"Log-likelihood: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.3f} \n"
    f"Score = {gpr.score(X_train,y_train)}"
)

n_samples = 10
y_samples = gpr.sample_y(X_pred, n_samples)

for idx, single_prior in enumerate(y_samples.T):
    plt.plot(
            X_pred,
            single_prior,
            linestyle="--",
            alpha=0.7,
            label=f"Sampled function #{idx + 1}",
    )
plt.title('Sample from posterior distribution')
plt.show()

y_pred, sigma = gpr.predict(X_pred, return_std=True)

plt.figure(figsize=(10,6))
plt.plot(X_train, y_train, 'r.', markersize=3, label='Observations')
plt.plot(X_pred, y_pred, 'b-', label='Prediction',)
plt.fill_between(X_pred[:,0], y_pred-1*sigma, y_pred+1*sigma,
         alpha=.4, fc='b', ec='None', label='68% confidence interval')
plt.fill_between(X_pred[:,0], y_pred-2*sigma, y_pred+2*sigma,
         alpha=.3, fc='b', ec='None', label='95% confidence interval')
plt.fill_between(X_pred[:,0], y_pred-3*sigma, y_pred+3*sigma,
         alpha=.1, fc='b', ec='None', label='99% confidence interval')
plt.legend()
plt.show()