Uncertainty Quantification of Next Generation Reservoir Computing with Applications to Memory-Driven Dynamical Systems

Livia Popa; Sumanta Basu; Martin T. Wells

stat.MLSep 26, 2026

Uncertainty Quantification of Next Generation Reservoir Computing with Applications to Memory-Driven Dynamical Systems

Authors: Livia Popa, Sumanta Basu, Martin T. Wells

Organizations: Department of Statistics and Data Science, Cornell University, New York, USA

Abstract

Nonlinear dynamical systems with memory arise across science and engineering, yet uncertainty quantification for efficient forecasting methods such as Next Generation Reservoir Computing (NGRC) remains underdeveloped. We study Bayesian ridge and conformal prediction intervals for NGRC and characterize when their uncertainty estimates agree or differ. In low dimensions, their asymptotic widths are governed by different summaries of the residual distribution, so agreement depends on residual shape rather than dimensionality alone. In high dimensions, regularization introduces a further tradeoff between estimation variance, shrinkage bias, and posterior uncertainty, leading to an explicit transition between regimes where Bayesian intervals are wider or narrower than conformal intervals. We extend these results to quadratic NGRC feature maps and give sufficient conditions for transferring the analysis to temporally dependent forecast windows. Simulations and real-data experiments support the theoretical predictions and illustrate how residual distribution, regularization, dimensionality, and distribution shift affect interval calibration and efficiency. These results provide a principled framework for choosing and interpreting uncertainty quantification methods in reservoir-based forecasting.

Figures & tables

Figure 1: Fixed-dimensional residual-shape effects. Left: Bayesian–conformal width difference and theoretical limit. Middle: Bayesian coverage and its theoretical limit. Right: split-conformal coverage; dashed lines denote nominal coverage.

Figure 2: High-dimensional width and coverage phase transition. Top: Marchenko–Pastur limits for $W_{B}^{2}-W_{C}^{2}$ and the equal-width boundary. Middle: empirical squared-width differences and finite-trace approximations. Bottom: empirical and theoretical Bayesian coverage; dotted lines mark 95% coverage.

Figure 3: Quadratic NGRC diffuseness and temporal transfer. Top: ridge-transformed diffuseness, residual Gaussian approximation, and finite-trace width error. Bottom: feature-design comparisons and convergence toward the independent-window benchmark with increasing forecast-origin separation.

Figure 4: Online response to an observation-noise variance change in the nonlinear Volterra system. The vertical dashed line marks the increase in observation-noise standard deviation from $0.15$ to $0.35$ at test step 160, and the horizontal dotted line denotes nominal 95% coverage. Panel (a) reports trailing-window coverage and panel (b) reports mean interval width. Curves are averages over 100 paired replications and shaded regions are 95% Monte Carlo confidence intervals.

Method	Constant coverage	Constant width	First-50 high-noise coverage	Full high-noise coverage	High-noise width	Recovery fraction	Observations from noise increase to recovery
Bayesian (frozen scale)	$0.941\pm 0.001$	$0.624\pm 0.001$	$0.702\pm 0.005$	$0.726\pm 0.001$	$0.858\pm 0.003$	$0/100$	–
SCP (frozen quantile)	$0.949\pm 0.001$	$0.645\pm 0.002$	$0.606\pm 0.008$	$0.595\pm 0.003$	$0.645\pm 0.002$	$0/100$	–
Rolling SCP	$0.950\pm 0.000$	$0.656\pm 0.002$	$0.794\pm 0.004$	$0.935\pm 0.001$	$1.502\pm 0.006$	$62/100$	$592.6\pm 5.1$
Updated Gaussian	$0.947\pm 0.001$	$0.643\pm 0.002$	$0.756\pm 0.005$	$0.929\pm 0.001$	$1.458\pm 0.005$	$58/100$	$587.2\pm 5.5$
Adaptive CP	$0.949\pm 0.000$	$0.650\pm 0.002$	$0.832\pm 0.003$	$0.938\pm 0.000$	$1.544\pm 0.007$	$80/100$	$546.6\pm 9.8$
Time-weighted CP	$0.945\pm 0.000$	$0.642\pm 0.002$	$0.828\pm 0.004$	$0.933\pm 0.000$	$1.476\pm 0.006$	$76/100$	$574.4\pm 7.6$

