By default it seems GPy doesn’t subtract the mean of the data prior to fitting.

So it’s worth including a mean function that does this:

 m = GPy.models.GPRegression(
       X,
       y,
       mean_function=GPy.mappings.Constant(
              input_dim=1,output_dim=1,value=np.mean(y))
     )

Then one needs to fix this value:

m.mean_function.C.fix()