  Fully coupled McKean-Vlasov forward-backward stochastic differential
equations (MV-FBSDEs) arise naturally from large population optimization
problems. Judging the quality of given numerical solutions for MV-FBSDEs, which
usually require Picard iterations and approximations of nested conditional
expectations, is typically difficult. This paper proposes an a posteriori error
estimator to quantify the $L^2$-approximation error of an arbitrarily generated
approximation on a time grid. We establish that the error estimator is
equivalent to the global approximation error between the given numerical
solution and the solution of a forward Euler discretized MV-FBSDE. A crucial
and challenging step in the analysis is the proof of stability of this Euler
approximation to the MV-FBSDE, which is of independent interest. We further
demonstrate that, for sufficiently fine time grids, the accuracy of numerical
solutions for solving the continuous MV-FBSDE can also be measured by the error
estimator. The error estimates justify the use of residual-based algorithms for
solving MV-FBSDEs. Numerical experiments for MV-FBSDEs arising from mean field
control and games confirm the effectiveness and practical applicability of the
error estimator.
