Within this research, we tested the efficiency of subtraction of the loading models i.e. atmospheric, hydrologic and non-tidal ocean loadings from the Global Navigation Satellite System (GNSS) vertical records. It is widely acknowledged, that the loading effects should be removed from the GNSS time series prior to the estimates of the horizontal and vertical velocities. On the one hand, the loading models can be removed directly from the GNSS data with no interest in leaving the spectral properties of GNSS time series intact. Other ways, the loading models can be removed in more sophisticated way than a direct subtraction. We employed data from more than 350 permanent IGS (International GNSS Service) stations, derived as the official contribution to ITRF2014 (International Terrestrial Reference Frame). We proposed the method of optimal subtraction of loading effects by applying the Improved Singular Spectrum Analysis. The advantage in the proposed method lies in no artificial loss of energy in the Power Spectrum Density between 4 and 80 cpy observed comparing to the case when environmental loadings have been subtracted directly from the position time series. Finally, we estimated Dilution of Precision (DP) of the vertical velocities of the ITRF2014 series. For 81% of analysed stations, the DP was lower than 1, which means the improvement in the vertical velocity uncertainty due to more efficient modelling of time-varying seasonal signals in the loadings.