There has been considerable scientific interest in third-order nonlinear optical materials for photonic applications. In particular, materials exhibiting a strong electronic optical Kerr effect serve as essential components in the ultrafast nonlinear photonic devices and are instrumental in the development of all-optical signal processing technologies. Therefore, the accurate prediction of material-relevant properties, such as second hyperpolarizabilities, remains a key topic in the search for efficient photonic materials. However, the field standards in quantum chemical computation are still inconsistent, as studies often lack a firm statistical foundation. This work presents a comprehensive in silico investigation based on multiple full-factorial experiments, aiming to clarify the strengths and limitations of various computational approaches. Our results indicate that the coupled-cluster approach at the CCSD level in its current response-equation implementations is not yet able to outperform the range-separated hybrid density functionals, such as LC-BLYP(0.33). The exceptional performance of the specifically tailored basis set Sadlej-pVTZ is also described. Not only was the presence of diffuse functions found to be mandatory, but also adding ample polarization functions is shown to be inefficient resource-wise. HF/Sadlej-pVTZ is proven to be reliable enough to use in molecular screening. Meta functionals are confirmed to produce poorly consistent results, and specific guidelines for constructing range-separated functionals for polarizability calculations are drawn out. Additionally, it was shown that many of the contemporary solvation models exhibit significant limitations in accurately capturing nonlinear optical properties. Therefore, further refinement in the current methods is pending. This extends to the statistical description as well: the mean absolute deviation descriptor is found to be deficient in rating various computational methods and should rather be replaced with the parameters of the linear correlation (the slope, the intercept, and the R2).