A data-driven algorithm recently proposed to solve the problem of model reduction by moment matching is extended to multi-input, multi-output systems. The algorithm is exploited for the model reduction of large-scale interconnected power systems and it offers, simultaneously, a low computational complexity approximation of the moments and the possibility to easily enforce constraints on the reduced order model. This advantage is used to preserve selected slow and poorly damped modes. The preservation of these modes has been shown to be important from a physical point of view and in obtaining an overall good approximation. The problem of the choice of the socalled tangential directions is also analyzed. The algorithm and the resulting reduced order model are validated with the study of the dynamic response of the NETS-NYPS benchmark system (68-Bus, 16-Machine, 5-Area) to multiple fault scenarios.Index Terms-Dynamic equivalents, model reduction of power systems, coherency, model reduction from data.