We extend the recently developed theory of bulk orbital magnetization to finite electric fields, and use it to calculate the orbital magnetoelectric (ME) response of periodic insulators. Working in the independent-particle framework, we find that the finite-field orbital magnetization can be written as a sum of three gauge-invariant contributions, one of which has no counterpart at zero field. The extra contribution is collinear with and explicitly dependent on the electric field. The expression for the orbital magnetization is suitable for first-principles implementations, allowing one to calculate the ME response coefficients by numerical differentiation. Alternatively, perturbation-theory techniques may be used, and for that purpose we derive an expression directly for the linear ME tensor by taking the first field-derivative analytically. Two types of terms are obtained. One, the 'Chern–Simons' term, depends only on the unperturbed occupied orbitals and is purely isotropic. The other, 'Kubo' terms, involve the first-order change in the orbitals and give isotropic as well as anisotropic contributions to the response. In ordinary ME insulators all terms are generally present, while in strong Z₂ topological insulators only the Chern–Simons term is allowed, and is quantized. In order to validate the theory, we have calculated under periodic boundary conditions the linear ME susceptibility for a 3D tight-binding model of an ordinary ME insulator, using both the finite-field and perturbation-theory expressions. The results are in excellent agreement with calculations on bounded samples.