Advanced statistical methods are applied on long-term data sets from Zagreb-Grič centennial observatory in Croatia to address recent climate change and variabilities; annual averages are used. Hilbert-Huang transform (HHT), consisting of empirical mode decomposition (EMD) and Hilbert spectral analysis (HSA), is an empirically based data-analysis method for extracting periodic components embedded within generally nonlinear and non-stationary data. The EMD splits the original series into a so-called intrinsic mode functions (IMFs). First, using the EMD algorithm, intrinsic mode functions (IMFs) are obtained. The analysis of low frequency IMFs indicates significant influence of the North-Atlantic Oscillation on the air temperature, cloudiness and precipitation series in the part of northern Croatia. For the considered climatic elements, the analysis revealed mild deviations in natural fluctuations of the analyzed signal for the last 30 to 40 years, which are most likely caused by anthropogenic activities. Next, HSA is applied to each IMF obtained for the series of annual averages of sea level air pressure. In order to validate this approach and the results, an associated Hilbert spectrum (HS) is compared with the continuous wavelet analysis, i.e., this is another corresponding checkup. The comparison shows significantly improved time and frequency resolution in favor of HS. Moreover, HS provides a unique capability of displaying intra-wave frequency modulation, i.e., changes in frequency that occur within one cycle of oscillation.