Table 1: Performance under constant observation-noise variance and after the noise increase in the nonlinear Volterra system. First-50 high-noise coverage is calculated over test observations 160–209. Full high-noise coverage and high-noise width use test observations 160–799. Recovery fraction is the number of replications satisfying the sustained-recovery criterion. Observations to recovery counts the number of test observations from the noise increase until sustained recovery and is reported as mean

\pm

standard error among replications that recover. All other entries are means

\pm

standard errors over 100 paired replications.

Figure 5: Fitting–calibration tradeoff under the fixed total budget $B=400$ . Bayesian ridge uses all observations for fitting, while split conformal varies the fitting fraction $\rho$ . Columns show coverage, interval width, Winkler score, and test RMSE for $p=10$ and $p=300$ . Results are means over 200 paired replications with 95% confidence intervals.

Dataset	$n_{\rm fit}$	$p/n_{\rm fit}$	$d_{\rm eff}/n_{\rm fit}$	$\widehat{\tau}_{\rm fit}$	$\widehat{\tau}_{\rm cal}$	$\widehat{\tau}_{\rm test}$	$\widehat{q}_{.95,\rm cal}$	$\widehat{q}/(z\widehat{\tau}_{\rm fit})$
LA ozone	6420	0.0394	0.0287	0.521	0.392	0.360	0.787	0.771
Solar	21009	0.0338	0.0038	0.029	0.026	0.031	0.057	1.001
Beijing PM 10	14016	0.0238	0.0116	0.231	0.248	0.182	0.388	0.858
Exchange returns	3029	0.0396	0.0038	0.993	1.445	1.284	2.704	1.389

Table 2: Theory-facing diagnostics for the one-step, 95% fixed-origin analysis. RMS scales and residual quantiles are on the standardized-response scale.

Dataset	Shape factor	Scale factor	$\rho_{q}$
LA ozone	1.024	0.752	0.771
Solar	1.119	0.897	1.001
Beijing PM 10	0.798	1.074	0.858
Exchange returns	0.955	1.455	1.389

Table 3: Decomposition of the residual diagnostic at 95% nominal coverage.

Dataset	Method	RMSE	Coverage	$\Delta_{\rm cov}$ (pp)	Width	Winkler	UQ time (s)
LA ozone	Bayesian ridge	0.360	0.993	4.25	2.047	2.083	0.015
	SCP		0.963	1.29	1.573	1.795	$<0.001$
	Asymmetric SCP		0.963	1.29	1.574	1.791	0.001
	Adaptive CP		0.954	0.36	1.485	1.784	6.523
	Time-weighted CP		0.954	0.36	1.477	1.782	8.677
	ResCP		0.929	-2.13	1.369	1.797	18.831

Table 4: Fixed-origin one-step NGRC uncertainty-quantification results for the 40% fit/40% calibration/20% test split at 95% nominal coverage. All methods within a dataset use the same NGRC point forecasts, so RMSE is reported once per dataset. The signed coverage error is

\Delta_{\rm cov}=100(\widehat{\mathrm{Cov}}-0.95)

percentage points, computed from unrounded empirical coverage. UQ time is the wall-clock time for method-specific calibration, tuning, sequential updating, and interval construction; it excludes the shared NGRC fit. Boldface marks the coverage closest to 0.95 and the lowest Winkler score within each dataset.

Figure 6: Coverage error versus mean interval width (top) and Winkler score (bottom) at 95% nominal coverage. The vertical line marks exact marginal coverage. All methods use the same NGRC point forecasts.

Dataset	Forecasting pipeline	RMSE	Coverage	PI width	Winkler	Total runtime (s)
LA ozone	ARIMA + ResCP	0.359	0.941	1.439	1.834	26.190
	GRU + ResCP	0.365	0.956	1.524	1.810	60.960
	Transformer + ResCP	0.364	0.935	1.391	1.841	141.942
	NGRC + ResCP	0.360	0.929	1.369	1.797	19.390
Solar	ARIMA + ResCP	0.021	0.904	0.038	0.069	156.337
	GRU + ResCP	0.020	0.894	0.041	0.070	442.007

Table 5: Comparison of point-forecasting architectures under the common official ResCP uncertainty procedure for one-step-ahead forecasting at 95% nominal coverage. Runtime includes preprocessing, point-model fitting, forecasting, ResCP tuning, and interval construction.

Appendix figures & tables14 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 7: Individual finite-trace width verification for Experiment 2. Points compare cell-average empirical squared widths with their finite-trace approximations for Bayesian intervals (left) and split-conformal intervals (right). Error bars are 95% confidence intervals for the cell means over 100 replications. Color denotes the aspect ratio $\gamma=p/n$ and markers in the split-conformal panel denote signal-to-noise ratio. The cells span $\lambda\in\{0.05,0.2,1\}$ . The dashed line denotes equality.

