To assess the impact of implementing energy efficiency and renewable energy measures, urban building energy models are emerging. In these models, due to the lack of data, the natural variability of the existing building stock is often highly underestimated and uncertainty on the simulated energy use arises. Therefore, this work proposes a probabilistic building characterization method to model the variability of the existing residential building stock. The method estimates realistic distributions of five input variables: U-values of the floor, external walls, windows and roof as well as window-to-wall ratio, based on known data (location, geometry and construction year). First, quantile regression has been implemented to generate the uncorrelated distributions based on the Flemish energy performance certificates database. The accuracy of the marginal distributions is good, as the empirical coverage on the 50%, 80%, 90% and 98% prediction interval deviates 0.6% at most. However, it is needed to include the correlations between these variables. Hence, three main methods to build multivariate distributions from marginal distributions and to draw correlated samples are implemented and extensively compared. The Gaussian copula method is put forward as the preferred method. Considering the mean-maximum discrepancy (MMD), this method performs eight times better than the uncorrelated case (MMD of 0.0027 versus 0.0228).