The development of an advanced Reduced Order Model (ROM) including four heated channels and meant to study global and regional Boiling Water Reactor (BWR) instabilities is described. The ROM contains three sub-models: a neutron-kinetic model (describing neutron transport), a thermal-hydraulic model (describing fluid transport) and a heat transfer model (describing heat transfer between the fuel and the coolant). All these three models are coupled to each other using two feedback mechanisms: the void feedback and the doppler feedback mechanisms. Each of the sub-models is described by a set of reduced ordinary differential equations, derived from the corresponding time- and space-dependent partial differential equations, by using different types of approximations and mathematical techniques that are explained in this paper. One of the novelties of the present ROM is that it takes the effect of the first three neutronic modes into account, namely the fundamental, first, and second azimuthal modes. In order to have a proper representation of both azimuthal modes and of their dependence on the thermal-hydraulic conditions in the heated channels, a four heated channel ROM was constructed. Another novelty of the present work is to develop a special methodology which guarantees the full consistency between the spatial discretization procedures used in the dynamical calculations and the ones implemented in the static case. Accordingly, a re-computation of the static solution based on the CORE SIM tool was embedded into the ROM in such a way that the balance equations expressing the conservation of neutron balance, heat generation, and mass, momentum, enthalpy for the flow, could be fulfilled for the steady-state solution of the coupled neutron-kinetic/thermal-hydraulic problem. Once the static problem is solved, the time-dependent solution in case of a perturbed system can be determined. Moreover, a non-uniform power profile representing different heat production rates in the one- and two-phase regions was introduced into the ROM. Careful attention was paid to the determination of the coupling coefficients for the reactivity effects related to both void fraction and fuel temperature, so that such coefficients correspond to the re-computed static solution. The evaluation of these coefficients was based on the cross-section perturbations estimated by the SIMULATE-3 code, and on the different neutronic eigenmodes of the heterogeneous core determined by the CORE SIM tool. The time-evolution of the above-mentioned modes was generated by the ROM for various operational conditions in order to demonstrate the feasibility and capability of the developed model.