This work is aimed at developing a density functional theory (DFT) approach that can be used to calculate 13C NMR shifts in any organometallic diamagnetic Pd complexes. Comparative analysis of calculated (GIAO method, DFT level) and experimental 13C NMR shifts for a wide range of diamagnetic palladium complexes (62 complexes in total) showed that the theory reproduces the experimental data well. A number of different basis sets, as well as quasi-relativistic and full-relativistic approximations, were tested. On the whole, the chemical shifts of carbon atoms directly bonded to Pd can be calculated within the framework of the Kohn–Sham theory level for most complexes with classical coordination bonds. The exceptions are complexes with carbons covalently bonded to metal and some NHC carbons due to relativistic effects. To summarize, in practice, the PBE0/{6-311G(2d,2p); Pd(SDD)}//PBE0/{6-31+G(d); Pd(SDD)} approximation can be recommended as a first step for most cases. Then, for complexes with the NHCs and covalently bonded ligands, 13C shifts should be calculated at a fully relativistic matrix Dirac–Kohn–Sham (mDKS) level, at least for atoms directly bonded to Pd (RMSE = 5.0 ppm). In all cases, a linear scaling procedure is necessary to minimize systematic errors.