GP with sigmoid link function for time-series data

So no particular reason jumps out to me why this model won’t work, but with a model like this it’s going to be super touchy. Have you tried using fixed values for all the GP lengthscales? Also, don’t forget to scale the GPs overall. You’ll get very different GPs with

cov1 = pm.gp.cov.ExpQuad(1, lengthscale)
# versus (with scale parameter eta)
eta = 100
cov2 = eta**2 * pm.gp.cov.ExpQuad(1, lengthscale)

I would try to start by get something that gives good prior predictive samples with fixed lengthscales and GP scales (eta above), and then gradually try and incorporate more uncertainty in those parameters.