The stochastic model, together with the functional model, form the mathematical model of observation that enables the estimation of the unknown parameters. In Global Navigation Satellite Systems (GNSS), the stochastic model is an especially important element as it affects not only the accuracy of the positioning model solution, but also the reliability of the carrier-phase ambiguity resolution (AR). In this paper, we study in detail the stochastic modeling problem for Multi-GNSS positioning models, for which the standard approach used so far was to adopt stochastic parameters from the Global Positioning System (GPS). The aim of this work is to develop an individual, empirical stochastic model for each signal and each satellite block for GPS, GLONASS, Galileo and BeiDou systems. The realistic stochastic model is created in the form of a fully populated variance-covariance (VC) matrix that takes into account, in addition to the Carrier-to-Noise density Ratio (C/N0)-dependent variance function, also the cross- and time-correlations between the observations. The weekly measurements from a zero-length and very short baseline are utilized to derive stochastic parameters. The impact on the AR and solution accuracy is analyzed for different positioning scenarios using the modified Kalman Filter. Comparing the positioning results obtained for the created model with respect to the results for the standard elevation-dependent model allows to conclude that the individual empirical stochastic model increases the accuracy of positioning solution and the efficiency of AR. The optimal solution is achieved for four-system Multi-GNSS solution using fully populated empirical model individual for satellite blocks, which provides a 2% increase in the effectiveness of the AR (up to 100%), an increase in the number of solutions with errors below 5 mm by 37% and a reduction in the maximum error by 6 mm compared to the Multi-GNSS solution using the elevation-dependent model with neglected measurements correlations.