  We discuss the numerical implementation of two related representations of
fermionic density matrices which have been introduced in Annals of Physics 370,
12 (2016). In both of them, the density matrix is expanded in a basis of
Bargmann coherent states with weights given by the two phase space
distributions. We derive the equations of motion for the distributions when
imaginary time evolution is generated by the Hubbard Hamiltonian. One of them
is a Grassmann Fokker-Planck equation that can be re-cast into a remarkably
simple It\^{o} form involving solely complex variables. In spite of this simple
form, we demonstrate that complications arise in numerically computing the
expectation value of any observable. These are due to exponential growth in the
matrix elements of the stochastic propagator, delicate numerical sensitivity in
performing primitive linear algebra operations, and the re-appearance of a sign
problem.
