The interpretation of nuclear magnetic resonance (NMR) spectra of paramagnetic molecules is complex for experimentalists and theoreticians alike. The magnetic interaction between the unpaired electron(s) and the nucleus of interest can cause profound and unpredictable changes to the observed spectra when compared with an analogous diamagnetic system. Many of the lanthanides and actinides in their most stable oxidation states are characterised as paramagnetic, meaning the assignment of NMR chemical shifts in f-element systems is challenging. In this work the computational and theoretical methodology required to reliably calculate 1H, 13C and 29Si NMR shifts in paramagnetic f-element compounds (PNMR) has been explored. By arriving at a computational protocol we aim to aid the interpretation of experimental data, provide understanding of the factors that govern the observed spectra and develop improved computational approaches for handling the f-elements. Relativistic effects are an important consideration where the f-elements are concerned and two relativistic approximations have been implemented in this work for assessment. These are the exact-2-component (X2C) and atomic zeroth order regular approximation (aZORA) schemes. aZORA proved an attractive method and was employed to perform geometry optimisations, using density functional theory (DFT), on all the molecules studied. The calculation of PNMR chemical shifts requires spin hamiltonian parameters: the electronic g tensor, the hyperfine coupling tensor A and the NMR shielding tensor Ï. The computation of each of these terms is not routine for the f-elements and therefore calibrations against detailed experimental data were performed in order to establish a suitable computational scheme for each. DFT linear response theory was shown to calculate A and Ï. However g is highly dependent on a reliable description of low-lying excited states, and therefore state averaged complete active space self-consistent field (SA-CASSCF) calculations were preferred for this term. The computational protocol determined was applied to a selection of lanthanide and actinide complexes with f1 and f3 electronic configurations. When compared with experimental data the accuracy of the calculations were found to be sensitive to multiple factors including the relativistic scheme employed, the level of state averaging included in the g tensor calculation, the geometry used and the consideration of experimental conditions. Zero field splitting effects were not incorporated into this work and this omission proved crucial for the calculation of PNMR chemical shifts in f3 systems. However, useful and reliable accuracy was obtained for the f1 systems giving confidence that this protocol can be used to aid the interpretation of PNMR spectra of similar complexes.