According to experimental observations, the temperature dependence of self-diffusion coefficient in most body-centered cubic metals (bcc) exhibits non-Arrhenius behavior. The origin of this behavior is likely related to anharmonic vibrational effects at elevated temperatures. However, it is still debated whether anharmonicity affects more the formation or migration of monovacancies, which are known to govern the self-diffusion. In this extensive atomistic simulation study we investigated thermodynamic properties of monovacancies in bcc molybdenum, here taken as a representative model system, from zero temperature to the melting point. We combined first-principles calculations and classical simulations based on three widely used interatomic potentials for Mo. In our analysis we employ static and dynamic atomistic calculations as well as statistical sampling techniques and thermodynamic integration to achieve thorough information about temperature variations of vacancy formation and migration free energies and diffusivities. In addition, we carry out large-scale molecular dynamics simulations that enable direct observation of high-temperature self-diffusion at the atomic scale. By scrutinizing the results obtained by different models and methods, we conclude that the peculiar self-diffusion behavior is likely caused by strong temperature dependence of the vacancy formation energy.