$d$	Signal structure	$\\|M_{d}\\|_{\rm op}$	$\\|D_{n}\\|_{\rm op}$	Residual KS
8	Concentrated	1.414	$0.973\pm 0.013$	$0.140\pm 0.001$
8	Diffuse	0.500	$0.545\pm 0.004$	$0.079\pm 0.001$
12	Concentrated	1.414	$0.941\pm 0.010$	$0.133\pm 0.001$
12	Diffuse	0.408	$0.444\pm 0.002$	$0.065\pm 0.001$
16	Concentrated	1.414	$0.923\pm 0.008$	$0.130\pm 0.001$
16	Diffuse	0.354	$0.386\pm 0.002$	$0.059\pm 0.001$

Appendix

Table 6: Raw and ridge-transformed quadratic diffuseness diagnostics for Experiment 3. Entries are means

\pm

standard errors over 150 replications. The generating-matrix norm refers to the unnormalized matrix

M_{d}

G^{\top}M_{d}G-\operatorname{tr}(M_{d})

Feature design	Signal	Empirical $W_{B}$	Trace $W_{B}$	Empirical $W_{C}$	Trace $W_{C}$	Bayesian cov.	SCP cov.	Residual KS
Independent Gaussian	Concentrated	$1.165\pm 0.000$	$1.165\pm 0.000$	$3.732\pm 0.016$	$3.715\pm 0.010$	$0.462\pm 0.001$	$0.950\pm 0.001$	$0.019\pm 0.000$
Independent Gaussian	Diffuse	$1.165\pm 0.000$	$1.166\pm 0.000$	$3.743\pm 0.015$	$3.747\pm 0.011$	$0.459\pm 0.002$	$0.949\pm 0.001$	$0.020\pm 0.000$
Quadratic Gaussian	Concentrated	$1.169\pm 0.000$	$1.173\pm 0.000$	$3.843\pm 0.035$	$3.870\pm 0.023$	$0.542\pm 0.002$	$0.950\pm 0.001$	$0.125\pm 0.001$
Quadratic Gaussian	Diffuse	$1.169\pm 0.000$	$1.173\pm 0.000$	$3.830\pm 0.015$	$3.872\pm 0.010$	$0.448\pm 0.001$	$0.951\pm 0.001$	$0.050\pm 0.001$
Temporal quadratic	Concentrated	$1.054\pm 0.001$	$1.055\pm 0.001$	$2.889\pm 0.024$	$2.886\pm 0.013$	$0.637\pm 0.002$	$0.949\pm 0.001$	$0.121\pm 0.001$
Temporal quadratic	Diffuse	$1.055\pm 0.001$	$1.057\pm 0.001$	$2.336\pm 0.032$	$2.373\pm 0.025$	$0.633\pm 0.005$	$0.942\pm 0.002$	$0.049\pm 0.001$

Appendix

Table 7: Width, coverage, and residual-normality results for Experiment 3 at

d=24

. Entries are means

\pm

standard errors over 150 replications. Trace denotes the finite-trace approximation.

Figure 8: Temporal quadratic-window bridge for a stationary Gaussian AR(1) process with $\rho=0.8$ and base dimension $d=20$ . Panel (a) compares the empirical split-conformal width with its finite-trace approximation. Panel (b) reports the KS distance between the conformal-width distribution and its independent-window benchmark, together with the residual-normality KS diagnostic. Panel (c) reports split-conformal coverage, with the dashed line denoting the nominal level $0.95$ . Error bars, where shown, are 95% confidence intervals for cell means over 150 replications.

Figure 9: Absolute-regularity corrections for the temporal quadratic-window experiment with $\rho=0.8$ and $d=20$ . The shaded region contains origin separations $\Delta<20$ , for which consecutive lag windows overlap and the disjoint-window transfer bound is not applied. Panel (a) shows the untruncated corrections $(N_{W}-1)\beta(h)$ and $(N_{C}-1)\beta(h)$ ; the horizontal dotted line marks one, above which the corresponding probability correction is vacuous. Panel (b) shows the resulting truncated upper bounds for the width-distribution discrepancy and split-conformal miscoverage. The horizontal dotted line in panel (b) marks the nominal miscoverage level $\alpha=0.05$ . Vertical dotted lines identify the four forecast-origin separations used in the simulation.

