Bayesian Modeling of PM₁₀ in Kraków Using Meteorological and Seasonal Factors

1. Introduction

A central difficulty in evaluating short-term ${\mathrm{PM}}_{10}$ models is that successive measurements are strongly dependent. A model containing the previous-hour concentration can appear markedly better than a meteorology-only regression even when the added meteorological variables contribute little beyond persistence. Consequently, a scientifically useful comparison should distinguish three questions: whether lagged ${\mathrm{PM}}_{10}$ dominates one-step-ahead forecasts, whether meteorology and seasonality add information beyond that lag, and whether any gain is visible in point errors, predictive distributions, or both. It should also compare the proposed model with a direct persistence forecast and with a conventional dynamic-regression benchmark under the same rolling one-step evaluation protocol.

This study addresses these questions using five years of hourly observations from Kraków. Four nested Bayesian models are considered: a lag-only model (M0), a meteorology-only model (M1), a meteorology-plus-seasonality model (M2), and a full model combining lagged ${\mathrm{PM}}_{10}$, meteorology, and cyclic seasonality (M3). They are compared with a historical-median baseline, persistence, and SARIMAX. Chronological rolling-origin validation, an untouched final test period, probabilistic scores, PSIS-LOO, posterior predictive checks, and moving-block bootstrap comparisons are used to avoid relying on a single favorable metric.

The principal contributions are:

• an interpretable ablation design that separates the contribution of lagged ${\mathrm{PM}}_{10}$, meteorological covariates, and cyclic annual and daily effects;
• a fair rolling one-step-ahead comparison with persistence and SARIMAX;
• joint assessment of point accuracy, predictive-distribution quality, interval calibration, and high-pollution performance;
• an analysis of how meteorological coefficient estimates change after short-term persistence is represented explicitly; and
• a fully reproducible Stan/Python workflow with explicit time-split, scaling, diagnostic, and sensitivity checks.

2. Materials and Methods

2.1. Data Sources and Study Period

The analysis used hourly observations from 1 January 2020 through 31 December 2024 at the Kraków Złoty Róg air-quality station (GIOŚ sensor ID 16786; station code MpKrakZloRog-PM10-1g; coordinates 50.081197^∘ N, 19.895358^∘ E). Air-quality measurements were obtained through the archival API of the Polish Chief Inspectorate of Environmental Protection (GIOŚ) [15]. Meteorological data were queried at the same coordinates from the Open-Meteo Historical Weather API [14]. The meteorological variables were temperature at 2 m, relative humidity at 2 m, wind speed at 10 m, and surface pressure.

Using the same coordinates for both sources avoided an additional spatial offset between the pollution monitor and the meteorological query point. Table 1 summarizes the study period, variables, and final split.

2.2. Data Preparation and Forecasting Information Set

Both sources were placed on a complete hourly grid and ordered chronologically. Missing ${\mathrm{PM}}_{10}$ values were not interpolated; rows without an observed target were excluded. Meteorological gaps could be forward-filled for at most two consecutive hours. Ten air-quality rows associated with five duplicated, DST-ambiguous timestamps were explicitly removed rather than silently deduplicated. The resulting modeling table contained 42,823 observations after lag construction.

The target was the one-hour-ahead concentration $y_t>0$. Forecasts for hour t used only information available by hour $t-1$: the observed $y_{t-1}$ where applicable, meteorological variables lagged by one hour, and deterministic calendar terms. Let

$$z_t=log{y}_t,x_{j,t-1}^{*}={\displaystyle \frac{x_{j,t-1}-{\overline{x}}_{j,\mathrm{train}}}{s_{j,\mathrm{train}}}},$$

(1)

where each meteorological predictor was standardized using the mean and standard deviation of the corresponding training window only. The centered lag variable was

$$z_{t-1}^c=z_{t-1}-{\overline{z}}_{\mathrm{train}}.$$

(2)

All scaling parameters were fitted separately within each cross-validation fold and once on the final 80% training period. No transformation parameter was estimated from validation or test observations.

Cyclic seasonality was represented by two annual and two daily Fourier harmonics. For period P and harmonic order $K=2$,

$$f_P\left(t\right)={\displaystyle \sum _{k=1}^2}\left[a_{k}sin\left({\displaystyle \frac{2\pi kt}{P}}\right)+b_{k}cos\left({\displaystyle \frac{2\pi kt}{P}}\right)\right],$$

(3)

with $P=24\times 365.25$ hours for the annual component and $P=24$ hours for the daily component.

2.3. Bayesian Regression Models

