Phase stability of stable and metastable vanadium nitrides is studied using density functional theory ͑DFT͒ based total-energy calculations combined with cluster expansion Monte Carlo simulation and supercell methods. We have computed the formation enthalpy of the various stable and metastable vanadium nitride phases considering the available structural models and found that the formation enthalpies of the different phases decrease in the same order as they appear in the experimental aging sequence. DFT calculations are known to show stoichiometric V 2 N to be polymorphic in ⑀-Fe 2 N and -Fe 2 N structures within a few meV and VN to be more stable in WC͑B h ͒ phase than in the experimentally observed NaCl͑B1͒ structure. As these nitrides are known to be generally nonstoichiometric due to presence of nitrogen vacancies, we used cluster expansion and supercell methods for examining the effect of nitrogen vacancies on the phase stability. It is found that nitrogen vacancies, represented by ᮀ, stabilize ⑀-Fe 2 N phase of V 2 N 1−x ᮀ x and NaCl͑B1͒ phase of VN 1−x ᮀ x compared to -Fe 2 N and WC͑B h ͒ phases respectively, rendering the computed phase stability scenario to be in agreement with experiments. Analysis of supercell calculated electronic density of states ͑DOS͒ of VN 1−x ᮀ x with varying x, shows that the nitrogen vacancies increase the DOS at Fermi level in WC phase, whereas they decrease the DOS in NaCl phase. And this serves as the mechanism of enhancement of the stability of the NaCl phase. Monte Carlo simulations were used for computing the finite temperature formation enthalpies of these phases as a function of nitrogen-vacancy concentration and found close agreement for NaCl͑B1͒ phase of VN 1−x ᮀ x for which measured values are available.