We present a simple derivation of a Feynman-Kac type formula to study fermionic systems. In this approach the real time or the imaginary time dynamics is expressed in terms of the evolution of a collection of Poisson processes. A computer implementation of this formula leads to a family of algorithms parametrized by the values of the jump rates of the Poisson processes. From these an optimal algorithm can be chosen which coincides with the Green Function Monte Carlo method in the limit when the latter becomes exact.A crucial issue in quantum Monte Carlo methods is the choice of the most convenient stochastic process to be used in the simulation of the dynamics of the system. This aspect is particularly evident in the case of fermion systems [1][2][3][4] due to the anticommutativity of the variables involved which for long time have not lent themselves to direct numerical evaluation. In a recent paper [5] progress has been made in this direction by providing exact probabilistic expressions for quantities involving variables belonging to Grassmann or Clifford algebras. In particular, the real time or the imaginary time evolution operator of a Fermi system or a Berezin integral can be expressed in terms of an associated stochastic dynamics of a collection of Poisson processes. This approach depends on an older general formula to represent probabilistically the solution of a system of ordinary differential equations (ODE) in terms of Poisson processes [8]. In this paper, we present a simple derivation of a similar formula to study fermionic systems, in particular, the Hubbard model. However, the fermionic nature of the Hamiltonian plays no special role and similar representations can be written down for any system described by a finite Hamiltonian matrix. A computer implementation of this formula leads to a family of algorithms parametrized by the values of the jump rates of the Poisson processes. For an optimal choice of these parameters we obtain an algorithm which coincides with the well known Green Function Monte Carlo method in the limit when the latter becomes exact [6]. In this way we provide a theoretical characterization of GFMC.Let us consider the Hubbard Hamiltonianwhere Λ ⊂ Z d is a finite d-dimensional lattice with cardinality |Λ|, {1, . . . , |Λ|} some total ordering of the lattice points, and c iσ the usual anticommuting destruction operators at site i and spin index σ. In this paper, we are interested in evaluating the matrix elements n ′ |e −Ht |n where n = (n 1↑ , n 1↓ , . . . , n |Λ|↑ , n |Λ|↓ ) are the occupation numbers taking the values 0 or 1 [7]. The total number of fermions per spin component is a conserved quantity, therefore we consider only configurations n and n ′ such thatIn the following we shall use the mod 2 addition n ⊕ n ′ = (n + n ′ ) mod 2. Let Γ = {(i, j), 1 ≤ i < j ≤ |Λ| : η ij = 0} and |Γ| its cardinality. For simplicity, we will assume that η ij = η if (i, j) ∈ Γ and γ i = γ. By introducingwhere 1 iσ = (0, . . . , 0, 1 iσ , 0, . . . , 0), andthe following representation holds...