Figure 10: Sensitivity to unequal total data budgets. Split conformal uses $N_{C}=400$ , while Bayesian ridge varies over $N_{B}/N_{C}\in\{0.5,0.75,1,1.5,2\}$ . The vertical dotted line marks equal budgets, and the dashed line in the high-dimensional row marks $N_{B}=p=300$ . Results are means over 200 paired replications with 95% confidence intervals.

Dataset	Observations	$H=1$	$H=3$	$H=7$
LA ozone	16,071	6420/6420/3210	6419/6419/3210	6417/6417/3210
Solar	52,560	21009/21009/10506	21008/21008/10506	21007/21007/10504
Beijing PM 10	35,064	14016/14016/7008	14015/14015/7008	14013/14013/7008
Exchange returns	7,587	3029/3029/1515	3028/3028/1515	3026/3026/1515

Appendix

Table 8: Processed series lengths and chronological fitting/calibration/test sample sizes. Each triplet is reported as

n_{\rm fit}/n_{\rm cal}/n_{\rm test}

Dataset and horizon	$\lambda$	Adaptive $\eta$	TWCP $\rho$	ResCP $T$	ResCP $W$
Beijing PM 10 , $H=1$	0.003	0.05	0.95	0.10	all
Beijing PM 10 , $H=3$	0.001	0.02	0.95	0.15	all
Beijing PM 10 , $H=7$	0.001	0.01	0.99	0.15	all
Exchange returns, $H=1$	10	0.01	0.95	0.15	1500
Exchange returns, $H=3$	10	0.02	0.95	0.25	100
Exchange returns, $H=7$	10	0.01	0.95	0.15	1000

Appendix

Table 9: Hyperparameters selected for the primary fixed-origin experiments at 95% nominal coverage. Here

T

is the ResCP sampling temperature and

W

is its residual-history window.

Dataset	Method	RMSE	Coverage	Width	Winkler
LA ozone	Bayesian ridge	$0.383\pm 0.003$	$0.981\pm 0.006$	$1.927\pm 0.123$	$2.054\pm 0.087$
	SCP		$0.959\pm 0.008$	$1.639\pm 0.057$	$1.929\pm 0.008$
	Asymmetric SCP		$0.958\pm 0.006$	$1.630\pm 0.047$	$1.927\pm 0.011$
	Adaptive CP		$0.953\pm 0.001$	$1.587\pm 0.005$	$1.921\pm 0.011$
	Time-weighted CP		$0.953\pm 0.001$	$1.586\pm 0.002$	$1.921\pm 0.011$
	ResCP		$0.947\pm 0.005$	$1.542\pm 0.032$	$1.936\pm 0.010$

Appendix

Table 10: Repeated chronological one-step results at 95% nominal coverage. Entries are means

\pm

standard deviations over three shifted forecast periods. RMSE is reported once per dataset because all UQ methods use the same NGRC point predictor within each origin. Because the windows overlap, the standard deviations quantify forecast-period sensitivity rather than sampling uncertainty or standard errors.

$H$	Method	Cov.	Width	Winkler
1	Bayesian ridge	0.993	2.047	2.083
	SCP	0.963	1.573	1.795
	Adaptive CP	0.954	1.485	1.784
	Time-weighted CP	0.954	1.477	1.782
	ResCP	0.929	1.369	1.797
3	Bayesian ridge	0.997	2.669	2.691

Appendix

Table 11: Complete availability-safe multi-horizon NGRC uncertainty-quantification results at nominal 95% coverage.

Figure 11: Fixed-origin empirical versus nominal coverage for all six NGRC UQ procedures. The dashed diagonal denotes exact marginal calibration.

			Bayesian ridge			SCP
Dataset	Features	RMSE	Coverage	Width	Winkler	Coverage	Width	Winkler
LA ozone	Linear lags	0.357	0.993	2.098	2.132	0.964	1.583	1.806
	Quadratic NGRC	0.360	0.993	2.047	2.083	0.963	1.573	1.795
	Short-memory NGRC	0.360	0.993	2.047	2.082	0.963	1.574	1.796
	Medium-memory NGRC	0.360	0.993	2.047	2.082	0.963	1.572	1.795
	Long-memory NGRC	0.355	0.992	2.040	2.079	0.963	1.555	1.783

Appendix

Table 12: Feature-construction robustness at

H=1

