The rich phenomenology of twisted bilayer graphene (TBG) near the magic angle is believed to arise from electron correlations in topological flat bands. An unbiased approach to this problem is highly desirable, but also particularly challenging, given the multiple electron flavors, the topological obstruction to defining tight binding models and the long-ranged Coulomb interactions. While numerical simulations of realistic models have thus far been confined to zero temperature, typically excluding some spin or valley species, analytic progress has relied on fixed point models away from the realistic limit. Here we present for the first time unbiased Monte Carlo simulations of realistic models of magic angle TBG at charge-neutrality. We establish the absence of a sign problem for this model in a momentum space approach, and describe a computationally tractable formulation that applies even on breaking chiral symmetry and including band dispersion. Our results include (i) the emergence of an insulating Kramers inter-valley coherent ground state in competition with a correlated semi-metal phase, (ii) detailed temperature evolution of order parameters and electronic spectral functions which reveal a `pseudogap regime, in which gap features are established at a higher temperature than the onset of order and (iii) predictions for electronic tunneling spectra and their evolution with temperature. Our results pave the way towards uncovering the physics of magic angle graphene through exact simulations of over a hundred electrons across a wide temperature range.