By using the variational Monte Carlo method, we study the magnetization process of the Kitaev honeycomb model in a magnetic field. Our trial wavefunction is a generalized Bardeen-Cooper-Schrieffer wave function with the Jastrow correlation factor, which exactly describes the ground state of the Kitaev model at zero magnetic field using the Jordan-Wigner (JW) transformation. We find that two phase transitions occur for the antiferromagnetic Kitaev coupling, while only one phase transition occurs for the ferromagnetic Kitaev coupling. For the antiferromagnetic Kitaev coupling, we also find that the topology of the momentum distribution of the JW fermions changes at the transition point from the Kitaev spin liquid to an intermediate state. Our numerical results indicate that the intermediate state between the Kitaev spin liquid and the fully polarized phases stably exists in the bulk system on two dimensions for the antiferromagnetic Kitaev coupling against many-body correlations.