We developed a new formulation and derived a set of equations for two-time correlation functions in discrete state space for general Markov chains. Numerical solutions of these equations lead to exact results of correlation functions, which is more efficient and accurate for complex systems than normal Gillespie simulations. The formalism is applied to two simple systems as examples: a two state system and self-regulators with results compared with Monte Carlo simulations and mean-field theory. For self-regulators, from non-adiabatic to adiabatic regimes, we observed both monotonic and turnover kinetic behavior of response time, which can be can tested by experiments. (C) 2010 Elsevier B.V. All rights reserved.