For each Bayesian model,

$$y_t\sim LogNormal({\eta }_t,\sigma ),$$

(4)

where ${\eta }_t$ is the model-specific predictor on the log scale and $\sigma >0$ is the residual scale.

The four Bayesian specifications were:

$$\begin{array}{cc}\hfill \mathrm{M}0:{\eta }_t& =\alpha +\varphi z_{t-1}^c,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \mathrm{M}1:{\eta }_t& =\alpha +{\beta }^{\top }{\mathbf{x}}_{t-1}^{*},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \mathrm{M}2:{\eta }_t& =\alpha +{\beta }^{\top }{\mathbf{x}}_{t-1}^{*}+f_{\mathrm{year}}\left(t\right)+f_{\mathrm{day}}\left(t\right),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathrm{M}3:{\eta }_t& =\alpha +{\beta }^{\top }{\mathbf{x}}_{t-1}^{*}+f_{\mathrm{year}}\left(t\right)+f_{\mathrm{day}}\left(t\right)+\varphi z_{t-1}^c.\hfill \end{array}$$

(8)

The persistence coefficient was constrained to $(-1,1)$ by $\varphi =tanh\left({\varphi }_{\mathrm{raw}}\right)$.

Table 2 summarizes all Bayesian and reference models.

The priors were

$$\begin{array}{cccc}\hfill \alpha & \sim \mathcal{N}(log30,0.7),\hfill & \hfill {\beta }_j& \sim \mathcal{N}(0,0.3),\hfill \end{array}$$

(9)

$$\begin{array}{cccc}\hfill a_k,b_k,c_k,d_k& \sim \mathcal{N}(0,0.3),\hfill & \hfill {\varphi }_{\mathrm{raw}}& \sim \mathcal{N}(0,0.5),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \sigma & \sim HalfNormal(0,0.5).\hfill \end{array}$$

(11)

Prior predictive simulations were used to check that these distributions permitted a broad but plausible range of concentrations.

2.4. Reference Models

B0 predicted the training-period median for every evaluation hour. B1 was deterministic persistence, ${\widehat{y}}_t=y_{t-1}$. B2 was a SARIMAX dynamic regression fitted to $z_t$ with an AR(1) component and the same lagged meteorological and Fourier regressors used by the Bayesian specifications, implemented with Statsmodels [27]. SARIMAX parameters were fitted once within each training window. Evaluation then proceeded one hour at a time: a one-step forecast was generated, the observed value was appended to the state without parameter refitting, and the next forecast was produced. This matched the information available to persistence and the lag-containing Bayesian models.

Figure 1 summarizes the nested Bayesian model sequence through directed dependency diagrams. The diagrams are explanatory graphical summaries of each model specification and are not intended as causal-identification assumptions.

2.5. Inference and Diagnostics

[21,22]. Four chains were run, each with 1000 warm-up and 1000 retained iterations, target acceptance probability 0.90, and maximum tree depth 12. Sampling quality was assessed using rank-normalized $\widehat{R}$, bulk and tail effective sample sizes, divergent transitions, maximum-tree-depth hits, and E-BFMI. Posterior processing and PSIS-LOO calculations used ArviZ [23].

2.6. Chronological Validation and Predictive Metrics

The data were split chronologically. The first 80% formed the final training set and the last 20% an untouched held-out test set. Model development on the first 80% used four expanding-window rolling-origin folds. The initial training window covered the first 60% of all observations, and each validation window covered the next 5%; later folds incorporated earlier validation windows into training. Figure 2 shows the time series and all split boundaries.

Point forecasts were evaluated using MAE and RMSE. Probabilistic forecasts were evaluated using the continuous ranked probability score (CRPS), empirical coverage of the 50% and 90% predictive intervals, and mean interval width. CRPS is a proper scoring rule that rewards both calibration and concentration of a predictive distribution [25]. The deterministic B0 and B1 forecasts were additionally treated as degenerate distributions when reporting CRPS; they do not provide genuine interval forecasts.

For Bayesian models, PSIS-LOO and WAIC were calculated as complementary measures of full predictive-density fit [24]. Because ordinary pointwise LOO does not reproduce a future-block forecasting task, chronological cross-validation and the final held-out test remained the primary evidence for out-of-time performance.

Pairwise uncertainty in metric differences was assessed using a paired moving-block bootstrap [26]. Two thousand resamples were generated with a primary block length of 24 hours and a sensitivity block length of 168 hours. High-pollution performance was assessed for test observations above the 90th percentile of the training-period ${\mathrm{PM}}_{10}$ distribution (58.7 $\mu \mathrm{g}{\mathrm{m}}^{-3}$).

