The energy consumption of cities is increasing fast due to growing global population and rapid urbanization. Urban Building Energy Models (UBEMs) are promising tools to simulate the energy demand of buildings under different urban scenarios. Nowadays, the major barriers to the effective use of UBEMs are the uncertainty related to their input parameters and the lack of good-quality, open energy consumption data. The latter make deterministic UBEM simulations unreliable, and calibration approaches unsuitable for most cities in the world. The present work proposes to combine physics-based UBEMs with Uncertainty and Sensitivity Analysis on the main input parameters using aggregated data on energy use from regional/national statistics. The proposed procedure selects the most influential input parameters and characterizes their uncertainty through Forward Uncertainty Analysis and Sensitivity Analysis to obtain stochastic load profiles for space heating and cooling. The method was first tested against hourly thermal power profiles metered on a heterogeneous sample of buildings in Verona (Italy). The average heating load profile obtained is significantly improved compared to deterministic, archetype-based simulations in terms of energy needs and peak loads. The overestimation of residential buildings peak load is reduced from 80% to 25%, and the deviation in the energy needs calculation drops from 18% to 10%. The proposed simulation procedure was then applied to a district of Milan (Italy), including more than 600 buildings, resulting in similar variations. Overall, the results demonstrate that considering the uncertainty of operational, geometrical and physical parameters is of the utmost importance to obtain reliable urban simulations.