For a number of applications like low-source reactor start-up or neutron coincidence counting, it is necessary to take into account the stochastic nature of neutron transport and go beyond the average neutron density, which is the solution to the linear Boltzmann equation. In this work, we are particularly interested in calculating the moments and probabilities of the number of neutrons detected during a time window. It is known that in the case of a single initial neutron, these quantities are the solution of a system of coupled adjoint transport equations and that a neutron source can be taken into account a second time using summations with the source strength. The purpose of the current work is first to present the derivation of these equations in a form where they can be solved up to an arbitrary order and where the fission terms are expanded according to the moments of the fission multiplicity. For simplicity, the derivation is first presented in a lumped, also called point, model approximation before considering the full phase-space case. Thereafter, we describe the implementation of the solution algorithm in the deterministic, discrete ordinates code PANDA. Finally, we present some numerical results and intercode comparisons for verification purposes and to illustrate the applicability of the method.