In insulating lanthanides, unquenched orbital momentum and weak crystal-field (CF) splitting of the atomic $J$ multiplet at lanthanide ions result in a highly ranked (multipolar) exchange interaction between them and a complex low-temperature magnetic order not fully uncovered by experiment. Explicitly correlated ab initio methods proved to be highly efficient for an accurate description of CF multiplets and magnetism of individual lanthanide ions in such materials. Here we extend this ab initio methodology and develop a first-principles microscopic theory of multipolar exchange interaction between $J$ multiplets in $f$-metal compounds. The key point of the approach is a complete account of Goodenough's exchange mechanism along with traditional Anderson's superexchange and other contributions, the former being dominant in many lanthanide materials. Application of this methodology to the description of the ground-state order in the neodymium nitride with rocksalt structure reveals the multipolar nature of its ferromagnetic order. We found that the primary and secondary order parameters (of ${T}_{1u}$ and ${E}_{g}$ symmetry, respectively) contain non-negligible $J$-tensorial contributions up to the ninth order. The calculated spin-wave dispersion and magnetic and thermodynamic properties show that they cannot be simulated quantitatively by confining to the ground CF multiplet on the Nd sites. Our results demonstrate that the ab initio approach to the low-energy Hamiltonian represents a powerful tool for the study of materials with complex magnetic order.