We have compared the performance of density functional theory (DFT) using five different exchange-correlation functionals with four coupled cluster theory based wave function methods in the calculation of geometrical derivatives of the polarizability tensor of methane. The polarizability gradients of hydrocarbons are important ingredients in the simulation of their electron energy loss spectra and reliable but cost-effective methods for obtaining the gradients need to be found. In the present work we present results of a systematic investigation on methane as a prototype molecule with special focus on DFT methods. The KT3, B3LYP, CAM-B3LYP, B97-2 and PBE0 DFT exchange-correlation functionals and the highly correlated wave function methods SOPPA(CCSD), CCSD-LR, CCSD and CCSD(T) were employed in combination with a series of eleven basis sets. Comparison of the DFT results with CCSD(T)/daug-cc-pVQZ reference values reveals that none of the investigated DFT approaches reaches the accuracy of correlated wave function based methods and that the best DFT results are obtained with the PBE0 exchange-correlation functional and Sadlej's polarized valence triple zeta basis set. The SOPPA(CCSD) method, on the other hand, produces results in close agreement with the more expensive pure coupled cluster methods.