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 $zsim6.5$ to $zsim 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 $zsim5.4$, with MIPS $24mu$m detection suggesting the presence of a powerful obscured AGN.