Few analytical methods exist for quantitative studies of large fluctuations in stochastic systems. In this article, we develop a simple diagrammatic approach to the chemical master equation that allows us to calculate multi-time correlation functions which are accurate to any desired order in van Kampenʼs system size expansion. Specifically, we present a set of Feynman rules from which this diagrammatic perturbation expansion can be constructed algorithmically. We then apply the methodology to derive in closed form the leading order corrections to the linear noise approximation of the intrinsic noise power spectrum for general biochemical reaction networks. Finally, we illustrate our results by describing noise-induced oscillations in the Brusselator reaction scheme which are not captured by the common linear noise approximation.