We present a theory for the temperature-dependent nuclear magnetic shielding tensor of molecules with arbitrary electronic structure. The theory is a generalization of Ramsey's theory for closed-shell molecules. The shielding tensor is defined as a second derivative of the Helmholtz free energy of the electron system in equilibrium with the applied magnetic field and the nuclear magnetic moments. This derivative is analytically evaluated and expressed as a sum over states formula. Special consideration is given to a system with an isolated degenerate ground state for which the size of the degeneracy and the composition of the wave functions are arbitrary. In this case, the paramagnetic part of the shielding tensor is expressed in terms of the g and A tensors of the electron paramagnetic resonance spin Hamiltonian of the degenerate state. As an illustration of the proposed theory, we provide an explicit formula for the paramagnetic shift of the central lanthanide ion in endofullerenes Ln@C(60), with Ln = Ce(3+), Nd(3+), Sm(3+), Dy(3+), Er(3+), and Yb(3+), where the ground state can be a strongly spin-orbit coupled icosahedral sextet for which the paramagnetic shift cannot be described by previous theories.