conformal prediction, uncertainty quantification, bootstrap, time series, distribution-free inference

Introduction and Statement of Need

Practitioners pair point forecasts with calibrated intervals, but the two distribution-free tools that would supply one assume away temporal structure: split conformal prediction assumes exchangeable calibration data, the ordinary bootstrap assumes IID observations, and both undercover on autoregressive streams. The time series bootstrap restores validity under dependence, and conformal calibration supplies the finite-sample coverage guarantee. Prior libraries supply one piece: forecasting toolkits (skforecast, darts, statsforecast, neuralforecast) attach residual-bootstrap or split-conformal intervals to their own forecasters, conformal libraries (MAPIE, TorchCP, puncc) calibrate a fitted predictor, and resampling libraries (arch) stop at classical intervals. To our knowledge, no existing library combines a general dependence-aware resampling engine (block, residual, sieve, and wild methods with typed specs, for arbitrary statistics) with an adaptive conformal layer and a streaming reduce in one API; tsbootstrap provides this combination. Derivations appear in the companion methods manuscript (arXiv:2404.15227).

Design and API

The single typed entry point bootstrap(X, *, method=spec, ...) selects the method by a concrete specification object whose type identifies it, visible to editors and type checkers. A typical analysis consists of three calls: diagnose(X) to recommend a spec, bootstrap(X, method=spec, n_bootstraps=999, random_state=0), and conf_int(res, statistic="mean").

The call returns a structured BootstrapResult: samples (the replicate array), provenance (spec, seed, backend), out-of-bag (the held-out index EnbPI calibrates on), and in-bag (the drawn index, for jackknife-after-bootstrap and studentized errors). Inputs are handled by a narwhals layer (a lightweight compatibility layer over dataframe libraries): NumPy arrays, Python lists, and frames or series from pandas, Polars, and PyArrow. Each spec carries a MethodMetadata record (metadata_for(spec)) declaring assumptions, capabilities, references, cost, and failure modes; diagnose(X) reads the lag-one autocorrelation and an ADF stationarity test and returns a frozen Diagnosis naming recommended specs with an auto-selected block length.

Results are deterministic per backend for a fixed seed: each replicate draws from its own RNG stream, so output is invariant to replicate ordering, parallelism, and thread count. The default numpy backend spawns per-replicate PCG64 streams via SeedSequence and is bit-for-bit reproducible for a pinned seed and environment. The optional compiled numba backend ([accel] extra), the source of the Section [sec:sec:perf] throughput, draws from a distinct counter-based stream, equal in distribution but not bit-identical to the numpy path; its byte stream is pinned within a release, not across releases.

Methods and UQ Coverage

Block methods cover moving [1], circular [2], stationary [3], non-overlapping [4], and tapered blocks with energy-normalized windows [5], with block length defaulting to the automatic spectral-density plug-in rule of [6], with the correction of [7]. Model-based bootstraps (ResidualBootstrap over AR, ARIMA, and VAR) regenerate series recursively from the fitted dynamics and resampled centered innovations [8], rather than the less accurate fitted-plus-residuals scheme; SieveAR implements the sieve bootstrap [9], [10], with exogenous support (ARX, VARX, ARIMAX). Innovations may be resampled IID (default), wild (Rademacher by default) [11]–[14], or block-wild for dependent innovations [15].

Classical intervals via conf_int cover percentile [16], basic, studentized bootstrap-t with a dependence-aware block-jackknife standard error [1], and BCa [17], the last gated to the IID spec because the studentized route is the second-order-correct choice under dependence [18]. The uq API adds conformal and adaptive intervals: EnbPI out-of-bag intervals for sklearn-style regressors [19], AR forecast bands, and static, sliding-window, ACI [20], NexCP [21], and AgACI [22] calibrators for stationary, volatility-clustered, and drifting regimes. Table 1 reports the empirical coverage these methods achieve.

Coverage under dependence. The study crosses five DGPs, 5,000 Monte Carlo datasets each, eight methods, and four functionals: 800,000 percentile intervals from \(B{=}999\) replicates at \(\alpha{=}0.10\); the table shows the mean for four DGPs and four methods. DGPs: white noise \(N(0,1)\); AR(1), \(\phi{=}0.9\), Gaussian innovations; an MA(2) control, \(\theta{=}(0.6, 0.3)\) (not displayed); AR(1)+ARCH(1), \(x_t = 0.5x_{t-1} + \sigma_t\varepsilon_t\) with \(\sigma_t^2 = 0.2 + 0.4x_{t-1}^2\); ARFIMA(\(0, 0.4, 0\)); series lengths in the caption. The ARFIMA series are generated exactly by the Davies-Harte circulant-embedding method; the white-noise series are drawn directly, the MA(2) series by a direct causal filter, and the AR(1) and AR+ARCH series recursively from their defining equations, with burn-ins of 400 (AR(1)), 200 (MA(2)), and 300 (AR+ARCH). Coverage targets are exact-analytic from each DGP’s closed-form autocovariance (for AR+ARCH via its martingale-difference property), corroborated by independent Monte Carlo to within 1%. Block lengths use the automatic rule above; the sieve selects its order by BIC; the design is paired (one shared bootstrap seed per dataset across methods); the maximum Monte Carlo standard error is 0.7 percentage points. Dependence-aware methods correct much of the IID bootstrap’s undercoverage: on AR(1) with \(\phi{=}0.9\) at \(n{=}200\), coverage recovers from 27.8% to roughly 70% (blocks) and 83.1% (sieve), nearest to nominal under short-memory linear dependence; the blocks fall short because the dependence is strong relative to the series length. Under AR(1)+ARCH(1) the sieve reaches 89.1%, within a point of nominal. No method considered attains acceptable coverage under long-memory ARFIMA (\(d{=}0.4\)), so the library exposes the method choice rather than imposing a single default.

