Hi @essicolo, it seems that the current pymc.gp.cov.Coregion does not support multi-outputs yet. The current implementation does not have the output_dim parameter compared to GPy Coregionalize
As the multitasks GP seems not fully supported in pymc yet, you may also try other libraries like Gpy, GPflow (tensorflow backend), GPytorch (pytorch backend), or recently mogptk (pytorch backend).
Here is my simple implementation using GPytorch with your data set. Although, I don’t like it as we need to write more boilerplate codes.