We apply the theory of the nuclear magnetic resonance (NMR) chemical shift for paramagnetic systems to demanding cobalt(II) complexes. Paramagnetic NMR (pNMR) chemical shift results by density-functional theory (DFT) can be very far from the experimental values. Therefore, it is of interest to investigate the applicability of electron-correlated ab initio computational methods to achieve useful accuracy. Here, we use ab initio wave function based electronic structure methods to calculate the pNMR chemical shift within the theoretical framework established recently. We applied the N-electron valence-state perturbation theory (NEVPT2) on three Co(II) systems, where the active space of the underlying complete active space self-consistent field (CASSCF) wave function consists of seven electrons in the five metal 3d orbitals. These complexes have the S = 3/2 electronic ground state consisting of two doublets separated by zero-field splitting (ZFS). To calculate the hyperfine coupling tensor A, DFT was used, while the g- and ZFS-tensors were calculated using the ab initio CASSCF and NEVPT2 methods. These results were combined to obtain the total chemical shifts. The shifts obtained from these calculations are in generally good agreement with the experimental results, in some cases suggesting a reassignment of the signals. The accuracy of this mixed ab initio/DFT approach is very promising for further applications to demanding pNMR problems involving transition metals.