3. Results

3.1. Data Characteristics and Prior Predictive Assessment

The final modeling dataset contained 42,823 hourly observations. The final training set contained 34,258 observations, and the held-out test set contained 8565 observations from 31 December 2023 19:00 through 31 December 2024 23:00. The observed series displayed pronounced right skewness, seasonal variation, and episodic peaks, including concentrations above 200 $\mu \mathrm{g}{\mathrm{m}}^{-3}$ in several winter periods (Figure 2).

Detailed prior predictive checks are provided in Supplementary Figure S1.

3.2. Sampling Diagnostics

All final Bayesian fits satisfied the diagnostic criteria. Maximum $\widehat{R}$ was 1.000 for M0–M3, minimum bulk ESS exceeded 2000, minimum tail ESS exceeded 2200, and no divergent transitions or maximum-tree-depth hits occurred. Minimum E-BFMI ranged from 0.868 to 1.008. Full values are reported in Supplementary Table S1.

3.3. Overall Predictive Performance

Table 3 and Figure 3 summarize final test-set performance. The historical-median baseline was clearly inadequate (MAE 13.013 $\mu \mathrm{g}{\mathrm{m}}^{-3}$). Models without lagged ${\mathrm{PM}}_{10}$ also performed poorly: M1 and M2 had MAE values of 11.605 and 11.374 $\mu \mathrm{g}{\mathrm{m}}^{-3}$, respectively. Adding cyclic seasonality to meteorology alone yielded only a modest improvement.

The four persistence-aware models were tightly grouped. M3 achieved the lowest MAE (4.865 $\mu \mathrm{g}{\mathrm{m}}^{-3}$) and CRPS (3.602 $\mu \mathrm{g}{\mathrm{m}}^{-3}$). Persistence had MAE 4.886 $\mu \mathrm{g}{\mathrm{m}}^{-3}$, SARIMAX had MAE 4.884 $\mu \mathrm{g}{\mathrm{m}}^{-3}$ and CRPS 3.622 $\mu \mathrm{g}{\mathrm{m}}^{-3}$, and M0 had MAE 4.930 $\mu \mathrm{g}{\mathrm{m}}^{-3}$ and CRPS 3.659 $\mu \mathrm{g}{\mathrm{m}}^{-3}$. The lowest RMSE was obtained by persistence (7.320 $\mu \mathrm{g}{\mathrm{m}}^{-3}$), closely followed by M0 (7.332), SARIMAX (7.426), and M3 (7.453). Thus, the ranking differed slightly between absolute-error, squared-error, and distributional criteria.

The probabilistic intervals were reasonably calibrated. Coverage of nominal 90% intervals was 0.885 for SARIMAX, 0.889 for M0, 0.873 for M1, 0.872 for M2, and 0.884 for M3. M3 produced the narrowest 90% interval among the probabilistic persistence-aware models (22.957 $\mu \mathrm{g}{\mathrm{m}}^{-3}$), marginally narrower than SARIMAX and M0.

The rolling-origin validation folds showed the same qualitative pattern: M3 had the best mean MAE, RMSE, and CRPS among the probabilistic models, followed closely by SARIMAX and M0, while M1 and M2 were substantially worse. This agreement reduces the likelihood that the final ranking resulted solely from a favorable final-year realization.

3.4. Uncertainty in Differences Between Models

The point differences among M3, persistence, M0, and SARIMAX were small relative to temporally dependent variation in forecast errors. Table 4 reports paired moving-block bootstrap comparisons for the 24-hour blocks. None of the 95% intervals for MAE or RMSE differences excluded zero. For example, the observed M3-minus-persistence MAE difference was $-0.020$ $\mu \mathrm{g}{\mathrm{m}}^{-3}$, with a 95% interval of $[-0.141,0.124]$.

For CRPS, the evidence favored M3 more consistently but remained slightly below a strict 95% threshold on the single held-out year. Relative to M0, the observed CRPS difference was $-0.057$ $\mu \mathrm{g}{\mathrm{m}}^{-3}$, with $P(\Delta <0)=0.941$ and a 95% interval of $[-0.118,0.015]$. Relative to SARIMAX, the difference was $-0.020$ $\mu \mathrm{g}{\mathrm{m}}^{-3}$, with $P(\Delta <0)=0.923$ and interval $[-0.046,0.008]$. The 168-hour block-length sensitivity analysis led to the same substantive conclusion.

