  Upcoming surveys of cosmic structures will probe scales close to the
cosmological horizon, which opens up new opportunities for testing the
cosmological concordance model to high accuracy. In particular, constraints on
the squeezed bispectrum could rule out the single-field hypothesis during
inflation. However, the squeezed bispectrum is also sensitive to dynamical
effects of general relativity as well as interactions of matter with residual
radiation from the early Universe. In this paper, we present a relativistic
simulation pipeline that includes these relativistic effects consistently. We
produce light cones and calculate the observed number counts of cold dark
matter for five redshift bins between $z=0.55$-$2.25$. We compare the
relativistic results against reference Newtonian simulations by means of
angular power- and bispectra. We find that the dynamical relativistic effects
scale roughly inversely proportional to the multipole in the angular power
spectrum, with a maximum amplitude of $10\%$ for $\ell \lesssim 5$. By using a
smoothing method applied to the binned bispectrum we detect the Newtonian
bispectrum with very high significance. The purely relativistic part of the
matter bispectrum is detected with a significance of $\sim 3\,\sigma$, mostly
limited by cosmic variance. We find that the pure dynamical relativistic
effects accounts for up to $3\%$ and $10\%$ of the total amplitude,
respectively in the squeezed and equilateral limits. Our relativistic pipeline
for modelling ultra-large scales yields gauge-independent results as we compute
observables consistently on the past light cone, while the Newtonian treatment
employs approximations that leave some residual gauge dependence. A
gauge-invariant approach is required in order to meet the expected level of
precision of forthcoming probes of cosmic structures on ultra-large scales.
