We describe the computational ingredients for an approach to treat interacting fermion systems in the presence of pairing fields, based on path-integrals in the space of Hartree-Fock-Bogoliubov (HFB) wave functions. The path-integrals can be evaluated by Monte Carlo, via random walks of HFB wave functions whose orbitals evolve stochastically. The approach combines the advantage of HFB theory in paired fermion systems and many-body quantum Monte Carlo (QMC) techniques. The properties of HFB states, written in the form of either product states or Thouless states, are discussed. The states preserve forms when propagated by generalized one-body operators. They can be stabilized for numerical iteration. Overlaps and one-body Green's functions between two such states can be computed. A constrained-path or phaseless approximation can be applied to the random walks of the HFB states if a sign problem or phase problem is present. The method is illustrated with an exact numerical projection in the Kitaev model, and in the Hubbard model with attractive interaction under an external pairing field.