  We build a Spitzer IRAC complete catalog of objects, obtained by
complementing the $K_\mathrm{s}$-band selected UltraVISTA catalog with objects
detected in IRAC only. With the aim of identifying massive (i.e.,
$\log(M_*/M_\odot)>11$) galaxies at $4<z<7$, we consider the systematic effects
on the measured photometric redshifts from the introduction of an old and dusty
SED template and from the introduction of a bayesian prior taking into account
the brightness of the objects, as well as the systematic effects from different
star formation histories (SFHs) and from nebular emission lines in the recovery
of stellar population parameters. We show that our results are most affected by
the bayesian luminosity prior, while nebular emission lines and SFHs only
introduce a small dispersion in the measurements. Specifically, the number of
$4<z<7$ galaxies ranges from 52 to 382 depending on the adopted configuration.
Using these results we investigate, for the first time, the evolution of the
massive end of the stellar mass functions (SMFs) at $4<z<7$. Given the rarity
of very massive galaxies in the early universe, major contributions to the
total error budget come from cosmic variance and poisson noise. The SMF
obtained without the introduction of the bayesian luminosity prior does not
show any evolution from $z\sim6.5$ to $z\sim 3.5$, implying that massive
galaxies could already be present when the Universe was $\sim0.9$~Gyr old.
However, the introduction of the bayesian luminosity prior reduces the number
of $z>4$ galaxies with best fit masses $\log(M_*/M_\odot)>11$ by 83%, implying
a rapid growth of very massive galaxies in the first 1.5 Gyr of cosmic history.
From the stellar-mass complete sample, we identify one candidate of a very
massive ($\log(M_*/M_\odot)\sim11.5$), quiescent galaxy at $z\sim5.4$, with
MIPS $24\mu$m detection suggesting the presence of a powerful obscured AGN.
