We show how consistency relations can be used to robustly extract the
amplitude of local primordial non-Gaussianity (${f}_{\mathrm{N}\mathrm{L}}$ ) from the squeezed
limit of the matter bispectrum, well into the non-linear regime. First, we
derive a non-perturbative relation between primordial non-Gaussianity and the
leading term in the squeezed bispectrum, revising some results present in the
literature. This relation is then used to successfully measure ${f}_{\mathrm{N}\mathrm{L}}$
from $N$ -body simulations. We discuss the dependence of our results on
different scale cuts and redshifts. Specifically, the analysis is strongly
dependent on the choice of the smallest soft momentum, ${q}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ , which is
the most sensitive to primordial bispectrum contributions, but is largely
independent of the choice of the largest hard momentum, ${k}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ , due to
the non-Gaussian nature of the covariance. We also show how the constraints on
${f}_{\mathrm{N}\mathrm{L}}$ improve at higher redshift, due to a reduced off-diagonal
covariance. In particular, for a simulation with ${f}_{\mathrm{N}\mathrm{L}}=100$ and a
volume of $(2.4\text{Gpc}/h{)}^{3}$ , we measure ${f}_{\mathrm{N}\mathrm{L}}=98\pm 12$ at
redshift $z=0$ and ${f}_{\mathrm{N}\mathrm{L}}=97\pm 8$ at $z=0.97$ . Finally, we compare our
results with a Fisher forecast, showing that the current version of the
analysis is satisfactorily close to the Fisher error. We regard this as a first
step towards the realistic application of consistency relations to constrain
primordial non-Gaussianity using observations.