## 1 Introduction  
Precise Point Positioning (PPP) is an absolute Global Navigation Satellite System (GNSS) positioning method that estimates very precise positions using observations from one receiver and does not require a nearby reference station or network of stations. It requires code and phase observations on multiple frequencies and sophisticated algorithms and requires the use of precise products that are globally validated and allow homogeneous positioning across and around the Earth. Currently, it can achieve an accuracy of dm-cm in kinematic measurements and cm-mm in static measurements. However, the biggest limitation is the long convergence time, which significantly limits the method in time critical applications. Nevertheless, today PPP has proven itself as a substantial GNSS positioning method, which in addition to the position also allows for estimating parameters such as the tropospheric delay, ionospheric delay and receiver clock correction. Therefore, PPP has found an application in various scientific and commercial applications [1]. The first concepts of PPP were presented in the works [2; 3; 4]. Then, [PERSON] and [PERSON] [5] and [PERSON] and [PERSON] [6] proposed the use of precise combined International GNSS products and ionosphere-free linear combination, which eliminates the first order ionospheric delay, known as the conventional model. All these works were based on GPS observations. For the first time, the use of GPS and GLOASS was presented in [PERSON] and [PERSON] [7] and the first use all of available GNSS systems (GPS, GLONASS, Galileo, BDS) was presented in [PERSON] et al. [8]. Since then, the PPP method has developed significantly, and the most important goal is to strive for the best possible multi-GNSS, multi-frequency measurements [9; 10; 11; 12; 13; 14; 15; 16; 17]. The mathematical model of GNSS positioning is defined by means of deterministic and stochastic description [18]. The functional (deterministic) model describes the mathematical relationship between GNSS observations and the estimated parameters already known. In contrast, the stochastic model is not fully known and still requires a number of studies, especially for the PPP method. The stochastic model describes the statistical properties of observations and can be presented as variance-covariance (VC) matrix \(\mathrm{VC}=\sigma_{0}^{2}\mathrm{Q}\), where \(\sigma_{0}\) is the value of precision of observations and \(\mathrm{Q}\) is the cofactor matrix which depends on the used weighting functions [19]. Its inverse is the weight matrix \(\mathrm{P}=\mathrm{Q}^{-1}\). The full population of VC matrix contains variance and covariance elements. Variance defines the precision of observations and covariance defines their physical correlations. Defining the full VC matrix is very difficult since it depends on used models, used GNSS systems (also blocks of satellites) and their observations (frequencies and signals), used equipment (receiver and antenna types), noise, multipath and other mis-modeled errors. VC matrix also needs spatial- and time-correlation between observations [20]. Consequently, it is impossible to determine one universal model that fits all GNSS measurements and equipment. Therefore, a simple approach is followed using only the variance element of VC matrix, defining the precision of observations and weighting functions. One of the methods used for stochastic modeling of GNSS observations is the Variance Component Estimation (VCE) method [21; 22]. This method is used for the estimation of unknown variance components of the VC matrix. There are different variations of the VCE method: minimum norm quadratic unbiased estimator (MINQUE), the best invariant quadratic unbiased estimator (BIQUE), the least-squares variance component estimator (LS-VCE), the restricted maximum likelihood estimator (REML) or the Bayesian approach to VCE [23]. One of the more commonly used in GNSS measurement analysis is the LS-VCE methods [24; 25; 26; 27].