We extend the framework of Neuts' matrix analytic approach to a reflected process generated by a discrete time multidimensional Markov additive process. This Markov additive process has a general background state space and a real vector valued additive component, and generates a multidimensional reflected process. Our major interest is to derive a closed form formula for the stationary distribution of this reflected process. To this end, we introduce a real valued level, and derive new versions of the Wiener-Hopf factorization for the Markov additive process with the multidimensional additive component. In particular, it is represented by moment generating functions, and we consider the domain for it to be valid. Our framework is general enough to include multi-server queues and/or queueing networks as well as non-linear time series which are currently popular in financial and actuarial mathematics. Our results yield structural results for such models. As an illustration, we apply our results to extend existing results on the tail behavior of reflected processes. A major theme of this work is to connect recent work on matrix analytic methods to classical probabilistic studies on Markov additive processes. Indeed, using purely probabilistic methods such as censoring, duality, level crossing and time-reversal (which are known in the matrix analytic methods community but date back to Arjas & Speed (1973) and Pitman (1972), we extend and unify existing results in both areas.