Comparison

Table 2 compares the five nearest libraries; tsbootstrap is the only column with entries in both the resampling and conformal rows. TorchCP [23] and MAPIE target post-hoc conformal calibration, tsbootstrap dependence-aware uncertainty quantification. Because the adaptive calibrators are pure functions over user-supplied nonconformity scores, residuals from any prediction pipeline can be calibrated. Beyond the table, puncc ships EnbPI; statsforecast and neuralforecast attach split-conformal intervals to their own forecasters; none has a general resampling engine.

3.5pt

Performance

Our harness (bench_vs_arch.py) runs the compiled reducer and arch’s bs.apply(np.mean, reps=B), each library’s own streaming reducer and the closest like-for-like path, on the same warmed workload on a dedicated Hetzner ccx33 (eight vCPU, AMD EPYC-Milan; Ubuntu 24.04, Python 3.12.3, NumPy 2.4.6, arch 8.0.0, numba 0.65.1), measured at tsbootstrap 0.4.0 (commit 53aa089) and confirmed directionally by a pinned v0.6.1 re-run. The comparison is not symmetric in generality: arch pays a per-replicate Python callback for an arbitrary statistic, while the compiled kernel fixes the statistic, fuses the loop, and runs replicate-parallel on all eight threads against arch’s single-threaded apply loop; the speedups below quantify this fixed-statistic path only.

The sixteen-cell grid crosses the four methods both libraries implement (IID, MovingBlock, CircularBlock, StationaryBlock) with \(n \in \{200, 2000\}\) and \(B \in \{999, 10000\}\) on an AR(1) input (\(\phi{=}0.6\)), block length 20, mean statistic; each cell is the median of three warmed repetitions, JIT and import excluded (full grid in benchmarks/). The compiled backend is faster in all sixteen cells: at \(n{=}2000\) it runs \(3.9\times\) faster on IID at \(B{=}999\) (3.0 versus 11.7 ms) and \(34.0\times\) faster on CircularBlock at \(B{=}10000\) (4.9 versus 167.0 ms); the \(n{=}200\) cells are larger still, from arch’s per-replicate Python overhead; pinned to one thread it still leads by roughly \(2\times\) to \(8\times\) (v0.6.1 re-measurement). For a generality-matched comparison, the default numpy backend (which, like arch, accepts an arbitrary statistic) runs \(1.3\times\) to \(2.6\times\) slower than arch on the same cells; the compiled speedup comes entirely from fusing a fixed statistic into the kernel. On memory, bootstrap_reduce avoids materializing the \(O(Bn)\) replicate tensor (the numpy backend reduces in 2,048-replicate chunks; the compiled kernel holds only the \((B, |\theta|)\) statistic array). At \(n{=}2000\), \(B{=}50000\) (MovingBlock mean), peak resident memory above the interpreter floor is 1,944 MB materializing every replicate versus 20 MB streaming, about \(96\times\) (\(103\times\) at \(B{=}10000\); fresh subprocesses, warmed JIT, benchmarks/results/). The streaming footprint grows only with \(B\), about 0.4 KB per replicate, not with \(n \times B\).

Reproducibility, Quality, and Availability

tsbootstrap is MIT licensed (github.com/astrogilda/tsbootstrap, PyPI v0.6.1, Zenodo DOI 10.5281/zenodo.8226495) and supports Python 3.10–3.13 with five core dependencies (numpy, scipy, pydantic, scikit-base, narwhals); statsmodels, scikit-learn, and numba are optional extras. Test coverage is 93% overall with a blocking 80% per-file CI gate; GitHub Actions CI is green on main across a 12-cell matrix (Python 3.10–3.13 on Linux, macOS, and Windows) and further gates lint, zero-error mypy and pyright, executed notebooks, and the docs build; a nightly mutation ratchet (about 1,100 mutants over the engine and model core) fails on any new survivor. The core also ships as sktime/skbase estimator adapters (validated by check_estimator) plus a read-only MCP server (tsbootstrap-mcp, stdio, uvx-runnable) for LLM agents; 14 CI-executed notebooks are at tsbootstrap.readthedocs.io.

References↩︎