Complex reaction networks commonly appear in natural systems. Some examples are chemical, ecological, metabolic and gene-expression networks.These networks can be described by graphs, in which the nodesrepresent reactive species and the edges represent reactions.Computer simulations of these networks are commonly done using rate equations,which provide the time dependent concentrations of thereactive species as well as the reaction rates.These equations are easy to construct and efficient to integrate numerically. However, they are applicable only whenthe populations of reactive species are large andfluctuations are negligible.In case that the reactive species appear in low copy numbers,stochastic methods, based on the master equation, are needed.However, the master equation is not feasible for complex networks,because the number of equations proliferatesexponentially with the number of reactive species.Here we present a new computational method, based on moment equations, which dramatically reduces the number of equations.It enables to efficiently simulate fluctuating reactionnetworks of any level of complexity.The method requires only one equation for each reactive speciesand one equation for each reaction,thus reducing the number of equtaions to the absolute minimum for a stochastic simulation. The reduction is achieved with no compromise in the accuracy of the results.