This paper concerns the application of Grassmann phase space theory (GSPT) to treat the dynamical evolution of systems of identical fermions, such as ultracold gases of fermionic atoms. Phase space theory (which originated from quantum optics) is increasing in importance since it overcomes certain issues associated with other theoretical methods, such as Greens functions, variational methods, quantum-Monte-Carlo equations, etc. In phase-space theory quantum states are represented by quasi-probability distribution functions of phase space variables associated with canonical system operators—such as annihilation, creation operators. Evolution is described via a Fokker-Planck equation for the distribution function, which is equivalent to Ito stochastic equations for (time dependent) stochastic phase space variables. Quantum correlation functions given as averages of products of phase space variables over the quasi-probability distributions then become stochastic averages of products of stochastic phase space variables. In GSPT, the phase space variables are Grassmann numbers, but as computer representation of g-numbers is difficult, Grassmann phase space methods were regarded as being computationally inaccessible. However, previous work using the un-normalised B distribution shows that computer representation of Grassmann variables is unnecessary. Stochastic averages of products for quantum correlation functions at later times are related linearly to stochastic averages at earlier times via stochastic matrices only involving c-numbers. Thus, GSPT calculations of quantum correlation functions now only involve c-number computations. This paper presents the first correct numerical calculation of a quantum correlation function for a fermionic system using stochastic methods based on Grassmann phase space theory, namely the time dependence of the coherence between two Cooper pair states in a four-mode fermion system, where the short and finite time solutions can be compared to known exact results. Good agreement between the stochastic and exact results is found, showing that GPST is a valid approach for treating fermionic systems. The treatment of time evolution involves a novel use of the eigenvalues and biorthogonal column eigenvectors of a stochastically determined c-number matrix M and its transpose. Other topics of interest in ultra-cold fermi gases for which the GSPT could be applied are highlighted, such as the strong interaction regime for the BEC/BCS crossover achieved using magnetically tuned Feshbach resonance techniques.