3.5. Predictive Density and Posterior Predictive Behavior

PSIS-LOO and WAIC provided a clearer separation of the predictive distributions than the final-year point metrics. M3 had the highest expected log pointwise predictive density, exceeding M0 by 1321.8 elpd units with a standard error of the difference of 52.1. M1 and M2 were more than 24,000 elpd units behind M3. Figure 4 shows the PSIS-LOO ranking relative to M3 together with the associated Pareto-k diagnostics. All maximum Pareto-k values were below 0.24, indicating stable PSIS approximations. Full numerical results are provided in Supplementary Table S2.

Posterior predictive checks also favored M3. The observed training-data 99th percentile was 100.6 $\mu \mathrm{g}{\mathrm{m}}^{-3}$; the median replicated value was 99.9 for M3, compared with 102.2 for M0, 108.7 for M1, and 111.3 for M2. The lag-containing models also left substantially less short-lag autocorrelation in test residuals than the meteorology-only models (Supplementary Figure S2).

Figure 5 illustrates M3 and persistence during a representative winter test interval. The predictive median tracked the evolution of the episode closely, while the predictive interval widened during the sharp increase. The largest observed peak, 227 $\mu \mathrm{g}{\mathrm{m}}^{-3}$, remained difficult to anticipate precisely, reflecting the challenge of sudden extreme episodes.

3.6. Meteorological Effects and Temporal Persistence

$$100[exp(c)-1]\%.$$

(12)

Posterior intervals were calculated after transforming every posterior draw rather than by transforming only a point estimate.

The contrast between M1 and M3 was substantial (Table 5, Figure 6). In M1, a one-standard-deviation increase in lagged temperature was associated with a 22.6% decrease in median ${\mathrm{PM}}_{10}$, and wind speed with a 23.4% decrease. After lagged ${\mathrm{PM}}_{10}$ was included, these effects shrank to 2.0% and 4.2%, respectively. The humidity effect shrank from $-1.6$% to $-0.7$%, and the surface-pressure effect from $+11.3$% to $+1.7$%.

The autoregressive coefficient was 0.914 (90% interval 0.910–0.917) in M0 and 0.865 (0.860–0.869) in M3. The lower M3 coefficient is consistent with meteorology and cyclic terms explaining a portion of variation that otherwise appears as pure persistence, although the lag term remains dominant.

3.7. Performance During High-Pollution Hours

The high-pollution subset contained 636 test hours above 58.7 $\mu \mathrm{g}{\mathrm{m}}^{-3}$. All models had larger errors in this regime, but the relative ranking changed. Persistence was the best point predictor, with MAE 10.4 and RMSE 14.6 $\mu \mathrm{g}{\mathrm{m}}^{-3}$. M0, M3, and SARIMAX had MAE values of 12.8, 13.0, and 13.2 $\mu \mathrm{g}{\mathrm{m}}^{-3}$, respectively (Figure 7). The meteorology-only models deteriorated much more strongly, with MAE above 41 $\mu \mathrm{g}{\mathrm{m}}^{-3}$. Therefore, M3’s overall distributional advantage should not be interpreted as uniformly superior point prediction during extreme episodes.

3.8. Sensitivity Analysis

Widening the meteorological coefficient prior from $\mathcal{N}(0,0.3)$ to $\mathcal{N}(0,0.6)$ had negligible effects on predictions and coefficient intervals. Replacing the lognormal model with a Student-t likelihood on log-${\mathrm{PM}}_{10}$ reduced MAE and CRPS slightly but increased RMSE by 28.5%, indicating greater tolerance of typical errors at the cost of occasional large misses. These results did not overturn the main ranking but show that tail-sensitive conclusions depend on the likelihood family (Supplementary Figure S4 and Tables S3–S4).

4. Discussion

4.1. Temporal Persistence as the Dominant Short-Horizon Signal

The primary empirical result is that one-hour-ahead point forecasting is dominated by the most recent ${\mathrm{PM}}_{10}$ observation. M1 and M2, which excluded the lag term, were markedly inferior to persistence-aware models. By contrast, M3, M0, persistence, and SARIMAX differed by only hundredths of a $\mu \mathrm{g}{\mathrm{m}}^{-3}$ in MAE. The block bootstrap showed that these point-error differences were not distinguishable from zero at the 95% level.

In the present analysis, however, the explicit M0 and persistence comparisons show that this improvement should not be attributed automatically to the meteorological part of a dynamic model. Much of the attainable one-hour point accuracy is already captured by persistence.

