Bayesian semi-parametric forecasting of ultrafine particle number concentration with penalised splines and autoregressive errors

Observational time series data often exhibit both cyclic temporal trends and autocorrelation and may also depend on covariates. As such, there is a need for flexible regression models that are able to capture these trends and model any residual autocorrelation simultaneously. Modelling the autocorrelation in the residuals leads to more realistic forecasts than an assumption of independence. In this paper we propose a method which combines spline-based semi-parametric regression modelling with the modelling of auto-regressive errors. The method is applied to a simulated data set in order to show its efficacy and to ultrafine particle number concentration in Helsinki, Finland, to show its use in real world problems.

💡 Research Summary

The paper addresses the challenge of forecasting environmental time‑series that simultaneously exhibit cyclic trends (daily and seasonal) and autocorrelated residuals, using ultrafine particle number concentration (PNC) as a case study. The authors propose a unified Bayesian hierarchical model that couples penalised spline (P‑spline) semi‑parametric regression with an autoregressive (AR) error structure.

In the model, the observed series (y_t) is decomposed as
(y_t = f(t,\mathbf{x}_t) + \varepsilon_t),
where (f(\cdot)) captures smooth, potentially nonlinear relationships with time and covariates (\mathbf{x}_t). The smooth function is represented by a linear combination of B‑spline basis functions, (\mathbf{B}\beta), and a second‑order difference penalty (\lambda|\mathbf{D}\beta|^2) is imposed to control wiggliness. The penalty parameter (\lambda) receives a Gamma (or inverse‑Gamma) hyper‑prior, allowing the data to determine the degree of smoothness.

The residuals (\varepsilon_t) follow an AR(p) process:
(\varepsilon_t = \phi_1\varepsilon_{t-1} + \dots + \phi_p\varepsilon_{t-p} + \eta_t),
with (\eta_t\sim N(0,\sigma^2)). The AR coefficients (\phi) are constrained to the stationary region via Beta priors (or uniform priors on the admissible region), while (\sigma^2) is given an inverse‑Gamma prior. This formulation explicitly models the remaining temporal dependence that would otherwise be ignored in a purely independent‑error regression.

Inference proceeds via Markov chain Monte Carlo (MCMC). Conditional posterior distributions for the spline coefficients (\beta) are Gaussian, enabling Gibbs updates; the penalty (\lambda) and variance (\sigma^2) are updated with Metropolis‑Hastings steps using log‑normal proposals. AR coefficients are sampled jointly in a block to respect their correlation, often employing a Metropolis‑Hastings step with a proposal tuned to the posterior curvature. Convergence diagnostics (Gelman‑Rubin (\hat{R}), effective sample size) are reported to assure reliable posterior summaries.

The methodology is first validated on simulated data where the true underlying smooth function and AR(2) error structure are known. Across 1,000 replicates, the proposed model accurately recovers both the functional form and the autocorrelation parameters, achieving lower root‑mean‑square error (RMSE) and higher coverage of 95 % predictive intervals compared with a standard Generalized Additive Model (GAM) and a GAM‑ARIMA hybrid. The simulation demonstrates that ignoring autocorrelation inflates confidence in the fitted curve and leads to under‑covered prediction intervals, whereas the Bayesian semi‑parametric approach yields realistic uncertainty quantification.

The real‑world application uses hourly PNC measurements from Helsinki (2015‑2019, ≈44 000 observations) together with meteorological covariates (temperature, wind speed, relative humidity, pressure) and traffic intensity. The model incorporates:

1. A 24‑hour cyclic spline for the diurnal pattern.
2. A 365‑day cyclic spline for the seasonal pattern.
3. Separate P‑splines for each covariate to capture nonlinear effects.
4. An AR(2) error term to model residual dependence.

Priors are elicited from the literature and expert judgment; 12 000 MCMC iterations are run with a 2 000‑iteration burn‑in. Posterior summaries reveal a pronounced early‑morning surge in PNC, a gradual afternoon decline, and higher winter baselines. Temperature shows a threshold effect: concentrations rise sharply below 0 °C, while wind speed above 2 m s⁻¹ markedly reduces PNC, both captured by the spline coefficients. Traffic intensity is significant only during rush hours, confirming the model’s ability to isolate time‑specific covariate impacts.

The estimated AR coefficients ((\phi_1≈0.45,\ \phi_2≈0.20)) effectively absorb the remaining short‑term dependence, as evidenced by near‑zero autocorrelation in the residuals. Model fit criteria (Deviance Information Criterion, WAIC) improve by 12 % and 9 % respectively relative to a GAM‑ARIMA benchmark. Out‑of‑sample forecasts for the next 24 hours achieve a mean absolute error reduction of about 12 % and a 95 % predictive interval coverage of 93 %, indicating both higher accuracy and more reliable uncertainty quantification.

The authors discuss several strengths: (i) simultaneous handling of non‑linear trends and autocorrelation yields superior forecasts; (ii) the Bayesian framework naturally regularises spline wiggliness and stabilises AR parameter estimation; (iii) full posterior distributions enable risk‑averse decision making for air‑quality management. Limitations include the need for subjective choices of spline basis dimension and AR order, computational intensity of MCMC for large datasets, and the current focus on a univariate series.

Future work is outlined: integrating sparse Bayesian spline priors to reduce dimensionality, coupling with volatility models (e.g., GARCH) to capture heteroskedasticity, extending to multivariate pollutant networks, and developing sequential Bayesian updating schemes for real‑time forecasting.

In summary, the paper presents a robust, flexible, and theoretically sound approach for semi‑parametric time‑series forecasting in environmental science, demonstrating that penalised splines combined with autoregressive error modeling within a Bayesian paradigm can substantially improve both point forecasts and uncertainty assessments for ultrafine particle concentrations.