We study the numerical behavior of stationary one-step or two-step matrix splitting iteration methods for solving large sparse systems of linear equations. We show that inexact solutions of inner linear systems associated with the matrix splittings may considerably influence the accuracy of the approximate solutions computed in finite precision arithmetic. For a general stationary matrix splitting iteration method, we analyze two mathematically equivalent implementations and discuss the conditions when they are componentwise or normwise forward or backward stable. We distinguish two different forms of matrix splitting iteration methods and prove that one of them is significantly more accurate than the other when employing inexact inner solves. The theoretical results are illustrated by numerical experiments with an inexact one-step and an inexact two-step splitting iteration method.