  We present a quantum Monte-Carlo algorithm for computing the perturbative
expansion in power of the coupling constant $U$ of the out-of-equilibrium
Green's functions of interacting Hamiltonians of fermions. The algorithm
extends the one presented in Phys. Rev. B 91 245154 (2015), and inherits its
main property: it can reach the infinite time (steady state) limit since the
computational cost to compute order $U^n$ is uniform versus time; the computing
time increases as $2^n$. The algorithm is based on the Schwinger-Keldysh
formalism and can be used for both equilibrium and out-of-equilibrium
calculations. It is stable at both small and long real times including in the
stationary regime, because of its automatic cancellation of the disconnected
Feynman diagrams. We apply this technique to the Anderson quantum impurity
model in the quantum dot geometry to obtain the Green's function and
self-energy expansion up to order $U^{10}$ at very low temperature. We
benchmark our results at weak and intermediate coupling with high precision
Numerical Renormalization Group (NRG) computations as well as analytical
results.
