Averaging magnitude squared coherence across multiple time series
I believe I might have a solution... Editing A Donda's work above, I think a simple rearrangement of the last two lines is necessary:

sxya = (abs(sxy1) + abs(sxy2)) / 2; 
ca = sxya.^ 2 ./ sxa ./ sya;

So here, the average of the magnitudes is taken, rather than the average of the complex CPSDs. Then, this average can be squared.