The fine-scale stochastic behavior of genetic regulatory networks is often modeled using stochastic master equations. The inherently high computational complexity of the stochastic master equation simulation presents a challenge in its application to biological system modeling even when the model parameters can be properly estimated. In this article, we present a new approach to stochastic model simulation based on Kronecker product analysis and approximation of Zassenhaus formula for matrix exponentials. Simulation results on model biological systems illustrate the comparative performance of our modeling approach to stochastic master equations with significantly lower computational complexity. We also provide a stochastic upper bound on the deviation of the steady state distribution of our model from the steady state distribution of the stochastic master equation.