Abstract-In the stochastic formulation of chemical reactions, the dynamics of the the first M -order moments of the species populations generally do not form a closed system of differential equations, in the sense that the time-derivatives of first M -order moments generally depend on moments of order higher thanM . However, for analysis purposes, these dynamics are often made to be closed by approximating the needed derivatives of the first M -order moments by nonlinear functions of the same moments. These functions are called the moment closure functions.Recent results have introduced the technique of derivativematching, where the moment closure functions are obtained by first assuming that they exhibit a certain separable form, and then matching time derivatives of the exact (not closed) moment equations with that of the approximate (closed) equations for some initial time and set of initial conditions. However, for multi-species reactions these results have been restricted to second order truncations, i.e, M = 2. This paper extends these results by providing explicit formulas to construct moment closure functions for any arbitrary order of truncation M . We show that with increasing M the closed moment equations provide more accurate approximations to the exact moment equations. Striking features of these moment closure functions are that they are independent of the reaction parameters (reaction rates and stoichiometry) and moreover the dependence of higher-order moment on lower order ones is consistent with the population being jointly lognormally distributed. To illustrate the applicability of our results we consider a simple bi-molecular reaction. Moment estimates from a third order truncation are compared with estimates obtained from a large number of Monte Carlo simulations.