We present a complete theoretical treatment of Stark effects in bulk doped silicon, whose predictions are supported by experimental measurements. A multi-valley effective mass theory, dealing non-perturbatively with valley-orbit interactions induced by a donor-dependent central cell potential, allows us to obtain a very reliable picture of the donor wave function within a relatively simple framework. Variational optimization of the 1s donor binding energies calculated with a new trial wave function, in a pseudopotential with two fitting parameters, allows an accurate match of the experimentally determined donor energy levels, while the correct limiting behavior for the electronic density, both close to and far from each impurity nucleus, is captured by fitting the measured contact hyperfine coupling between the donor nuclear and electron spin.We go on to include an external uniform electric field in order to model Stark physics: With no extra ad hoc parameters, variational minimization of the complete donor ground energy allows a quantitative description of the field-induced reduction of electronic density at each impurity nucleus. Detailed comparisons with experimental values for the shifts of the contact hyperfine coupling reveal very close agreement for all the donors measured (P, As, Sb and Bi). Finally, we estimate field ionization thresholds for the donor ground states, thus setting upper limits to the gate manipulation times for single qubit operations in Kane-like architectures: the Si:Bi system is shown to allow for A gates as fast as ≈10 MHz.