A smoothed aggregation multigrid method is presented for the numerical calculation of the stationary probability vector of irreducible Markov chains. It is shown how smoothing the interpolation and restriction operators can dramatically increase the efficiency of aggregation multigrid methods for Markov chains that have been proposed in the literature. The proposed smoothing approach is inspired by smoothed aggregation multigrid for linear systems, supplemented with a new lumping technique that assures well-posedness of the coarse-level problems: the coarse-level operators are singular M-matrices on all levels, resulting in strictly positive coarse-level corrections on all levels. Numerical results show how these methods lead to nearly-optimal multigrid efficiency for a representative set of test problems, both when geometric and algebraic aggregation strategies are used.