Abstract. We consider a family of matrix valued orthogonal polynomials obtained by Pacharoni and Tirao in connection with spherical functions for the pair (SU(N + 1), U(N )); see [I. Pacharoni and J. A. Tirao, Constr. Approx., 25 (2007), pp. 177-192]. After an appropriate conjugation, we obtain a new family of matrix valued orthogonal polynomials where the corresponding block Jacobi matrix is stochastic and has special probabilistic properties. This gives a highly nontrivial example of a nonhomogeneous quasi-birth-and-death process for which we can explicitly compute its "nstep transition probability matrix" and its invariant distribution. The richness of the mathematical structures involved here allows us to give these explicit results for a several parameter family of quasi-birth-and-death processes with an arbitrary (finite) number of phases. Some of these results are plotted to show the effect that choices of the parameter values have on the invariant distribution.