  We study primordial non-Gaussian signatures in the redshift-space halo field
on non-linear scales, using a quasi-maximum likelihood estimator based on
optimally compressed power spectrum and modal bispectrum statistics. We train
and validate the estimator on a suite of halo catalogues constructed from the
Quijote-PNG N-body simulations, which we release to accompany this paper. We
verify its unbiasedness and near optimality, for the three main types of
primordial non-Gaussianity (PNG): local, equilateral, and orthogonal. We
compare the modal bispectrum expansion with a $k$-binning approach, showing
that the former allows for faster convergence of numerical derivatives in the
computation of the score-function, thus leading to better final constraints. We
find, in agreement with previous studies, that the local PNG signal in the
halo-field is dominated by the scale-dependent bias signature on large scales
and saturates at $k \sim 0.2~h\,\mathrm{Mpc}^{-1}$, whereas the small-scale
bispectrum is the main source of information for equilateral and orthogonal
PNG. Combining power spectrum and bispectrum on non-linear scales plays an
important role in breaking degeneracies between cosmological and PNG
parameters; such degeneracies remain however strong for equilateral PNG. We
forecast that PNG parameters can be constrained with $\Delta
f_\mathrm{NL}^\mathrm{local} = 45$, $\Delta f_\mathrm{NL}^\mathrm{equil} =
570$, $\Delta f_\mathrm{NL}^\mathrm{ortho} = 110$, on a cubic volume of $1
\left({ {\rm Gpc}/{ {\rm h}}} \right)^3$, at $z = 1$, considering scales up to
$k_\mathrm{max} = 0.5~h\,\mathrm{Mpc}^{-1}$.
