We present a numerical simulation of the formation of a massive star using MonteCarlo-based radiation hydrodynamics (RHD). The star forms via stochastic disc accretion and produces fast, radiation-driven bipolar cavities. We find that the evolution of the infall rate (considered to be the mass flux across a 1500 au spherical boundary), and the accretion rate onto the protostar, are broadly consistent with observational constraints. After 35 kyr the star has a mass of 25 M and is surrounded by a disc of mass 7 M and 1500 au radius, and we find that the velocity field of the disc is close to Keplerian. Once again these results are consistent with those from recent high-resolution studies of discs around forming massive stars. Synthetic imaging of the RHD model shows good agreement with observations in the near-and far-IR, but may be in conflict with observations that suggests that MYSOs are typically circularly symmetric on the sky at 24.5 µm. Molecular line simulations of a CH 3 CN transition compare well with observations in terms of surface brightness and line width, and indicate that it should be possible to reliably extract the protostellar mass from such observations.