A novel algorithmic method is proposed to fit matrix geometric distributions of desired order to empirical data or arbitrary discrete distributions. The proposed method effectively combines two existing approaches from two different disciplines: well-established model reduction methods used in system theory and moment matching methods of applied probability that employ second-order discrete phase-type distributions. The proposed approach is validated with exhaustive numerical examples including well-known statistical data. Copyright