and 95% nominal coverage. RMSE is reported once per feature construction because Bayesian ridge and SCP use the same fitted point predictor.

Dataset	Horizon	Coverage	Coverage range	Width	Winkler score
LA ozone	$H=1$	$0.933\pm 0.003$	$[0.929,0.936]$	$1.375\pm 0.006$	$1.818\pm 0.017$
	$H=3$	$0.954\pm 0.001$	$[0.953,0.955]$	$1.755\pm 0.012$	$2.095\pm 0.006$
	$H=7$	$0.952\pm 0.001$	$[0.950,0.953]$	$1.772\pm 0.008$	$2.111\pm 0.012$
Solar	$H=1$	$0.960\pm 0.003$	$[0.955,0.963]$	$0.059\pm 0.001$	$0.078\pm 0.002$
	$H=3$	$0.928\pm 0.010$	$[0.915,0.937]$	$0.183\pm 0.012$	$0.258\pm 0.009$
	$H=7$	$0.928\pm 0.009$	$[0.919,0.936]$	$0.419\pm 0.008$	$0.576\pm 0.010$

Appendix

Table 13: Sensitivity of ResCP to its stochastic seed. Entries are means

\pm

standard deviations over five seeds; the coverage range gives the minimum and maximum empirical coverage across those seeds. Standard deviations describe computational seed sensitivity, not sampling standard errors.

Dataset	Method	$H=3$	$H=7$
LA ozone	Bayesian ridge	0.0125	0.0107
	SCP	0.0003	0.0007
	Adaptive CP	6.4812	6.4685
	Time-weighted CP	9.5753	9.2434
	ResCP	22.2484	31.4085
Solar	Bayesian ridge	0.1008	0.0988

Appendix

Table 14: Uncertainty-quantification runtime in seconds for forecast horizons

H=3

and

H=7

. Times include pre-test uncertainty-quantification hyperparameter selection where applicable and exclude fitting of the common NGRC point predictor.

Explore similar work

Inference of Unknown Dynamical Components Using Next Generation Reservoir Computing: From Chaotic Systems to Climate Data
Sep 21, 2026Jule Budnick, Andrew Keane, Serhiy YanchukNonlinear System Identification Reservoir Computing
School of Mathematical Sciences, University College Cork, Cork, Ireland · Potsdam Institute for Climate Impact Research, Potsdam, Germany
RareCP: Regime-Aware Retrieval for Efficient Conformal Prediction
May 9, 2026Manuel Heurich, Maximilian Granz, Tim LandgrafWeighted Conformal Prediction Time Series Forecasting
Institute for Computer Science Freie Universität Berlin Arnimallee 7 14195 Berlin
Retrieval-Corrected Conformal Prediction for Time Series
Aug 11, 2026Sangjin Jin, Kangmin Kim, Junhyeong Lee +1Time Series Forecasting Conformal Prediction
Ulsan National Institute of Science and Technology Ulsan, Republic of Korea

Sep 21, 2026J/K move · Enter open · S save

Inference of Unknown Dynamical Components Using Next Generation Reservoir Computing: From Chaotic Systems to Climate Data

Jule Budnick, Andrew Keane, Serhiy Yanchuk

School of Mathematical Sciences, University College Cork, Cork, Ireland · Potsdam Institute for Climate Impact Research, Potsdam, Germany

Nonlinear System Identification Reservoir Computing

Open paper AskPDF arXiv

We investigate next generation reservoir computing (NGRC) as a data-driven approach for inferring unseen components of dynamical systems. We compare NGRC with traditional reservoir computing (RC) using the Lorenz and Rössler system, where two unknown components are inferred from one given component. For both systems, NGRC achieves accurate results while requiring fewer training data and less computational time than RC. We identified an inverse proportional behavior between the number of time-delayed steps needed for NGRC and the temporal resolution, indicating that the physical time span covered by the delay interval is an important factor in determining the required number of delayed steps. Finally, we apply NGRC to the observational climate data of ENSO (El Niño--Southern Oscillation) and infer one observable from the remaining variables. Despite the noise and complexity of the real-world data, the NGRC shows promising results. Our findings demonstrate the potential of NGRC for efficient inference of unseen components in both controlled dynamical systems and real-world data.

Abstract

Explore similar work

Inference of Unknown Dynamical Components Using Next Generation Reservoir Computing: From Chaotic Systems to Climate Data

RareCP: Regime-Aware Retrieval for Efficient Conformal Prediction

Retrieval-Corrected Conformal Prediction for Time Series