The distinction is important methodologically. Comparing M3 only with M1 or M2 would suggest a dramatic improvement from the full Bayesian specification. Comparing it with persistence and M0 reveals a more modest incremental contribution. This illustrates why a direct persistence baseline and matched rolling one-step protocol are essential in short-horizon environmental forecasting.

4.2. Where Meteorology and Seasonality Add Value

Meteorology and cyclic seasonality provided limited incremental value in MAE and RMSE but more visible value in the full predictive distribution. M3 achieved the lowest test CRPS, the highest PSIS-LOO elpd, good interval calibration, and the closest posterior reproduction of the observed upper-tail quantile. These gains are meaningful for applications that require exceedance probabilities or uncertainty intervals rather than only a single central prediction.

The bootstrap evidence for CRPS was directionally strong but not conclusive under a strict 95% criterion: $P(\Delta <0)$ was approximately 0.94 against M0 and 0.92 against SARIMAX. PSIS-LOO produced a much clearer separation because it aggregated full log predictive densities across the training observations. Nevertheless, the chronological test and rolling-origin results remain the more direct evidence for future-time forecasting; LOO should be interpreted as complementary evidence about predictive-density fit rather than a replacement for time-respecting validation.

4.3. Interpretation of Meteorological Associations

The magnitude of these associations depended strongly on whether temporal persistence was included. The apparent M1 temperature and wind effects were reduced by approximately an order of magnitude in M3. This attenuation indicates that meteorological variables in a static model partly proxy for persistent pollution episodes and their seasonal context. The M3 effects remain conditional associations, not causal effects: unmodeled emissions, boundary-layer dynamics, traffic, and heating demand can influence both meteorology and pollution.

Replacing four categorical seasons from the first manuscript version with cyclic Fourier terms also avoided artificial discontinuities at season boundaries and represented both daily and annual periodicity.

4.4. High-Pollution Episodes

The high-pollution analysis changes the practical interpretation. Persistence was the best point predictor above the training 90th percentile, while M3 had higher MAE and RMSE in that subset. Sudden peaks are difficult for smooth probabilistic regressions because the preceding hour may be more informative than meteorological variables that change gradually. This does not negate M3’s overall calibration and distributional performance, but it shows that an operational alert system should not assume uniform improvement across regimes.

Future operational extensions could model exceedance events or regime changes directly, incorporate emission-related variables, or combine calibrated probabilistic forecasts with a persistence-based safeguard during rapidly developing episodes. Such developments would constitute new experiments and are outside the scope of the present study.

4.5. Strengths and Limitations

The nested M0–M3 design further separates the contributions of persistence, meteorology, and seasonality, and the validation framework prevents the final test year from influencing preprocessing or model selection.

The study also has limitations. First, all observations come from a single Kraków station, so the estimates may reflect local emission sources, urban structure, topography, and sensor context. The results cannot be assumed to generalize directly to other cities or monitoring stations without recalibration. Second, the meteorological predictors were obtained from a gridded historical-weather service rather than a co-located physical weather station. Third, traffic, residential heating, industrial activity, boundary-layer height, precipitation, and gaseous co-pollutants were not included. Fourth, the moving-block bootstrap used one held-out year; it accounts for serial dependence within that year but does not provide multiple independent annual test replications. Finally, the Student-t sensitivity result indicates that tail-related conclusions depend partly on the likelihood specification.

5. Conclusions

The corrected five-year analysis shows that lagged ${\mathrm{PM}}_{10}$ is the dominant source of one-hour-ahead point accuracy. Models without the lag term were substantially worse, whereas M3, M0, persistence, and SARIMAX had statistically indistinguishable MAE and RMSE on the final test year.

The full M3 specification nevertheless provided the strongest overall probabilistic evidence: it had the lowest test CRPS, the best PSIS-LOO elpd, calibrated 90% intervals, and posterior predictive upper-tail behavior closest to the observed data. Meteorological effects became much smaller after lagged ${\mathrm{PM}}_{10}$ was introduced, demonstrating that static meteorology-only regressions can overstate associations by absorbing short-term persistence.

At the same time, persistence remained the best point predictor during high-pollution hours, so the full model should not be presented as uniformly superior across all regimes.

The main practical conclusion is therefore nuanced: for one-hour-ahead forecasts, meteorology and cyclic seasonality add little to central point accuracy beyond the most recent concentration, but they add clearer value to the quality and interpretation of the predictive distribution. External validation at additional stations and explicit modeling of high-pollution regimes are the most relevant directions for future work.