I've been playing around with this a bit today and I looked at your results. I get an R matrix of
1.041734 -0.00457 -0.01165
0.012072 0.992504 -0.00335
-0.00652 0.007879 0.979193
But (and it's a fairly big but) this matrix shouldn't be used as a simple adjustment matrix in the current HCFR. There is some renormalization that needs to be done (the sentence below equation 10 in the paper)
Say you start with meaured xyY, save the Y for later work out xyz then multiply by the above matrix to give abc then work out the adjusted xyY as x = a/(a+b+c) y = b/(a+b+c) Y=Y from before. Without the last bit that matrix would make things worse, say using your red example of measured 0.633, 0.334, using the matrix gives abc of (.6575, 0.3390, 0.0308) which then gives xy of 0.64, 0.33 as expected, but if you just applied that matrix to the XYZ values the error goes up.
Obviously I'm looking at fixing this in HCFR so that it all works properly.