Paper deep dive

SPACE: Sample-cloud Predictive Adaptive Conformal Ellipsoids for Multivariate Time-Series Forecasting

Baishi Li, Kelvin J. L. Koa, Ke-Wei Huang

Year: 2026Venue: arXiv preprintArea: stat.MLType: PreprintEmbeddings: 122

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%

Last extracted: 8/19/2026, 4:47:35 AM

Summary

The paper introduces SPACE (Sample-cloud Predictive Adaptive Conformal Ellipsoids), a conformal wrapper for multivariate time-series forecasting. It addresses limitations in existing methods by decoupling geometry estimation from calibration radius adjustment. SPACE estimates local covariance geometry directly from the current forecast sample cloud and selects an adaptive historical calibration window based on regime consistency, achieving superior coverage-efficiency trade-offs under distribution shifts.

Entities (8)

SPACE · method · 100%Asian Institute of Digital Finance · organization · 95%Conformal Prediction · technique · 95%Multivariate Time Series Forecasting · task · 95%National University of Singapore · organization · 95%Adaptive Conformal Inference · baseline-method · 90%Ellipsoidal Joint Prediction Regions · output-structure · 90%MultiDimSPCI · baseline-method · 90%

Relation Signals (7)

Kelvin J.L. Koa → affiliatedwith → National University of Singapore

confidence 95% · Kelvin J.L. Koa ... Affiliation: Asian Institute of Digital Finance, National University of Singapore

Ke-Wei Huang → affiliatedwith → National University of Singapore

confidence 95% · Ke-Wei Huang ... Affiliation: Asian Institute of Digital Finance, National University of Singapore

Baishi Li → affiliatedwith → National University of Singapore

confidence 95% · Baishi Li ... Affiliation: Asian Institute of Digital Finance, National University of Singapore

SPACE → isa → Conformal Wrapper

confidence 95% · we propose SPACE, a conformal wrapper for sample-generating multivariate forecasters.

SPACE → uses → Sample-cloud Predictive Adaptive Conformal Ellipsoids

confidence 95% · SPACE constructs ellipsoidal joint prediction regions by estimating time-local covariance geometry directly from the current forecast sample cloud

SPACE → addresseslimitationof → MultiDimSPCI

confidence 90% · MultiDimSPCI [40] is the closest benchmark to our work. It constructs sequential ellipsoidal regions by using historical residual covariance... Space differs by extracting the ellipsoidal geometry from the current predictive sample cloud

SPACE → outperforms → MultiDimSPCI

confidence 90% · reduces the overall mean coverage gap by over 80% relative to the strongest baseline (MultiDimSPCI)

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Modern probabilistic time-series forecasters often express uncertainty through forecast samples. While typically converted into nominal prediction regions using empirical quantiles, these model-implied sets lack formal coverage guarantees and frequently deviate from nominal targets under distribution shift. Existing multivariate conformal methods can calibrate these regions online, but they typically estimate geometry from historical residuals using fixed or accumulating look-back windows. This reliance on the past limits their ability to exploit the instantaneous dependence structure of current predictions and leaves them vulnerable to stale-regime contamination. To address this, we propose SPACE, a conformal wrapper for sample-generating multivariate forecasters. SPACE constructs ellipsoidal joint prediction regions by estimating time-local covariance geometry directly from the current forecast sample cloud, calibrating the region's radius via a dynamic backward window-selection scheme. Across diverse multivariate datasets, probabilistic forecasters, and conformal baselines, SPACE consistently brings realized joint and rolling coverage closer to the nominal target, achieving superior coverage-efficiency tradeoffs relative to competing wrappers.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

118,758 characters extracted from source content.

Expand or collapse full text

SPACE: Sample-cloud Predictive Adaptive Conformal Ellipsoids for Multivariate Time-Series Forecasting Baishi Li Kelvin J.L. Koa Affiliation: Asian Institute of Digital Finance, National University of Singapore[2pt] baishili@u.nus.edu, kelvin.koa@u.nus.edu, dishkw@nus.edu.sg Ke-Wei Huang Affiliation: Asian Institute of Digital Finance, National University of Singapore[2pt] baishili@u.nus.edu, kelvin.koa@u.nus.edu, dishkw@nus.edu.sg [4pt] Department of Information Systems Analytics National University of Singapore Abstract Modern probabilistic time-series forecasters often express uncertainty through forecast samples. While typically converted into nominal prediction regions using empirical quantiles, these model-implied sets lack formal coverage guarantees and frequently deviate from nominal targets under distribution shift. Existing multivariate conformal methods can calibrate these regions online, but they typically estimate geometry from historical residuals using fixed or accumulating look-back windows. This reliance on the past limits their ability to exploit the instantaneous dependence structure of current predictions and leaves them vulnerable to stale-regime contamination. To address this, we propose Space, a conformal wrapper for sample-generating multivariate forecasters. Space constructs ellipsoidal joint prediction regions by estimating time-local covariance geometry directly from the current forecast sample cloud, calibrating the region’s radius via a dynamic backward window-selection scheme. Across diverse multivariate datasets, probabilistic forecasters, and conformal baselines, Space consistently brings realized joint and rolling coverage closer to the nominal target, achieving superior coverage-efficiency tradeoffs relative to competing wrappers. 1 Introduction While modern deep forecasting architectures—spanning diffusion [27, 19, 43], autoregressive [31, 26], copula [3, 30], and flow-based models [28, 14]—excel at generating multivariate sample clouds, converting these rich representations into statistically calibrated joint prediction regions under non-stationarity remains an open challenge. Building on foundational conformal prediction work [36, 18, 2], recent methods have made major strides in sequential coverage calibration [41, 11, 1], often by using historical residual covariance to shape multivariate prediction sets [40, 24]. This approach is natural and effective, but it can adapt slowly when the instantaneous uncertainty geometry changes abruptly. In this paper, we build on these foundations by showing that robustness to sudden shifts can be improved by separating two roles that are often coupled: geometry estimation and conformal radius calibration. Our method shapes joint regions using the forecaster’s current predictive sample covariance, while calibrating the radius using an adaptively selected, same-regime historical window. The practical consequences of delayed geometric adaptation become visible under non-stationary dynamics. As illustrated in Figure 1(a), raw sample-based regions from modern forecasters can be substantially miscalibrated at the joint level across diverse datasets, often deviating due to model misspecification and finite-sample error [29, 20]. Furthermore, as shown in Figure 1(b), this realized coverage often drifts away from the nominal 90%90\% target over time. Many multivariate conformal baselines correct this coverage gap by relying on past residuals to estimate both region geometry and calibration thresholds. After abrupt regime changes, this coupling creates a difficult trade-off: expanding the radius may recover coverage but sacrifice efficiency, while maintaining a tight region with outdated geometry preserves efficiency at the cost of local coverage. Even online adaptive methods, such as Adaptive Conformal Inference (ACI) [11, 34], adjust only the scalar miscoverage rate, leaving the geometric shape inherited from historical residuals unchanged. The closest prior work, MultiDimSPCI [40], constructs ellipsoidal regions from historical residual covariance within a fixed-length sequential window. We address this limitation by ensuring that both the region’s geometry and its calibration window are chosen adaptively in response to the current forecast state. (a) (b) Figure 1: Joint coverage at the 90% target. Panel (a) shows the distribution of empirical joint coverage. The x-axis evaluates different constructions: Raw Forecaster denotes native sample-based regions without conformal calibration, Space is our method, and the rest are competing conformal baselines. The y-axis shows empirical joint coverage, with the dotted line marking the 90% target; better methods concentrate near this line. Panel (b) shows rolling joint coverage on the weather dataset across forecasters. Raw rolling coverage (gray) is unstable and systematically off-target, whereas Space (red) stays much closer to the desired level over time. To address these limitations, we introduce Space (Sample-cloud Predictive Adaptive Conformal Ellipsoids), a post-hoc conformal wrapper that directly exploits the probabilistic capabilities of the underlying forecaster. Our primary contributions are summarized as follows: • Decoupling geometry from calibration. We formalize a principle for multivariate conformal forecasting under non-stationarity: the forecaster’s predictive samples should determine the local shape of the joint region, while conformal calibration scores should determine its radius. Based on this principle, we introduce Space, a post-hoc wrapper that constructs ellipsoidal regions using the time-local covariance of the current predictive sample cloud, reducing the adaptation delay associated with residual-based conformal geometries. • Regime-local adaptive calibration. We propose a backward same-regime testing procedure for calibrating the radius of the forecast-conditioned ellipsoid. By adaptively selecting the longest recent calibration window whose conformity scores remain statistically compatible with the current regime, the method reduces the influence of outdated dynamics while balancing local coverage and statistical efficiency in sequential prediction. • Theoretical guarantees and error decomposition. We establish finite-sample joint coverage under an idealized same-regime assumption and derive a sequential local coverage-gap bound that separates the effects of sample-based geometry estimation, adaptive window selection, and finite calibration sample size. This analysis clarifies how predictive sample size and calibration-window quality influence local coverage error. • Extensive empirical validation. We evaluate Space across 7 multivariate datasets spanning energy, finance, weather, and traffic domains, 8 diverse probabilistic forecasters (including diffusion, flow, and autoregressive architectures), and 5 competing conformal baselines (MultiDimSPCI, local ellipsoidal CP, and copula-based CP). Across these diverse benchmark pairs, Space achieves a median realized joint coverage gap of just 0.2 percentage points from the nominal target. Furthermore, it reduces the overall mean coverage gap by over 80% relative to the strongest baseline (MultiDimSPCI) while consistently improving local calibration stability under distribution shift. 2 Related Work Sample-generating probabilistic forecasters. Modern forecasters—including autoregressive [31], diffusion [27, 33, 19, 43], copula [30, 3], and flow-based models [14]—encode joint dependencies through predictive sample clouds. However, these native sample regions remain model-implied rather than coverage-calibrated. We treat these models as black boxes, converting their raw predictive samples into formally calibrated joint prediction regions without altering the underlying architecture. Sequential, localized, and online conformal prediction. A growing literature studies conformal prediction beyond the exchangeable setting. Methods such as EnbPI [41], Adaptive Conformal Inference (ACI) [11], and Conformal PID [1] robustly adapt calibration thresholds or miscoverage levels online. Recent extensions have improved online tracking [6], rolling risk control [1], covariate shift and localized conformal prediction [34, 4, 12]. However, these methods do not specify a forecast-conditioned geometry for conformalized multivariate joint regions. Multivariate and ellipsoidal conformal prediction. Prior work has developed conformal prediction for cross-sectional multivariate regression using structured dependence via copulas [23, 32] and ellipsoidal sets [24, 13, 10]. For multivariate time series, MultiDimSPCI [40] is the closest benchmark to our work. It constructs sequential ellipsoidal regions by using historical residual covariance and handles temporal dependence through sequential calibration within a fixed-length window. Space differs by extracting the ellipsoidal geometry from the current predictive sample cloud—rather than from a historical surrogate—and by adaptively selecting the calibration window via same-regime testing. Another approach, FCP [17], targets multivariate time series by learning a context-conditioned residual flow, whereas Space requires no auxiliary generative model fitting. Positioning of Space. Space is a conformal wrapper that simultaneously (i) operates on sample-generating multivariate forecasters without auxiliary model fitting, (i) derives its multivariate geometry from the current predictive sample cloud rather than historical residuals, and (i) adapts its calibration window to the current regime. Table 1 summarizes these structural differences. Table 1: Positioning of Space relative to representative conformal approaches for multivariate and sequential prediction. Method Sequential TS Joint multivariate set Geometry source Requires auxiliary CP model Adaptive calibration window Adaptive TS CP [41, 11, 1] ✓ ✗ Scalar / marginal residuals ✗ ✗∗ MultiDimSPCI [40] ✓ ✓ Historical residual covariance ✗ ✗ Copula-based CP [23, 32] Mixed ✓ Fitted copula dependence ✓ ✗ FCP [17] ✓ ✓ Learned residual flow ✓ ✗ Space (ours) ✓ ✓ Current predictive samples ✗ ✓ ∗These methods perform adaptive online scalar updates, but do not adaptively select contiguous historical calibration windows. 3 Method: Space We consider a probabilistic multivariate forecaster that outputs a set of M predictive samples (e.g., generated trajectories from a conditional generative model) for a target yt∈ℝdy_t ^d, where d spans multiple variables, forecast horizons, or both. We denote this predictive sample cloud by ^t=y^t(m)m=1M,y^t(m)∈ℝd. Y_t=\ y_t^(m)\_m=1^M, y_t^(m) ^d. (1) Given a target miscoverage level α∈(0,1)α∈(0,1), our goal is to construct a joint prediction set t(α)⊆ℝdsuch thatPr⁡(yt∈t(α))≈1−α,C_t(α) ^d that \! (y_t _t(α) )≈ 1-α, (2) while maximizing statistical efficiency. Operating strictly causally over a chronological stream (^t,yt)t=1T\( Y_t,y_t)\_t=1^T, Space constructs this set sequentially. At prediction time t, it uses only the current sample cloud ^t Y_t and realized past data (s<ts<t). The method consists of two core components: First, Space extracts a time-local multivariate geometry directly from ^t Y_t to shape the region (Section 3.1). Second, it calibrates the ellipsoid’s radius via a sequential same-regime procedure that dynamically selects the optimal historical calibration window (Section 3.2). 3.1 Sample-derived local geometry For each time step t, we first construct a median-based center μt _t from ^t=y^t(m)m=1M Y_t=\ y_t^(m)\_m=1^M, and form predictive residual samples rt(m)=y^t(m)−μtr_t^(m)= y_t^(m)- _t. To account for varying marginal scales, we compute the empirical covariance of the standardized residuals: C^t=Cov⁡(r~t(1),…,r~t(M)), C_t=Cov\! ( r_t^(1),…, r_t^(M) ), (3) where r~t(m)=D−1rt(m) r_t^(m)=D^-1r_t^(m) and D=diag⁡(s1,…,sd)D=diag(s_1,…,s_d) with sj>0s_j>0 being the coordinate-wise scale estimate fitted from the wrapper-training segment. Mapping this covariance back to the original scale gives the local multivariate geometry at time t: Σ^t=DC^tD. _t=D C_tD. (4) Crucially, as we establish in Proposition 1, this sample-derived covariance acts as a statistically consistent proxy that robustly concentrates around the forecaster’s true conditional geometry. We define the radial Mahalanobis-type nonconformity score at time t as et=(yt−μt)⊤Σ^t+(yt−μt).e_t= (y_t- _t) _t^+(y_t- _t). (5) Here, Σ^t+ _t^+ denotes the stabilized pseudo-inverse of the local covariance matrix Σ^t _t. Since Σ^t _t is symmetric positive semidefinite, it admits a spectral decomposition Σ^t=UtΛtUt⊤ _t=U_t _tU_t , where Λt=diag⁡(λt,1,…,λt,d) _t=diag( _t,1,…, _t,d) contains the nonnegative eigenvalues. The stabilized pseudo-inverse is formed by inverting only eigenvalues above an adaptive tolerance ρt _t: Σ^t+=Utdiag⁡(λt,i>ρtλt,i)Ut⊤. _t^+=U_tdiag\! ( 1\ _t,i> _t\ _t,i )U_t . (6) The score in (5) therefore measures how far the realized target lies from the center of the forecast cloud relative to its implied local covariance geometry. 3.2 Dynamic same-regime calibration window selection At forecast time t, let t=t−p,…,t−1P_t=\t-p,…,t-1\ denote the recent probe block of length p. Space searches backward over candidate calibration lengths L∈Lmin,Lmin+h,Lmin+2h,…,Lmax,L∈\L_ ,\,L_ +h,\,L_ +2h,…,L_ \, (7) where h is the backward extension step size. For each candidate length L, define the full candidate block twhole(L)=t−p−L,…,t−p−1C_t^whole(L)=\t-p-L,…,t-p-1\ and the newly added increment ℬtinc(L)=t−p−L,…,t−p−L+h−1B_t^inc(L)=\t-p-L,…,t-p-L+h-1\. While the radial score ete_t captures overall ellipsoidal distance, it may miss regime shifts that manifest purely through coordinate-wise extremes or directional changes. Therefore, for the conformal uniformity test, we deploy a bundle of three complementary scalar diagnostics. For any observed time s, let r~s=D−1(ys−μs) r_s=D^-1(y_s- _s) denote the realized standardized residual: s,1=es,s,2=‖r~s‖∞,s,3=|v1⊤r~s|,s_s,1=e_s, s_s,2=\| r_s\|_∞, s_s,3=|v_1 r_s|, (8) where v1v_1 is the top eigenvector of the training residual covariance in standardized space. Thus, s,1s_s,1 captures radial size, s,2s_s,2 captures extremeness, and s,3s_s,3 captures variation along a dominant direction. Conformal same-regime uniformity test. For each diagnostic j∈1,2,3j∈\1,2,3\ and candidate block ∈twhole(L),ℬtinc(L)A∈\C_t^whole(L),\,B_t^inc(L)\, Space compares the probe scores su,j:u∈t\s_u,j:u _t\ with the calibration scores sv,j:v∈\s_v,j:v \ using randomized conformal p-values: Uu,j()=1+#⁡v∈:sv,j>su,j+ξu,j#v∈:sv,j=su,j||+1,ξu,j∼Unif(0,1).U_u,j(A)= 1+\#\v :s_v,j>s_u,j\+ _u,j\,\#\v :s_v,j=s_u,j\|A|+1, _u,j (0,1). (9) The finite-sample validity of this test is formalized in Assumption (A2): clean same-regime blocks have controlled false-rejection probability, while contaminated extensions have controlled false-acceptance probability. Proposition 2 shows that these local test guarantees imply recovery of the oracle same-regime window with high probability. Accordingly, A passes the test if, for each diagnostic j: KSt,L,j=supx∈[0,1]|F^t,L,j(x)−x|<τt,L,j,j=1,2,3,KS_t,L,j^\,A= _x∈[0,1] | F_t,L,j^\,A(x)-x |< _t,L,j^\,A, j=1,2,3, (10) where F^t,L,j(x) F_t,L,j^\,A(x) is the empirical CDF of the probe p-values, KSt,L,jKS_t,L,j^\,A is the Kolmogorov–Smirnov discrepancy, and τt,L,j _t,L,j^\,A is the familywise-corrected Dvoretzky–Kiefer–Wolfowitz (DKW) threshold. Sequential extension rule. Starting from LminL_ , Space extends the calibration window from L to L+hL+h only if both twhole(L)C_t^whole(L) and ℬtinc(L)B_t^inc(L) pass (10). This dual-check guards against dilution, since a stale older increment could be masked by the previously accepted history if only the full block were tested. The selected calibration length LtL_t is the last accepted length before the first rejection. 3.3 Prediction Set Given LtL_t, Space calibrates the radius using the same-regime score history and the probe block. Let ℋt(Lt)=es:s∈twhole(Lt)∪tH_t(L_t)=\e_s:s _t^whole(L_t) _t\ (11) denote the final calibration score set available at time t. The radius is defined as rt=Q1−αt(ℋt(Lt)),r_t=Q_1- _t\! (H_t(L_t) ), (12) where Q1−αt(⋅)Q_1- _t(·) is the empirical quantile at level 1−αt1- _t, and αt _t is updated online using adaptive conformal inference (ACI) [11]. The final prediction set at time t is the ellipsoid t(α)=y∈ℝd:(y−μt)⊤Σ^t+(y−μt)≤rt2.C_t(α)= \y ^d:(y- _t) _t^+(y- _t)≤ r_t^2 \. (13) Figure 2 visualizes the local prediction geometry used by our method. The ellipses are not fixed in either orientation or aspect ratio: instead, they rotate, elongate, or contract over time to match the local structure of the forecast sample cloud. This behavior is a direct consequence of our time-local covariance estimated from the contemporaneous generative sample cloud. Figure 2: Illustration of local prediction geometry over time on the Weather--CSDI stream, shown on randomly selected target dimensions (d11d_11 and d12d_12) and time indices. Each ellipse is the 2D projection of Space’s ellipsoidal prediction region, while the gray points denote the forecast sample cloud. By rigorously decoupling the sample-based shape estimation from the same-regime scale calibration, Space transcends empirical heuristics to achieve the formal finite-sample coverage guarantees and theoretical error decomposition presented next in Section 4. Computational Efficiency. Despite performing a dynamic backward search at every time step, Space introduces minimal operational overhead. The backward search relies strictly on one-dimensional empirical CDF comparisons (KS test) over the pre-computed scalar st,js_t,j, thereby bypassing repeated matrix inversions. Algorithm 1 in Appendix B summarizes the full execution. 4 Theoretical Analysis We analyze through: (i) finite-sample concentration of the sample-cloud covariance, (i) recovery of a same-regime calibration window, and most importantly — (i) a sequential coverage-error decomposition. The four sets of formal assumptions and other proofs are deferred to Appendix A. Assumptions (informal). Let ℱtF_t be the filtration generated by realized targets, forecast clouds, and auxiliary randomness up to time t. We define ℙt(⋅)=ℙ(⋅∣ℱt−1)P_t(·)=P(· _t-1) as the probability conditional on the immediate past. For same-regime properties, we use ℚt(⋅)=ℙ(⋅∣ℱτt)Q_t(·)=P(· _ _t), conditioning only on information prior to the most recent change point τt _t. Our per-time coverage bounds are thus conditional; unconditional variants follow by integrating over the past. (A1) Sample-cloud regularity and stable score geometry: conditional on the past, the M forecast samples are i.i.d. sub-Gaussian draws from the forecaster-implied law, their standardized population covariance is well defined, and the resulting Mahalanobis conformal set is locally stable under small perturbations of the sample-cloud covariance [37, 40]. This assumption provides the regularity needed for the LΣL_ term below and rules out discontinuities caused by atoms in the ideal score distribution or eigenvalues crossing the pseudo-inverse threshold.11 1 Assumption A1 is natural for generative architectures whose output samples are obtained by passing light-tailed latent variables through sufficiently regular maps. For example, if a standard Gaussian latent variable is transformed by an L-Lipschitz generator fθf_θ, Gaussian concentration implies sub-Gaussian output behavior up to the Lipschitz scale [35]. Spectral normalization is one common way to encourage such Lipschitz control in neural generators [25]. (A2) Same-regime test calibration and shift detectability: clean candidate blocks are rejected by the same-regime test with probability at most δuni _uni, while every contaminated extension has at least one diagnostic/block comparison whose false-acceptance probability is at most βdet _det [39, 6]. For the KS-style version used by Space, a sufficient DKW-type form is βdet=2exp⁡(−2|t|Δsep2) _det=2 (-2|P_t| _sep^2) under a population separation margin Δsep>0 _sep>0. (A3) Piecewise augmented exchangeability: within each stationary regime, the augmented forecast–target objects t=(Yt,^t)Z_t=(Y_t, Y_t) are conditionally exchangeable, and the ideal conformal scores used for coverage inherit this same-regime exchangeability [4, 11]. (A4) Adaptive-level validity: the adaptive miscoverage level αt _t is causal and does not use the oracle test score. Formally, its effect is controlled by a tracking term εtACI _t^ACI together with a rank-validity condition for the oracle same-regime conformal predictor using αt _t, up to the usual finite calibration-sample error [11, 1]. Step 1: Geometry Concentration. First, we establish that the empirical sample-cloud covariance C^t C_t reliably estimates the forecaster-implied conditional geometry. Proposition 1 (Sample-cloud covariance concentration). Under (A1), for any fixed prediction time t and any δ∈(0,1)δ∈(0,1), there exists a universal constant c0>0c_0>0 such that, conditional on ℱt−1F_t-1, ℙt(‖C^t−Ct⋆‖op≤ηM,δ)≥1−δ,ηM,δ=c0K2(d+log⁡(1/δ)M+d+log⁡(1/δ)M).P_t\! ( \| C_t-C_t \|_op≤ _M,δ )≥ 1-δ, _M,δ=c_0K^2 ( d+ (1/δ)M+ d+ (1/δ)M ). (14) Proof intuition: This follows from the standard operator-norm concentration bound for sample covariance matrices of sub-Gaussian vectors (see [35]). Remark on Statistical Rate and Generative Capacity. The convergence rate ηM,δ _M,δ matches the standard dimension dependence for unstructured sub-Gaussian covariance estimation. Unlike classical settings constrained by historical data, the generative sample size M is fully controllable and can theoretically be arbitrarily large, driving the geometry estimation error to zero. Practically, M need not be large: once M≳dM d, the leading d/M d/M term decays rapidly. Faster estimation rates would require structural assumptions on Ct⋆C_t , such as sparsity or low rank. Step 2: Window Recovery. ρtwin=ℚt(L^t≠Lt⋆) _t^win=Q_t( L_t≠ L_t ) denotes the window-selection error probability, where Lt⋆L_t is the largest candidate length whose calibration and probe blocks lie in the current regime. Proposition 2 (Same-regime window recovery). Under (A2), if the same-regime tests have false-rejection probability at most δuni _uni and false-acceptance probability at most βdet _det, then ρtwin≤6||δuni+||βdet. _t^win≤ 6|G| _uni+|G| _det. (15) In the DKW-style separation case, βdet=2exp⁡−2|t|Δsep2. _det=2 \-2|P_t| _sep^2\.22 2 The Dvoretzky–Kiefer–Wolfowitz inequality, refined by Massart [9, 21], guarantees that the chance of missing a regime shift drops exponentially fast as we collect more recent data (|t||P_t|) or as the shift becomes more severe (Δsep _sep). Proof intuition. Because the window extends sequentially, a false rejection at any clean length prematurely halts the search before reaching the oracle length Lt⋆L_t . A union bound over the 6||6|G| clean tests (2 block types × 3 diagnostics) controls this premature stopping. Conversely, Assumption (A2) ensures that any invalid length L>Lt⋆L>L_t is detected and rejected with probability at least 1−βdet1- _det. A combined union bound over these candidate lengths yields the final recovery probability [39, 6]. Step 3: Coverage Gap Decomposition. Let ℋt⋆=twhole(Lt⋆)∪t,nt⋆=|ℋt⋆|,H_t =C_t^whole(L_t ) _t,n_t =|H_t |, and tACI(α)C_t^ACI(α) denote the final prediction set when the nominal target level is α, allowing the algorithm to use the adaptive internal level αt _t. The per-time coverage error can be decomposed into window-selection error, adaptive-level error, finite calibration-sample error, and sample-geometry perturbation error. Theorem 1. Assume (A1)–(A4). For any fixed prediction time t and any δt∈(0,1) _t∈(0,1), suppose ηM,δt/(nt⋆+1)≤η0 _M, _t/(n_t +1)≤ _0. Then |ℚt(Yt∈tACI(α))−(1−α)|≤ρtwin+εtACI+1nt⋆+1+LΣηM,δt/(nt⋆+1)+δt, |Q_t\! (Y_t _t^ACI(α) )-(1-α) |≤ _t^win+ _t^ACI+ 1n_t +1+L_ _M, _t/(n_t +1)+ _t, (16) where LΣL_ is the score-geometry stability constant from (A1), and ηM,δ _M,δ is the covariance concentration rate from Proposition 1. The window-selection term ρtwin _t^win satisfies Proposition 2. Proof intuition. The proof compares the empirical prediction set with an oracle same-regime set. The coverage gap decomposes into additive penalties for window-selection failure, ACI tracking, finite calibration size, covariance estimation, and the concentration failure probability. The term ρtwin _t^win accounts for selecting the wrong window without conditioning on data-dependent search paths [36, 18]. Corollary 1 follows by averaging over t and applying the tower property. Corollary 1 (Expected empirical average coverage). Under the conditions of Theorem 1, assume the per-time right-hand-side terms are deterministic or hold almost surely under the outer expectation. Let Cov^T=T−1∑t=1TYt∈tACI(α) Cov_T=T^-1 _t=1^T1\Y_t _t^ACI(α)\, and define ρ¯Twin=T−1∑t=1Tρtwin ρ_T^win=T^-1 _t=1^T _t^win, ε¯TACI=T−1∑t=1TεtACI _T^ACI=T^-1 _t=1^T _t^ACI, η¯T=T−1∑t=1TηM,δt/(nt⋆+1) η_T=T^-1 _t=1^T _M, _t/(n_t +1), and δ¯T=T−1∑t=1Tδt δ_T=T^-1 _t=1^T _t. Then |⁡[Cov^T]−(1−α)|≤ρ¯Twin+ε¯TACI+1T∑t=1T1nt⋆+1+LΣη¯T+δ¯T. |E[ Cov_T]-(1-α) |≤ ρ_T^win+ _T^ACI+ 1T _t=1^T 1n_t +1+L_ η_T+ δ_T. (17) If these terms are random, interpret the averaged right-hand-side terms as their outer expectations. 5 Experiments 5.1 Experimental setup and metrics We evaluate Space on multivariate forecasting streams produced by eight heterogeneous probabilistic forecasters (TimeGrad [27], CSDI [33], TSFlow [14], DeepAR [31], GP-Copula [30], NsDiff [43], TMDM [19], and TACTiS-2 [3]) across seven real-world datasets summarized in Table 8. We then compare the performance of Space against five conformal baselines (MultiDimSPCI [40], Local Ellipsoid [24], CopulaCPTS [32], Empirical Copula and Gaussian Copula [23]). The probabilistic forecasters are trained with a 7:1:2 train/validation/test split on the original dataset. Each conformal wrapper is evaluated with a 6:2:2 train/calibration/test split on the sample streams. In downstream decision settings, predictive uncertainty should be close to the nominal target, rather than merely exceed it [44, 15, 8]. We therefore report both the coverage gap and the rolling coverage gap. The coverage gap is the absolute deviation between empirical joint coverage and the target level 1−α1-α, and therefore penalizes both under- and over-coverage. To evaluate efficiency, we report the mean log-volume of the prediction regions, denoted LogVol, following prior multivariate conformal work that measures sharpness through region size or volume [40, 17, 42]. This section focuses on the target coverage level 1−α=0.901-α=0.90. The full results at different target coverage levels and other implementation details are deferred to Appendix C and F. 5.2 Aggregate comparison against conformal wrappers We first compare Space with five strong conformal wrappers from prior work. Table 2 provides the wrapper-level comparison. The full dataset-forecaster level comparison is in Appendix C. Table 2: Aggregated comparison over all dataset–forecaster pairs. Mean Gap and Median Gap summarize the pair-level absolute coverage gap from the target. Std. LogVol denotes the pair-level log-volume after standardization within each dataset and target level, so it measures size on a comparable scale across datasets. For Mean Gap and Mean Std. LogVol, the reported standard deviation is the across-pair standard deviation within each wrapper. Mean Std. LogVol and Median Std. LogVol are computed only over pairs whose coverage gap is at most 0.050.05. Valid Ratio reports the fraction of dataset–forecaster pairs satisfying this coverage-gap threshold. Wrapper Mean Gap ↓ Median Gap ↓ Mean Std. LogVol ↓ Median Std. LogVol ↓ Valid Ratio SPACE 0.003 ± 0.004 0.002 -0.111 ± 1.043 -0.334 1.000 MultiDimSPCI 0.032 ± 0.041 0.013 -0.435 ± 1.024 -0.804 0.759 Local Ellipsoid 0.083 ± 0.111 0.054 -0.580 ± 0.878 -0.753 0.463 Empirical Copula 0.165 ± 0.205 0.053 0.031 ± 0.630 0.012 0.481 Gaussian Copula 0.130 ± 0.184 0.042 0.003 ± 0.618 -0.037 0.593 CopulaCPTS 0.057 ± 0.112 0.034 1.021 ± 0.789 1.037 0.593 Space achieves the best Mean Gap and Median Gap among all compared wrappers, with the smallest standard deviation. This is the main wrapper-level result: the mean coverage gap of Space is 0.30.3 percentage points from the nominal target. The efficiency results are more nuanced. Space does not attain the smallest LogVol overall; MultiDimSPCI has the lowest median volume. However, this smaller volume comes with noticeably worse calibration and a lower validity ratio. In particular, the mean coverage gap of Space is only about one tenth of that of the strongest baseline, MultiDimSPCI, and Space is the only method with a validity ratio of 1.001.00. This indicates that Space is the most consistent wrapper for bringing realized coverage close to the nominal target across diverse datasets and forecasters. The key point is therefore coverage–efficiency trade-off at the target level. Space achieves much better calibration while remaining competitive in median log-volume, which is the more relevant comparison for a conformal wrapper whose purpose is to deliver regions that are both reliable and usable. 5.3 Rolling coverage under drift and regime change Rolling metrics show whether a wrapper stays close to the nominal target locally over time. This is a pattern that could be overlooked by aggregated mean or median measures. We use a rolling window length w=30w=30 in all evaluations. Figure 3 reports boxplots for Mean Rolling Gap, P90 Rolling Gap, and Frac. Bad Windows, where lower values indicate better local calibration stability. Mean Rolling Gap averages the absolute deviation between rolling joint coverage and the target level 1−α1-α over all valid sliding windows. P90 Rolling Gap reports the 9090th percentile of these deviations, capturing severe local miscalibration. Frac. Bad Windows is the proportion of windows whose rolling coverage gap exceeds δ=0.10δ=0.10. Table 3 complements the figure by providing the dataset-level breakdown. Figure 3: Distributional boxplots of rolling coverage diagnostics across wrapper benchmark pairs. The center line denotes the median. The three panels report Mean Rolling Gap, P90 Rolling Gap, and Frac. Bad Windows, respectively, with lower values indicating better local calibration stability. Table 3 gives the corresponding dataset-level breakdown. Table 3: Dataset-level rolling coverage diagnostics for conformal wrappers at the 90%90\% target. Lower is better for all metrics. These metrics evaluate local calibration stability under drift and regime change. Space performs best or near-best on all datasets. Dataset Metric SPACE MultiDimSPCI Local Ell. Emp. Copula Gauss. Copula CopulaCPTS electricity Mean Rolling Gap ↓ 0.041 0.109 0.278 0.326 0.125 0.082 P90 Rolling Gap ↓ 0.086 0.214 0.557 0.556 0.257 0.100 Frac. Bad Windows ↓ 0.021 0.360 0.584 0.689 0.298 0.010 etth1 Mean Rolling Gap ↓ 0.032 0.044 0.074 0.048 0.049 0.072 P90 Rolling Gap ↓ 0.067 0.092 0.133 0.096 0.096 0.112 Frac. Bad Windows ↓ 0.006 0.014 0.074 0.020 0.010 0.041 ettm1 Mean Rolling Gap ↓ 0.054 0.051 0.083 0.062 0.063 0.107 P90 Rolling Gap ↓ 0.108 0.100 0.146 0.112 0.108 0.204 Frac. Bad Windows ↓ 0.065 0.047 0.067 0.054 0.038 0.185 ettm2 Mean Rolling Gap ↓ 0.045 0.055 0.070 0.067 0.066 0.065 P90 Rolling Gap ↓ 0.096 0.100 0.146 0.112 0.108 0.117 Frac. Bad Windows ↓ 0.029 0.061 0.118 0.091 0.072 0.077 exchange Mean Rolling Gap ↓ 0.050 0.067 0.156 0.257 0.250 0.211 P90 Rolling Gap ↓ 0.100 0.163 0.263 0.444 0.436 0.313 Frac. Bad Windows ↓ 0.066 0.162 0.225 0.460 0.435 0.296 traffic Mean Rolling Gap ↓ 0.055 0.069 0.103 0.467 0.467 0.072 P90 Rolling Gap ↓ 0.100 0.157 0.205 0.676 0.676 0.105 Frac. Bad Windows ↓ 0.039 0.174 0.281 0.950 0.950 0.038 weather Mean Rolling Gap ↓ 0.052 0.073 0.121 0.198 0.167 0.111 P90 Rolling Gap ↓ 0.096 0.142 0.288 0.488 0.408 0.204 Frac. Bad Windows ↓ 0.051 0.136 0.173 0.441 0.344 0.166 Across datasets, Space shows its clear advantage over the best baseline in mean rolling gap, especially on the higher-dimensional, more nonstationary settings, where maintaining stable coverage over time is most difficult. In particular, the separation between Space and the best competing method is largest on datasets such as electricity and traffic, which have richer joint dependence structure and stronger temporal instability; this suggests that Space is especially effective when the calibration problem is genuinely multivariate and the underlying coverage behavior drifts over time. By contrast, on smaller and more stationary datasets, the gap between methods narrows substantially. In ettm1, Space and the best baseline perform similarly. A more detailed analysis in Appendix D shows that Space’s gains are more pronounced in datasets with stronger covariance drift and more challenging rolling coverage behavior. 6 Conclusion In this paper, we introduced Space, a post-hoc conformal wrapper that bridges the gap between the rich uncertainty representations of modern probabilistic forecasters and the rigorous calibration guarantees of conformal prediction. By decoupling the estimation of multivariate geometry—extracted directly from the predictive sample cloud—from the sequential calibration of the region’s radius, Space overcomes the adaptation delays inherent in historical residual-based methods. Our theoretical analysis establishes formal finite-sample coverage guarantees under stated assumptions, and extensive experiments demonstrate that Space consistently tightens realized coverage toward the nominal target across diverse forecasting architectures (e.g., diffusion, flow, and autoregressive models) under non-stationarity. Ultimately, this work demonstrates that predictive sample clouds contain valuable geometric signals that can be exploited for calibrated predictive inference. Natural extensions of this framework include generalizing the set geometry beyond ellipsoids to capture non-convex predictive topologies. 7 Limitations This work has two main limitations. First, Space depends on the quality of the predictive sample cloud produced by the generative forecasters. The conformal calibration can correct the overall scale of the prediction region, but the local covariance geometry is still estimated from generated samples. Therefore, if the forecaster produces biased or insufficiently diverse samples, the calibrated regions may remain valid only at the cost of reduced efficiency. The second limitation concerns adaptation under difficult regime dynamics. When changes are gradual, recent scores may not show a sharp enough diagnostic signal; when changes are frequent, little same-regime history may be available. In either case, Space may select a window that is not fully representative of the current regime, leading to coverage deviations. This challenge is shared by time-series uncertainty calibration methods without external change indicators, where adaptation must be inferred from observed errors. References [1] A. Angelopoulos, E. Candes, and R. J. Tibshirani (2023) Conformal PID control for time series prediction. Advances in Neural Information Processing Systems 36, p. 23047–23074. Cited by: item (A4), item (A4), §1, §2, Table 1, §4. [2] A. N. Angelopoulos and S. Bates (2021) A gentle introduction to conformal prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.07511. Cited by: item (A3), item (A3), item (A4), §1. [3] A. Ashok, É. Marcotte, V. Zantedeschi, N. Chapados, and A. Drouin (2024) TACTis-2: better, faster, simpler attentional copulas for multivariate time series. In The Twelfth International Conference on Learning Representations, Cited by: Table 4, Table 4, §1, §2, §5.1. [4] R. F. Barber, E. J. Candes, A. Ramdas, and R. J. Tibshirani (2023) Conformal prediction beyond exchangeability. The Annals of Statistics 51 (2), p. 816–845. Cited by: §2, §4. [5] M. Basseville I. V. Nikiforov et al. (1993) Detection of abrupt changes: theory and application. Vol. 104, Prentice Hall, Englewood Cliffs. Cited by: item (A2), item (A2). [6] A. Bhatnagar, H. Wang, C. Xiong, and Y. Bai (2023) Improved online conformal prediction via strongly adaptive online learning. In International Conference on Machine Learning, p. 2337–2363. Cited by: item (A4), item (A4), §2, §4, §4. [7] T. T. Cai, C. Zhang, and H. H. Zhou (2010) Optimal rates of convergence for covariance matrix estimation. The Annals of Statistics 38 (4), p. 2118–2144. Cited by: §A.2, §A.2. [8] L. Deng, H. Xiong, F. Wu, S. Kapoor, S. Gosh, Z. Shahn, and L. Lehman (2025) Uncertainty quantification for conditional treatment effect estimation under dynamic treatment regimes. In Proceedings of the 4th Machine Learning for Health Symposium, Vol. 259, p. 248–266. Cited by: §5.1. [9] A. Dvoretzky, J. Kiefer, and J. Wolfowitz (1956) Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator. The Annals of Mathematical Statistics, p. 642–669. Cited by: item (A2), item (A2), item (A2), footnote 2. [10] S. Feldman, S. Bates, and Y. Romano (2023) Calibrated multiple-output quantile regression with representation learning. Journal of Machine Learning Research 24 (24), p. 1–48. Cited by: §2. [11] I. Gibbs and E. Candes (2021) Adaptive conformal inference under distribution shift. Advances in Neural Information Processing Systems 34, p. 1660–1672. Cited by: item (A4), item (A4), §1, §1, §2, Table 1, §3.3, §4. [12] L. Guan (2023) Localized conformal prediction: a generalized inference framework for conformal prediction. Biometrika 110 (1), p. 33–50. Cited by: §2. [13] C. Johnstone and E. Ndiaye (2025) Exact and approximate conformal inference for multi-output regression. In Proceedings of the Fourteenth Symposium on Conformal and Probabilistic Prediction with Applications, Vol. 266, p. 153–172. Cited by: §2. [14] M. Kollovieh, M. Lienen, D. Lüdke, L. Schwinn, and S. Günnemann (2025) Flow matching with gaussian process priors for probabilistic time series forecasting. In The Thirteenth International Conference on Learning Representations, Cited by: Table 4, Table 4, §1, §2, §5.1. [15] V. Kuleshov and S. Deshpande (2022) Calibrated and sharp uncertainties in deep learning via density estimation. In International Conference on Machine Learning, p. 11683–11693. Cited by: §5.1. [16] G. Lai, W. Chang, Y. Yang, and H. Liu (2018) Modeling long- and short-term temporal patterns with deep neural networks. In The 41st international ACM SIGIR conference on research & development in information retrieval, p. 95–104. Cited by: item 3, item 4. [17] J. Lee, C. Xu, and Y. Xie (2026) Flow-based conformal prediction for multi-dimensional time series. In The Fourteenth International Conference on Learning Representations, Cited by: §2, Table 1, §5.1. [18] J. Lei, M. G’Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman (2018) Distribution-free predictive inference for regression. Journal of the American Statistical Association 113 (523), p. 1094–1111. Cited by: item (A3), item (A3), item (A4), §1, §4. [19] Y. Li, W. Chen, X. Hu, B. Chen, B. Sun, and M. Zhou (2024) Transformer-modulated diffusion models for probabilistic multivariate time series forecasting. In The Twelfth International Conference on Learning Representations, Cited by: Table 4, Table 4, §1, §2, §5.1. [20] C. Marx, S. Zalouk, and S. Ermon (2023) Calibration by distribution matching: trainable kernel calibration metrics. In Advances in Neural Information Processing Systems, Vol. 36, p. 25910–25928. Cited by: §1. [21] P. Massart (1990) The tight constant in the dvoretzky-kiefer-wolfowitz inequality. The Annals of Probability, p. 1269–1283. Cited by: item (A2), item (A2), item (A2), footnote 2. [22] D. S. Matteson and N. A. James (2014) A nonparametric approach for multiple change point analysis of multivariate data. Journal of the American Statistical Association 109 (505), p. 334–345. Cited by: Figure 5, Figure 5, Appendix E. [23] S. Messoudi, S. Destercke, and S. Rousseau (2021) Copula-based conformal prediction for multi-target regression. Pattern Recognition 120, p. 108101. Cited by: §2, Table 1, §5.1. [24] S. Messoudi, S. Destercke, and S. Rousseau (2022) Ellipsoidal conformal inference for multi-target regression. In Conformal and Probabilistic Prediction with Applications, p. 294–306. Cited by: §1, §2, §5.1. [25] T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida (2018) Spectral normalization for generative adversarial networks. In International Conference on Learning Representations, Cited by: footnote 1. [26] S. S. Rangapuram, M. W. Seeger, J. Gasthaus, L. Stella, Y. Wang, and T. Januschowski (2018) Deep state space models for time series forecasting. Advances in Neural Information Processing Systems 31. Cited by: §1. [27] K. Rasul, C. Seward, I. Schuster, and R. Vollgraf (2021) Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting. In International Conference on Machine Learning, p. 8857–8868. Cited by: Table 4, Table 4, §1, §2, §5.1. [28] K. Rasul, A. Sheikh, I. Schuster, U. M. Bergmann, and R. Vollgraf (2021) Multivariate probabilistic time series forecasting via conditioned normalizing flows. In International Conference on Learning Representations, Cited by: §1. [29] Y. Romano, E. Patterson, and E. Candes (2019) Conformalized quantile regression. Advances in Neural Information Processing Systems 32. Cited by: item (A3), §1. [30] D. Salinas, M. Bohlke-Schneider, L. Callot, R. Medico, and J. Gasthaus (2019) High-dimensional multivariate forecasting with low-rank gaussian copula processes. Advances in Neural Information Processing Systems 32. Cited by: Table 4, Table 4, §1, §2, §5.1. [31] D. Salinas, V. Flunkert, J. Gasthaus, and T. Januschowski (2020) DeepAR: probabilistic forecasting with autoregressive recurrent networks. International journal of forecasting 36 (3), p. 1181–1191. Cited by: Table 4, Table 4, §1, §2, §5.1. [32] S. H. Sun and R. Yu (2024) Copula conformal prediction for multi-step time series prediction. In The Twelfth International Conference on Learning Representations, Cited by: §2, Table 1, §5.1. [33] Y. Tashiro, J. Song, Y. Song, and S. Ermon (2021) CSDI: conditional score-based diffusion models for probabilistic time series imputation. Advances in Neural Information Processing Systems 34, p. 24804–24816. Cited by: Table 4, Table 4, §2, §5.1. [34] R. J. Tibshirani, R. Foygel Barber, E. Candes, and A. Ramdas (2019) Conformal prediction under covariate shift. Advances in Neural Information Processing Systems 32. Cited by: §1, §2. [35] R. Vershynin (2018) High-dimensional probability: an introduction with applications in data science. Vol. 47, Cambridge University Press. Cited by: item (A1), §A.2, §4, footnote 1. [36] V. Vovk, A. Gammerman, and G. Shafer (2005) Algorithmic learning in a random world. Springer. Cited by: item (A3), item (A3), item (A3), item (A4), §1, §4. [37] M. J. Wainwright (2019) High-dimensional statistics: a non-asymptotic viewpoint. Vol. 48, Cambridge University Press. Cited by: §4. [38] H. Wu, J. Xu, J. Wang, and M. Long (2021) Autoformer: decomposition transformers with auto-correlation for long-term series forecasting. Advances in Neural Information Processing Systems 34, p. 22419–22430. Cited by: item 5. [39] Y. Xie, J. Huang, and R. Willett (2013) Change-point detection for high-dimensional time series with missing data. IEEE Journal of Selected Topics in Signal Processing 7 (1), p. 12–27. Cited by: item (A2), item (A2), §4, §4. [40] C. Xu, H. Jiang, and Y. Xie (2024) Conformal prediction for multi-dimensional time series by ellipsoidal sets. In Proceedings of the 41st International Conference on Machine Learning, Vol. 235, p. 55076–55099. Cited by: §1, §1, §2, Table 1, §4, §5.1, §5.1. [41] C. Xu and Y. Xie (2021) Conformal prediction interval for dynamic time-series. In Proceedings of the 38th International Conference on Machine Learning, Vol. 139, p. 11559–11569. Cited by: §1, §2, Table 1. [42] Q. Yang, Q. J. Zhu, J. Giezendanner, Y. Marzouk, S. Bates, and S. Wang (2026) Conformal prediction for generative models via adaptive cluster-based density estimation. arXiv preprint arXiv:2601.22298. Cited by: §5.1. [43] W. Ye, Z. Xu, and N. Gui (2025) Non-stationary diffusion for probabilistic time series forecasting. In Proceedings of the 42nd International Conference on Machine Learning, Vol. 267, p. 72112–72130. Cited by: Table 4, Table 4, §1, §2, §5.1. [44] S. Zhao, M. Kim, R. Sahoo, T. Ma, and S. Ermon (2021) Calibrating predictions to decisions: a novel approach to multi-class calibration. Advances in Neural Information Processing Systems 34, p. 22313–22324. Cited by: §5.1. [45] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang (2021) Informer: beyond efficient transformer for long sequence time-series forecasting. In Proceedings of the AAAI conference on artificial intelligence, Vol. 35, p. 11106–11115. Cited by: item 2. Appendix A Proofs for Section 4 A.1 Assumptions. We work under the following four assumptions. They are stated in the same order and with the same meaning as in Section 4. Let ℱtF_t be the natural filtration generated by realized targets, forecast sample clouds, and all auxiliary randomization up to time t. Let τt<t _t<t denote the last change point before prediction time t. Let ℚt(⋅)Q_t(·) denote the probability measure conditional on the pre-regime information ℱτtF_ _t and the deterministic window-design parameters, but marginal over the realized augmented objects inside the current regime. Let ℚt[⋅]E_Q_t[·] be the corresponding expectation. (A1) Conditional sample-cloud regularity and stable score geometry. We assume that, at each prediction time t, the forecast sample cloud is conditionally i.i.d. and sufficiently light-tailed after standardization. We also assume that the Mahalanobis conformal set is locally stable with respect to covariance-estimation error, in a pathwise sense conditional on the realized covariance estimates. These regularity conditions are standard in high-dimensional probability and are adapted from the foundational framework established by Vershynin [35]. More formally: (A1a) Conditional sample-cloud law. Conditional on ℱt−1F_t-1, y^t(1),…,y^t(M)∼i.i.d.P^t(⋅∣ℱt−1). y_t^(1),…, y_t^(M) i.i.d. P_t(· _t-1). (18) Let D=diag⁡(s1,…,sd)D=diag(s_1,…,s_d) be fixed before time t, with 0<smin≤sj≤smax<∞.0<s_ ≤ s_j≤ s_ <∞. For an ℱt−1F_t-1-measurable center μt⋆ _t , define ξt(m)=D−1(y^t(m)−μt⋆),m=1,…,M. _t^(m)=D^-1 ( y_t^(m)- _t ), m=1,…,M. (A1b) Sub-Gaussian residuals and covariance estimation target. Conditional on ℱt−1F_t-1, the vectors ξt(1),…,ξt(M) _t^(1),…, _t^(M) are mean-zero, K-sub-Gaussian, and have conditional covariance Ct⋆=Cov⁡(ξt(m)∣ℱt−1).C_t =Cov\! ( _t^(m) _t-1 ). (19) The empirical covariance C^t C_t in (3) is the centered empirical covariance of the standardized sample cloud. Hence it is unchanged by replacing the ideal center μt⋆ _t with the median-based center μt _t used in the method. (A1c) Pathwise stability of the conformal set under covariance perturbation. For any deterministic or previsible same-regime calibration index set ℋtH_t, and in particular for the oracle set ℋt⋆H_t , let t⋆(αt)C_t ( _t) denote the conformal set constructed from ideal scores, and let t(αt)C_t( _t) denote the conformal set constructed from empirical scores. For such a set ℋtH_t, define the conditional law ℚ~t,ℋt(⋅)=ℙ(⋅∣ℱτt,(C^s)s∈ℋt∪t). Q_t,H_t(·)=P\! (· _ _t,( C_s)_s _t∪\t\ ). There exist constants η0>0 _0>0 and LΣ<∞L_ <∞ such that, on the event maxs∈ℋt∪t⁡‖C^s−Cs⋆‖op≤η≤η0, _s _t∪\t\\| C_s-C_s \|_op≤η≤ _0, we have |ℚ~t,ℋt(Yt∈t(αt))−ℚ~t,ℋt(Yt∈t⋆(αt))|≤LΣη. | Q_t,H_t\! (Y_t _t( _t) )- Q_t,H_t\! (Y_t _t ( _t) ) |≤ L_ η. (20) This condition is implied, for example, by local stability of the stabilized pseudo-inverse map together with an anti-concentration condition for the ideal radial score distribution. (A2) Same-regime test calibration and contaminated-window detectability. A2 is a finite-sample version of two standard requirements in change-point and sequential goodness-of-fit testing: size control under homogeneous samples and power under separated alternatives. Fix a prediction time t, and let τt _t denote the most recent regime change before time t. Thus, the current regime contains the observations indexed by τt+1,…,t−1.\ _t+1,…,t-1\. Let ℚtQ_t denote the probability law, or the relevant conditional probability law, governing the data and random test statistics at time t. This assumption states that same-regime candidate blocks are not rejected too often, while any candidate window that crosses a regime boundary is detected with high probability. The assumption is stated directly at the level of the finite-sample test errors used by the window-selection algorithm. This calibration-versus-detectability form is standard in change-point analysis: homogeneous samples require size control, while contaminated samples require a separation condition large enough to overcome sampling noise [5, 39]. More formally: (A2a) Oracle same-regime window. Let =Lmin,Lmin+h,…,LmaxG=\L_ ,L_ +h,…,L_ \ be the finite grid of candidate window lengths, where LminL_ and LmaxL_ are the minimum and maximum candidate lengths, respectively, and h is the grid step size. Let tP_t denote the probe block at time t, namely the recent block used as the reference block for testing whether a candidate block belongs to the same regime. Assume that tP_t lies inside the current regime and that there exists at least one same-regime candidate length. For each L∈L , let twhole(L)C_t^whole(L) denote the whole candidate block of length L used by the algorithm. Define the oracle same-regime length Lt⋆L_t as the largest candidate length whose whole candidate block and probe block are both contained in the current regime: Lt⋆=max⁡L∈:twhole(L)∪t⊆τt+1,…,t−1.L_t = \L :C_t^whole(L) _t \ _t+1,…,t-1\ \. (21) For each candidate block ∈twhole(L),ℬtinc(L),A∈\C_t^whole(L),B_t^inc(L)\, where ℬtinc(L)B_t^inc(L) denotes the incremental block associated with candidate length L, and for each diagnostic j∈1,2,3j∈\1,2,3\, let KSt,L,jKS_t,L,j^A denote the empirical KS-type discrepancy in (10). Here j indexes the diagnostic being tested, L indexes the candidate length, and A indexes the tested candidate block. Let τt,L,j _t,L,j^A denote the corresponding rejection threshold for this KS-type test. The thresholds are deterministic functions of the window-design parameters, sample sizes, and nominal test levels, or are measurable with respect to information outside the realized test statistics. For the KS-type diagnostics used here, such thresholds can be justified by the Dvoretzky–Kiefer–Wolfowitz inequality and Massart’s sharp constant, which give finite-sample uniform bounds for empirical distribution functions [9, 21]. (A2b) Clean-block calibration. For every clean tested block, meaning every tested block A whose union with the probe block tP_t lies entirely inside the current regime, ∪t⊆τt+1,…,t−1,A _t \ _t+1,…,t-1\, the probability of false rejection is controlled: ℚt(KSt,L,j>τt,L,j)≤δuni.Q_t (KS_t,L,j^A> _t,L,j^A )≤ _uni. (22) Here δuni _uni is the per-test, or union-adjusted, upper bound on the false rejection probability. This is the finite-sample size-control requirement for the same-regime test. When A and tP_t are drawn from the same regime, the population discrepancy is zero, and the empirical discrepancy is controlled by standard empirical-process concentration; the constant δuni _uni can also absorb a union bound over the finite set of candidate lengths, blocks, and diagnostics [9, 21]. (A2c) Contaminated-window detectability. For every L>Lt⋆L>L_t , the candidate length L is longer than the oracle same-regime length and hence its associated candidate window necessarily crosses a regime boundary. The assumption requires that at least one tested pair (,j),∈twhole(L),ℬtinc(L),j∈1,2,3,(A,j), ∈\C_t^whole(L),B_t^inc(L)\, j∈\1,2,3\, detects this contamination. Specifically, the probability that this contaminated tested pair fails to reject is bounded by βdet _det: ℚt(KSt,L,j≤τt,L,j)≤βdet.Q_t (KS_t,L,j^A≤ _t,L,j^A )≤ _det. (23) Here βdet _det is the upper bound on the missed-detection probability. This is the corresponding finite-sample power requirement. It is a local identifiability condition: if a candidate window crosses a regime boundary, then at least one of the diagnostics must see a population distributional shift large enough to be detected with probability at least 1−βdet1- _det, as in standard abrupt-change and high-dimensional change-point formulations [5, 39]. A sufficient, but not necessary, condition for (23) is the following. For some tested pair (,j)(A,j), let Ht,L,jH_t,L,j^A denote the population distribution function, or population probability-integral-transform target, of the diagnostic values comparing the candidate block A with the probe block tP_t. Let Id⁡(x)=xId(x)=x denote the identity distribution function on [0,1][0,1], which corresponds to the ideal same-regime case after probability-integral transformation. Let F^t,L,j F_t,L,j^A denote the empirical distribution function estimated from the probe block, and let Δsep>0 _sep>0 denote the separation margin between the population discrepancy and the test threshold. Suppose ‖Ht,L,j−Id‖∞≥τt,L,j+Δsep,\|H_t,L,j^A-Id\|_∞≥ _t,L,j^A+ _sep, where ∥⋅∥∞\|·\|_∞ denotes the uniform sup-norm distance, and suppose that the empirical distribution function satisfies the concentration bound ℚt(‖F^t,L,j−Ht,L,j‖∞>Δsep)≤2exp⁡−2|t|Δsep2.Q_t (\| F_t,L,j^A-H_t,L,j^A\|_∞> _sep )≤ 2 \-2|P_t| _sep^2\. Here |t||P_t| is the number of observations in the probe block. Then (23) holds with βdet=2exp⁡−2|t|Δsep2. _det=2 \-2|P_t| _sep^2\. Indeed, the population separation margin Δsep _sep ensures that the empirical KS statistic can fall below the threshold only if the empirical process deviates from its population target by more than Δsep _sep. The displayed exponential bound is therefore exactly the DKW/Massart-type tail probability for this failure event [9, 21]. (A3) Piecewise augmented exchangeability. This assumption imposes the standard exchangeability condition used in conformal prediction, but only locally within the current regime. Specifically, conditional on the pre-regime information ℱτtF_ _t, where τt _t is the most recent regime change before time t, the observations inside the current regime are treated as exchangeable. This is the classical condition under which conformal methods obtain finite-sample rank validity and distribution-free coverage [36, 18, 2]. Here the condition is imposed piecewise, after conditioning on the regime start, which is natural in nonstationary settings where exchangeability is expected to hold only locally rather than over the entire time series. More formally, for each s=τt+1,…,ts= _t+1,…,t, define the augmented forecast-target object s=(Ys,^s),Z_s=(Y_s, Y_s), where YsY_s is the realized target and ^s Y_s is the corresponding forecast object, such as a predictive distribution, trajectory forecast, or prediction set. The augmented objects τt+1,…,tZ_ _t+1,…,Z_t are conditionally exchangeable given ℱτtF_ _t. That is, for any finite permutation π of τt+1,…,t\ _t+1,…,t\, (τt+1,…,t)∣ℱτt=(π⁡(τt+1),…,π⁡(t))∣ℱτt.(Z_ _t+1,…,Z_t) _ _t d=(Z_π( _t+1),…,Z_π(t)) _ _t. This is the regime-local analogue of the exchangeability assumption commonly used in split conformal prediction and related finite-sample conformal procedures [36, 18, 29]. The diagnostics and ideal conformity scores are assumed to be generated from these augmented objects by common measurable permutation-equivariant maps, possibly together with auxiliary tie-breaking random variables that are independent of the augmented objects conditional on ℱτtF_ _t. Permutation equivariance means that if the augmented objects are permuted, the resulting diagnostics or scores are permuted in the same way, rather than changing their joint distribution. Therefore, by the standard closure property of exchangeability under common measurable equivariant transformations, the ideal score collection for any same-regime index set ℋtH_t, es⋆:s∈ℋt∪et⋆,\e_s :s _t\∪\e_t \, is conditionally exchangeable under ℚtQ_t. This is the exact property needed in the conformal coverage proof: the test score et⋆e_t has a uniform rank among the same-regime calibration scores, up to the usual tie-breaking convention [36, 2]. (A4) ACI adaptation error and rank validity. Let t,or⋆(αt)C_t,or ( _t) denote the oracle same-regime conformal set constructed from the ideal scores on ℋt⋆H_t using the adaptive miscoverage level αt∈[0,1] _t∈[0,1], and let nt⋆=|ℋt⋆|n_t =|H_t | be the number of oracle same-regime calibration scores. This assumption separates the effect of adaptive level selection from the usual conformal rank-validity argument. Such a separation is standard in adaptive conformal inference: the adaptive procedure controls or tracks the target miscoverage level over time, while the conformal coverage argument relies on the rank of the test score among the calibration scores [11, 1, 6]. (A4a) Mean tracking error. The adaptive level remains close to the nominal level in expectation: ℚt[|αt−α|]≤εtACI.E_Q_t\! [| _t-α| ]≤ _t^ACI. (24) Here α is the target miscoverage level, αt _t is the data-adaptive level selected by the ACI update rule at time t, and εtACI _t^ACI summarizes the finite-sample adaptation or tracking error of this update. This type of condition is the standard way to abstract the performance of an adaptive conformal update under distribution shift: rather than requiring αt=α _t=α exactly at every time point, one allows the online update to deviate from the target level by a controlled amount [11, 1]. Strongly adaptive online-learning formulations provide one route for bounding such tracking errors over nonstationary sequences [6]. (A4b) Rank validity under adaptivity. Conditional on αt _t, the rank of the test ideal score among the nt⋆+1n_t +1 oracle ideal scores es⋆:s∈ℋt⋆∪et⋆\e_s :s _t \∪\e_t \ is uniformly distributed. Equivalently, the oracle ideal conformal predictor satisfies the finite-sample rank bound |ℚt(Yt∈t,or⋆(αt))−(1−ℚt[αt])|≤1nt⋆+1. |Q_t\! (Y_t _t,or ( _t) )- (1-E_Q_t[ _t] ) |≤ 1n_t +1. (25) This is the usual finite-sample conformal rank argument, applied after conditioning on the adaptive level. It is satisfied, for example, when αt _t is predictable with respect to the score-ranking step, meaning that it is determined by past information, previous coverage errors, or other external randomness, but not by the relative rank of et⋆e_t among the same-regime calibration scores at time t. Under this previsible-level condition, conditioning on αt _t does not destroy the exchangeability of the oracle ideal scores, so the test score keeps the usual uniform-rank property [36, 18, 2]. The additional 1/(nt⋆+1)1/(n_t +1) term is the standard finite-sample discretization error from using nt⋆n_t calibration scores in conformal prediction. Combining (24) and (25), we have |ℚt(Yt∈t,or⋆(αt))−(1−α)|≤εtACI+1nt⋆+1. |Q_t\! (Y_t _t,or ( _t) )-(1-α) |≤ _t^ACI+ 1n_t +1. (26) Thus, A4 is not an additional distributional assumption beyond the local exchangeability condition in A3. Rather, it records two standard ingredients needed when conformal prediction is combined with an adaptive level: a bounded ACI tracking error and the preservation of the conformal rank argument under a previsible adaptive level. A.2 Proof of Proposition 1 Our goal is to establish the concentration of the empirical sample-cloud covariance. While this is a foundational step for our method, the proof relies on well-established techniques rather than novel probabilistic machinery. Specifically, we follow the standard high-dimensional probability pipeline detailed in Vershynin [35]. Proof. Fix a prediction time t and condition on ℱt−1F_t-1. Write ℙt(⋅)=ℙ(⋅∣ℱt−1)P_t(·)=P(· _t-1) and t[⋅]=[⋅∣ℱt−1]E_t[·]=E[· _t-1]. For brevity, drop the time subscript and write ξm=D−1(y^t(m)−μt⋆),m=1,…,M. _m=D^-1 ( y_t^(m)- _t ), m=1,…,M. By Assumption (A1), conditional on ℱt−1F_t-1, the vectors ξ1,…,ξM _1,…, _M are independent, mean-zero, and sub-Gaussian with sub-Gaussian parameter at most K, and Ct⋆=t[ξ1ξ1⊤].C_t =E_t[ _1 _1 ]. Let ξ¯=1M∑m=1Mξm,C^t=1M∑m=1M(ξm−ξ¯)(ξm−ξ¯)⊤. ξ= 1M _m=1^M _m, C_t= 1M _m=1^M( _m- ξ)( _m- ξ) . Then C^t=1M∑m=1Mξmξm⊤−ξ¯ξ¯⊤. C_t= 1M _m=1^M _m _m - ξ ξ . Therefore, ‖C^t−Ct⋆‖op≤‖1M∑m=1Mξmξm⊤−Ct⋆‖op+‖ξ¯‖22.\| C_t-C_t \|_op≤ \| 1M _m=1^M _m _m -C_t \|_op+\| ξ\|_2^2. (27) We first bound the uncentered empirical second-moment term. For any symmetric matrix A, if N is a 1/41/4-net of d−1S^d-1, then ‖A‖op≤2maxv∈|v⊤Av|.\|A\|_op≤ 2 _v |v Av|. Choose such a net with cardinality ||≤9d|N|≤ 9^d. Applying this to A=1M∑m=1Mξmξm⊤−Ct⋆,A= 1M _m=1^M _m _m -C_t , we obtain ‖1M∑m=1Mξmξm⊤−Ct⋆‖op≤2maxv∈|1M∑m=1M⟨ξm,v⟩2−t[⟨ξ1,v⟩2]|. \| 1M _m=1^M _m _m -C_t \|_op≤ 2 _v | 1M _m=1^M _m,v ^2-E_t[ _1,v ^2] |. (28) Fix v∈v . Since ξm _m is conditionally K-sub-Gaussian, ⟨ξm,v⟩ _m,v is conditionally K-sub-Gaussian. Hence Zm,v=⟨ξm,v⟩2−t[⟨ξm,v⟩2]Z_m,v= _m,v ^2-E_t[ _m,v ^2] is conditionally centered sub-exponential with sub-exponential norm bounded by c1K2c_1K^2, for an absolute constant c1>0c_1>0. Bernstein’s inequality for independent centered sub-exponential random variables gives, for every u>0u>0, ℙt(|1M∑m=1MZm,v|≥u)≤2exp⁡[−c2Mmin⁡(u2K4,uK2)],P_t ( | 1M _m=1^MZ_m,v |≥ u )≤ 2 [-c_2M ( u^2K^4, uK^2 ) ], where c2>0c_2>0 is an absolute constant. Taking a union bound over v∈v yields ℙt(maxv∈⁡|1M∑m=1MZm,v|≥u)≤2exp⁡[dlog⁡9−c2Mmin⁡(u2K4,uK2)].P_t ( _v | 1M _m=1^MZ_m,v |≥ u )≤ 2 [d 9-c_2M ( u^2K^4, uK^2 ) ]. Thus, for a sufficiently large absolute constant c3>0c_3>0, with ℙtP_t-probability at least 1−δ/21-δ/2, ‖1M∑m=1Mξmξm⊤−Ct⋆‖op≤c3K2(d+log⁡(1/δ)M+d+log⁡(1/δ)M). \| 1M _m=1^M _m _m -C_t \|_op≤ c_3K^2 ( d+ (1/δ)M+ d+ (1/δ)M ). (29) It remains to bound the centering term ‖ξ¯‖22\| ξ\|_2^2. For any fixed v∈d−1v ^d-1, the scalar random variable ⟨ξ¯,v⟩ ξ,v is conditionally sub-Gaussian with parameter at most K/MK/ M. Applying the same 1/41/4-net argument and a union bound gives that, with ℙtP_t-probability at least 1−δ/21-δ/2, ‖ξ¯‖22≤c4K2d+log⁡(1/δ)M,\| ξ\|_2^2≤ c_4K^2 d+ (1/δ)M, (30) for an absolute constant c4>0c_4>0. Combining (27), (29), and (30), and absorbing constants into a universal constant c0>0c_0>0, we obtain that, with conditional probability at least 1−δ1-δ, ‖C^t−Ct⋆‖op≤c0K2(d+log⁡(1/δ)M+d+log⁡(1/δ)M).\| C_t-C_t \|_op≤ c_0K^2 ( d+ (1/δ)M+ d+ (1/δ)M ). Since the argument holds conditionally on ℱt−1F_t-1, the same high-probability bound also holds marginally. This completes the proof. ∎ Discussion on the Tightness of the Covariance Bound. The high-probability bound established in Proposition 1 is statistically tight and cannot be substantially improved without imposing strong structural assumptions on the generative model’s output [7]. Specifically, estimating an arbitrary d×d× d covariance matrix in the operator norm requires Ω⁡(d/M) ( d/M) error, as established by standard minimax lower bounds via Fano’s inequality [7]. In the standard operational regime of Space, where the number of generated samples M is reasonably large relative to the feature dimension d (i.e., M≥dM≥ d), the linear dimension term ⁡(d/M)O(d/M) is negligible, and our estimator achieves the optimal ⁡(d/M)O( d/M) minimax rate. This bound would only be considered loose if the true conditional covariance Ct⋆C_t possessed known, exploitable structure—such as high sparsity or strict low-rank constraints. In such cases, thresholded or PCA-based estimators could replace C^t C_t to achieve rates depending on the effective intrinsic dimension rather than the ambient dimension d. However, because Space operates on black-box generative models where the latent geometry of the predictive sample cloud is not known a priori, the unstructured minimax optimal rate achieved here represents the fundamental statistical limit of the local geometry estimation step. A.3 Proof of Proposition 2 Proof. Fix a prediction time t. We work under the conditional probability measure ℚtQ_t specified in Assumption (A2). Let J=||J=|G|. Define the collection of clean tests clean=(L,,j):L≤Lt⋆,∈twhole(L),ℬtinc(L),j∈1,2,3.T_clean= \(L,A,j):L≤ L_t ,\,A∈\C_t^whole(L),B_t^inc(L)\,\,j∈\1,2,3\ \. By the definition of Lt⋆L_t and the nesting of the backward candidate windows, for every L≤Lt⋆L≤ L_t , twhole(L)∪ℬtinc(L)∪t⊆τt+1,…,t−1.C_t^whole(L) _t^inc(L) _t \ _t+1,…,t-1\. Therefore every test in cleanT_clean is a clean-block test. By the clean-block calibration part of Assumption (A2), ℚt(KSt,L,j>τt,L,j)≤δuniQ_t (KS_t,L,j^A> _t,L,j^A )≤ _uni for every (L,,j)∈clean(L,A,j) _clean. Since |clean|≤6J|T_clean|≤ 6J, the union bound gives ℚt(∃(L,,j)∈clean:KSt,L,j>τt,L,j)≤6Jδuni.Q_t (∃(L,A,j) _clean:KS_t,L,j^A> _t,L,j^A )≤ 6J _uni. Next consider any candidate length L>Lt⋆L>L_t . By the contaminated-window detectability part of Assumption (A2), there exists a tested pair (L,jL),L∈twhole(L),ℬtinc(L),jL∈1,2,3,(A_L,j_L), _L∈\C_t^whole(L),B_t^inc(L)\, j_L∈\1,2,3\, such that ℚt(KSt,L,jLL≤τt,L,jLL)≤βdet.Q_t (KS_t,L,j_L^A_L≤ _t,L,j_L^A_L )≤ _det. The selector accepts candidate length L only if all block-diagnostic checks for that length pass. Hence acceptance of L implies that the detectably contaminated test (L,jL)(A_L,j_L) passes. Therefore, ℚt(Space accepts candidate length L)≤βdet.Q_t ($ Space$ accepts candidate length L )≤ _det. Taking a union bound over all invalid candidate lengths L>Lt⋆L>L_t , of which there are at most J, yields ℚt(∃L>Lt⋆:Space accepts candidate length L)≤Jβdet.Q_t (∃ L>L_t :$ Space$ accepts candidate length L )≤ J _det. Let ℰtE_t be the event that no clean test is falsely rejected and no invalid candidate length is accepted. Combining the two union bounds, ℚt(ℰtc)≤6Jδuni+Jβdet.Q_t(E_t^c)≤ 6J _uni+J _det. On ℰtE_t, every candidate length L≤Lt⋆L≤ L_t passes all required checks, while every candidate length L>Lt⋆L>L_t is not accepted. Since the selector returns the longest accepted candidate length, L^t=Lt⋆on ℰt. L_t=L_t E_t. Therefore, ℚt(L^t=Lt⋆)≥1−6||δuni−||βdet.Q_t( L_t=L_t )≥ 1-6|G| _uni-|G| _det. If the sufficient empirical-process condition in Assumption (A2) holds with βdet=2exp⁡−2|t|Δsep2, _det=2 \-2|P_t| _sep^2\, then ℚt(L^t=Lt⋆)≥1−6||δuni−2||exp⁡−2|t|Δsep2.Q_t( L_t=L_t )≥ 1-6|G| _uni-2|G| \-2|P_t| _sep^2\. This completes the proof. ∎ A.4 Proof of Theorem 1 Lemma 1 (Exchangeability of oracle ideal radial scores). Let ℋt⋆H_t be the oracle same-regime calibration set associated with Lt⋆L_t . Under Assumption (A3), the ideal score collection es⋆:s∈ℋt⋆∪et⋆\e_s :s _t \∪\e_t \ is conditionally exchangeable under ℚtQ_t. Proof. By the definition of Lt⋆L_t , the indices ℋt⋆∪tH_t ∪\t\ lie in the same regime. Assumption (A3) directly gives conditional exchangeability of the ideal score collection for any same-regime index set. Taking ℋt=ℋt⋆H_t=H_t proves the claim. ∎ Proof of Theorem 1. Fix a prediction time t, and let nt⋆=|ℋt⋆|n_t =|H_t |. Let t,or⋆(αt)C_t,or ( _t) denote the oracle same-regime conformal set constructed from ideal scores on ℋt⋆H_t , and let t,or(αt)C_t,or( _t) denote the corresponding oracle set constructed from empirical scores. Let tACI(α)C_t^ACI(α) denote the actual set returned by the algorithm when the nominal target is α, internally using αt _t. Step 1: Oracle ideal conformal coverage. By Lemma 1 and Assumption (A4b), conditional on αt _t, the rank of et⋆e_t among the nt⋆+1n_t +1 oracle ideal scores is uniform. Hence the split-conformal rank argument gives |ℚt(Yt∈t,or⋆(αt))−(1−ℚt[αt])|≤1nt⋆+1. |Q_t\! (Y_t _t,or ( _t) )- (1-E_Q_t[ _t] ) |≤ 1n_t +1. Step 2: Covariance perturbation. Define tcov=maxs∈ℋt⋆∪t∥C^s−Cs⋆∥op≤ηM,δt/(nt⋆+1).G_t^cov= \ _s _t ∪\t\\| C_s-C_s \|_op≤ _M, _t/(n_t +1) \. By Proposition 1, iterated expectation, and a union bound over the nt⋆+1n_t +1 time points, ℚt((tcov)c)≤δt.Q_t\! ((G_t^cov)^c )≤ _t. Since ηM,δt/(nt⋆+1)≤η0 _M, _t/(n_t +1)≤ _0, Assumption (A1c) applies on tcovG_t^cov. Thus replacing ideal scores by empirical scores changes the oracle conformal coverage probability by at most LΣηM,δt/(nt⋆+1)L_ _M, _t/(n_t +1). Accounting for (tcov)c(G_t^cov)^c, |ℚt(Yt∈t,or(αt))−ℚt(Yt∈t,or⋆(αt))|≤LΣηM,δt/(nt⋆+1)+δt. |Q_t\! (Y_t _t,or( _t) )-Q_t\! (Y_t _t,or ( _t) ) |≤ L_ _M, _t/(n_t +1)+ _t. Step 3: Window-selection error. Let ℰt=L^t=Lt⋆E_t=\ L_t=L_t \. By Proposition 2, ℚt(ℰtc)≤ρtwin.Q_t(E_t^c)≤ _t^win. On ℰtE_t, the selected empirical calibration window equals the oracle empirical calibration window, so tACI(α)=t,or(αt).C_t^ACI(α)=C_t,or( _t). Therefore, |ℚt(Yt∈tACI(α))−ℚt(Yt∈t,or(αt))|≤ℚt(ℰtc)≤ρtwin. |Q_t\! (Y_t _t^ACI(α) )-Q_t\! (Y_t _t,or( _t) ) | _t(E_t^c)≤ _t^win. Step 4: ACI tracking. Combining Steps 1–3, |ℚt(Yt∈tACI(α))−(1−ℚt[αt])|≤ |Q_t\! (Y_t _t^ACI(α) )- (1-E_Q_t[ _t] ) |≤ ρtwin+1nt⋆+1 _t^win+ 1n_t +1 +LΣηM,δt/(nt⋆+1)+δt. +L_ _M, _t/(n_t +1)+ _t. By Assumption (A4a), |(1−ℚt[αt])−(1−α)|≤ℚt[|αt−α|]≤εtACI. | (1-E_Q_t[ _t] )-(1-α) | _Q_t[| _t-α|]≤ _t^ACI. The triangle inequality yields |ℚt(Yt∈tACI(α))−(1−α)|≤ρtwin+εtACI+1nt⋆+1+LΣηM,δt/(nt⋆+1)+δt. |Q_t\! (Y_t _t^ACI(α) )-(1-α) |≤ _t^win+ _t^ACI+ 1n_t +1+L_ _M, _t/(n_t +1)+ _t. This completes the proof. ∎ A.5 Proof of Corollary 1 Proof. The empirical average coverage is Cov^T=1T∑t=1TYt∈tACI(α). Cov_T= 1T _t=1^T1\Y_t _t^ACI(α)\. By the tower property, ⁡[Cov^T]=1T∑t=1T⁡[ℚt(Yt∈tACI(α))].E[ Cov_T]= 1T _t=1^TE\! [Q_t\! (Y_t _t^ACI(α) ) ]. Hence |⁡[Cov^T]−(1−α)| |E[ Cov_T]-(1-α) | ≤1T∑t=1T⁡[|ℚt(Yt∈tACI(α))−(1−α)|]. ≤ 1T _t=1^TE\! [ |Q_t\! (Y_t _t^ACI(α) )-(1-α) | ]. Applying Theorem 1 for each t gives |⁡[Cov^T]−(1−α)|≤1T∑t=1T(ρtwin+εtACI+1nt⋆+1+LΣηM,δt/(nt⋆+1)+δt), |E[ Cov_T]-(1-α) |≤ 1T _t=1^T ( _t^win+ _t^ACI+ 1n_t +1+L_ _M, _t/(n_t +1)+ _t ), provided these per-time terms are deterministic or hold almost surely under the outer expectation. If they are random, the same display holds with outer expectations applied to each right-hand-side term. Defining ρ¯Twin=1T∑t=1Tρtwin,ε¯TACI=1T∑t=1TεtACI, ρ_T^win= 1T _t=1^T _t^win, _T^ACI= 1T _t=1^T _t^ACI, η¯T=1T∑t=1TηM,δt/(nt⋆+1),δ¯T=1T∑t=1Tδt, η_T= 1T _t=1^T _M, _t/(n_t +1), δ_T= 1T _t=1^T _t, we obtain |⁡[Cov^T]−(1−α)|≤ρ¯Twin+ε¯TACI+1T∑t=1T1nt⋆+1+LΣη¯T+δ¯T. |E[ Cov_T]-(1-α) |≤ ρ_T^win+ _T^ACI+ 1T _t=1^T 1n_t +1+L_ η_T+ δ_T. This completes the proof. ∎ Appendix B Algorithm Algorithm 1 Sequential Space wrapper 1: Predictive sample stream ^tt=1T\ Y_t\_t=1^T, observations revealed sequentially, target miscoverage α, training scales D, probe length p, window range [Lmin,Lmax][L_ ,L_ ], step size h 2: Initialize internal level αt←α _t←α 3: Initialize score and diagnostic histories 4: for t=tcal,…,Tt=t_cal,…,T do 5: Compute center μt _t from ^t Y_t 6: Form predictive residual samples rt(m)=y^t(m)−μtr_t^(m)= y_t^(m)- _t 7: Standardize r~t(m)=D−1rt(m) r_t^(m)=D^-1r_t^(m) 8: Compute C^t=Cov⁡(r~t(1),…,r~t(M)) C_t=Cov( r_t^(1),…, r_t^(M)) and set Σ^t=DC^tD _t=D C_tD 9: Form the regularized inverse Σ^t+ _t^+ 10: Construct t=t−p,…,t−1P_t=\t-p,…,t-1\ from stored past diagnostics 11: Select LtL_t by backward same-regime search over L∈Lmin,Lmin+h,…,LmaxL∈\L_ ,L_ +h,…,L_ \ 12: Accept the extension from L−hL-h to L only if both twhole(L)C_t^whole(L) and ℬtinc(L)B_t^inc(L) pass the same-regime test. Stop at the first rejected extension and set LtL_t to the last accepted length 13: Form ℋt(Lt)=es:s∈twhole(Lt)∪tH_t(L_t)=\e_s:s _t^whole(L_t) _t\ 14: Compute rt=Q1−αt(ℋt(Lt))r_t=Q_1- _t(H_t(L_t)) 15: Output t(α)=y∈ℝd:(y−μt)⊤Σ^t+(y−μt)≤rt2C_t(α)= \y ^d:(y- _t) _t^+(y- _t)≤ r_t^2 \ 16: Observe yty_t and compute et=(yt−μt)⊤Σ^t+(yt−μt)e_t= (y_t- _t) _t^+(y_t- _t) 17: Compute r~t=D−1(yt−μt) r_t=D^-1(y_t- _t) and diagnostics st,1=et,st,2=‖r~t‖∞,st,3=|v1⊤r~t|s_t,1=e_t, s_t,2=\| r_t\|_∞, s_t,3=|v_1 r_t| 18: Store ete_t and (st,1,st,2,st,3)(s_t,1,s_t,2,s_t,3) 19: Update αt _t using ACI 20: end for Appendix C Result Tables C.1 Performance across forecasters In this subsection, we evaluate Space across heterogeneous probabilistic forecasters. Table 4 averages results over datasets for each forecaster, so it should be read as a robustness check on the intended application of the wrapper. Table 4: Forecaster summary under Space at the 90%90\% target. Forecasters are GP-Copula [30], TimeGrad [27], CSDI [33], TSFlow [14], DeepAR [31], NsDiff [43], TMDM [19], and TACTiS-2 [3]. Coverage columns report average empirical coverage, while Gap columns report the average absolute coverage gaps. Gap reduction is reported as the relative percentage reduction 100×(RawGap−SpaceGap)/RawGap100×( RawGap- SpaceGap)/ RawGap. Backbone Raw Cov Raw Gap ↓ SPACE Cov SPACE Gap ↓ Gap Red. CSDI 0.787 0.125 0.902 0.004 96.8% DeepAR 0.636 0.274 0.897 0.003 99.0% GP-Copula 0.400 0.500 0.898 0.003 99.4% NsDiff 0.824 0.189 0.899 0.003 98.4% TACTiS-2 0.376 0.524 0.903 0.004 99.3% TimeGrad 0.371 0.529 0.897 0.005 99.0% TMDM 0.644 0.308 0.899 0.002 99.4% TSFlow 0.509 0.403 0.901 0.002 99.4% Space improves calibration across a wide range of forecasters with very different raw uncertainty behavior. It substantially reduces the gap for clearly under-covered forecasters such as GP-Copula, TACTiS-2, and TimeGrad while also improving on raw forecasters that are already producing competitive results, such as NsDiff. This is consistent with the role of the wrapper: it is designed to pull realized coverage toward the nominal target regardless of whether the raw samples under- or over-cover. C.2 Full result tables We report the complete numerical results for all evaluated methods at the three target coverage levels in this subsection. These tables complement the aggregate summaries in the main text by providing the full dataset–forecaster level comparisons, including empirical joint coverage gap, and efficiency metrics for each conformal wrapper. Table 5: Detailed conformal wrapper comparison at the 90%90\% target. The table reports coverage gap and log-volume for each dataset–backbone pair across all conformal wrappers. Lower is better for both metrics. Coverage gap ↓ LogVol ↓ Dataset Backbone SPACE MultiDim SPCI Local Ell. Emp. Cop. CopulaCPTS Gauss. Cop. SPACE MultiDim SPCI Local Ell. Emp. Cop. CopulaCPTS Gauss. Cop. electricity DeepAR 0.004 0.139 0.139 0.121 0.083 0.002 1797.94 1743.91 1751.16 2072.63 2816.26 2166.29 GP-Copula 0.002 0.142 0.151 0.254 0.072 0.030 1827.45 1851.09 1842.54 2082.68 2970.49 2186.00 NsDiff 0.000 0.079 0.426 0.585 0.099 0.124 3436.61 1916.57 1879.74 2046.48 3047.53 2201.14 TACTiS-2 0.002 0.134 0.183 0.080 0.081 0.028 1867.36 1790.10 1791.68 2134.99 2848.78 2187.38 TimeGrad 0.003 0.073 0.384 0.580 0.075 0.299 3955.57 2066.06 2034.33 2126.94 3086.63 2236.49 TMDM 0.002 0.065 0.390 0.477 0.043 0.109 3446.15 1986.95 1956.58 2102.96 3140.46 2221.45 TSFlow 0.003 0.117 0.121 0.142 0.077 0.002 2787.87 1758.87 1756.52 2051.49 2803.15 2191.19 etth1 CSDI 0.003 0.007 0.054 0.024 0.099 0.027 10.52 7.53 7.51 11.49 20.79 11.69 DeepAR 0.002 0.015 0.041 0.028 0.058 0.041 9.79 7.76 7.86 11.13 15.68 11.56 GP-Copula 0.005 0.011 0.058 0.047 0.037 0.053 10.26 8.70 9.17 12.08 16.50 12.40 NsDiff 0.005 0.014 0.063 0.023 0.034 0.023 11.90 9.05 9.49 11.81 16.47 11.81 TACTiS-2 0.002 0.005 0.064 0.027 0.034 0.034 10.30 7.66 8.18 11.49 15.11 11.67 TimeGrad 0.008 0.025 0.165 0.033 0.058 0.009 16.19 11.96 13.22 14.35 19.49 14.55 TMDM 0.000 0.002 0.057 0.031 0.060 0.043 11.46 7.67 8.13 11.81 16.61 12.16 TSFlow 0.003 0.012 0.063 0.030 0.100 0.037 10.37 7.87 8.14 11.19 20.87 11.63 ettm1 CSDI 0.000 0.017 0.054 0.030 0.009 0.042 6.59 4.60 5.27 9.06 11.78 9.46 DeepAR 0.001 0.009 0.059 0.028 0.019 0.038 6.74 4.82 5.74 8.80 11.99 9.19 GP-Copula 0.000 0.005 0.049 0.021 0.006 0.038 5.53 5.95 6.59 9.74 12.44 10.21 NsDiff 0.005 0.002 0.046 0.016 0.009 0.029 9.04 7.28 7.57 9.68 13.78 10.00 TACTiS-2 0.001 0.004 0.054 0.027 0.001 0.037 6.55 4.91 5.47 9.14 11.84 9.55 TimeGrad 0.011 0.021 0.117 0.015 0.076 0.007 13.51 9.18 10.30 12.40 18.36 12.55 TMDM 0.001 0.005 0.053 0.026 0.003 0.038 7.26 4.70 5.33 8.94 11.50 9.31 TSFlow 0.000 0.018 0.057 0.025 0.011 0.037 7.22 5.44 5.94 9.16 11.69 9.53 ettm2 CSDI 0.000 0.009 0.019 0.004 0.004 0.019 10.59 6.82 6.33 10.02 12.37 10.33 DeepAR 0.000 0.006 0.036 0.004 0.007 0.010 11.17 7.50 7.68 10.43 12.61 10.64 GP-Copula 0.000 0.007 0.054 0.019 0.019 0.008 9.72 8.26 7.73 10.87 12.22 11.09 NsDiff 0.001 0.001 0.009 0.010 0.009 0.012 12.61 10.50 10.40 12.30 13.55 12.66 TACTiS-2 0.000 0.010 0.009 0.000 0.009 0.017 10.84 6.65 6.33 9.57 12.52 9.85 TimeGrad 0.000 0.010 0.039 0.011 0.074 0.019 14.24 13.97 14.12 15.15 18.35 15.25 TMDM 0.000 0.004 0.012 0.004 0.005 0.015 10.76 7.23 6.91 10.04 12.24 10.28 TSFlow 0.000 0.005 0.009 0.001 0.026 0.007 19.02 15.57 16.04 16.59 19.81 16.55 exchange CSDI 0.014 0.018 0.001 0.041 0.044 0.041 -36.70 -23.14 -35.28 -36.88 -29.04 -36.88 DeepAR 0.009 0.055 0.018 0.140 0.024 0.094 -29.75 -18.29 -28.40 -30.11 -22.54 -29.78 GP-Copula 0.012 0.025 0.002 0.055 0.014 0.055 -34.64 -20.83 -31.97 -32.78 -27.08 -32.78 NsDiff 0.005 0.048 0.044 0.134 0.051 0.134 -29.63 -21.82 -26.09 -29.38 -21.91 -29.38 TACTiS-2 0.014 0.018 0.024 0.048 0.064 0.048 -36.41 -23.43 -34.81 -36.66 -28.35 -36.66 TimeGrad 0.009 0.045 0.551 0.818 0.821 0.818 -11.24 -9.03 -12.69 -17.86 -6.30 -17.86 TMDM 0.009 0.012 0.047 0.439 0.196 0.439 -13.05 -9.87 -15.25 -18.56 -9.52 -18.56 TSFlow 0.005 0.028 0.166 0.117 0.101 0.084 -33.11 -20.75 -33.39 -35.09 -27.04 -34.92 traffic DeepAR 0.001 0.070 0.034 0.367 0.069 0.367 -2004.30 -1655.92 -1917.09 -1110.90 -1444.22 -1110.90 GP-Copula 0.003 0.060 0.018 0.417 0.057 0.417 -2043.70 -1485.45 -1843.90 -1047.90 -1498.27 -1047.90 NsDiff 0.004 0.007 0.058 0.503 0.027 0.503 -1666.08 -1049.60 -1774.45 -1090.52 -1513.42 -1090.52 TACTiS-2 0.006 0.054 0.023 0.419 0.057 0.419 -1947.96 -1612.59 -1837.45 -1075.69 -1471.75 -1075.69 TimeGrad 0.006 0.070 0.078 0.659 0.039 0.659 -1419.68 -1282.09 -1501.35 -955.67 -1165.37 -955.67 TMDM 0.000 0.020 0.035 0.436 0.015 0.436 -1524.18 -704.16 -1708.73 -1086.10 -1484.19 -1086.10 TSFlow 0.005 0.041 0.055 0.421 0.059 0.421 -2008.51 -1600.29 -1888.35 -1111.23 -1576.92 -1111.23 weather CSDI 0.002 0.007 0.059 0.115 0.021 0.045 41.87 5.26 2.86 25.90 100.70 28.92 DeepAR 0.003 0.004 0.036 0.157 0.027 0.125 54.96 35.45 39.73 49.29 113.17 50.42 GP-Copula 0.001 0.009 0.057 0.107 0.005 0.049 43.25 22.63 21.91 41.69 102.62 44.55 NsDiff 0.001 0.003 0.020 0.096 0.024 0.058 61.59 53.72 55.02 59.86 108.53 61.17 TACTiS-2 0.001 0.000 0.020 0.051 0.020 0.045 43.82 8.22 4.96 28.99 100.33 29.32 TimeGrad 0.001 0.001 0.085 0.271 0.010 0.229 57.09 41.14 42.66 46.96 106.63 48.13 TMDM 0.001 0.004 0.001 0.127 0.022 0.068 50.32 35.10 36.39 43.23 102.61 45.18 TSFlow 0.000 0.170 0.034 0.264 0.026 0.221 60.52 67.85 79.83 43.48 104.69 44.51 Table 6: Detailed conformal wrapper comparison at the 50%50\% target. The table reports coverage gap and log-volume for each dataset–backbone pair across all conformal wrappers. Lower is better for both metrics. Coverage gap ↓ LogVol ↓ Dataset Backbone SPACE MultiDim SPCI Local Ell. Emp. Cop. CopulaCPTS Gauss. Cop. SPACE MultiDim SPCI Local Ell. Emp. Cop. CopulaCPTS Gauss. Cop. electricity DeepAR 0.004 0.111 0.122 0.197 0.458 0.121 1716.61 1696.69 1684.79 1884.26 2702.49 1914.07 GP-Copula 0.001 0.100 0.158 0.158 0.226 0.088 1765.49 1802.28 1766.75 1990.40 2564.03 2017.31 NsDiff 0.002 0.115 0.305 0.385 0.378 0.287 2694.55 1876.96 1834.83 1972.64 2799.25 2010.21 TACTiS-2 0.005 0.106 0.179 0.140 0.159 0.007 1776.44 1743.45 1725.51 1964.06 2500.17 1993.67 TimeGrad 0.001 0.108 0.232 0.373 0.175 0.337 3304.03 2042.12 1995.94 2076.98 2819.30 2088.17 TMDM 0.002 0.106 0.262 0.331 0.095 0.195 2759.85 1960.32 1928.50 2040.14 2900.43 2071.59 TSFlow 0.001 0.105 0.103 0.193 0.424 0.098 2351.84 1714.78 1695.45 1907.09 2643.04 1933.06 etth1 CSDI 0.001 0.012 0.160 0.065 0.499 0.105 4.63 4.16 3.79 5.78 20.79 6.30 DeepAR 0.005 0.038 0.216 0.081 0.458 0.147 4.39 4.33 4.43 5.69 15.68 6.43 GP-Copula 0.008 0.027 0.189 0.095 0.047 0.144 5.33 5.58 5.54 7.15 8.72 7.77 NsDiff 0.004 0.011 0.174 0.077 0.199 0.115 6.18 6.02 6.06 6.95 12.25 7.45 TACTiS-2 0.005 0.001 0.183 0.113 0.094 0.156 4.68 4.18 4.18 6.26 8.28 6.72 TimeGrad 0.002 0.038 0.196 0.110 0.160 0.027 13.35 10.41 10.36 10.60 16.19 11.24 TMDM 0.008 0.006 0.128 0.060 0.065 0.115 6.05 4.85 4.77 6.77 9.49 7.26 TSFlow 0.006 0.001 0.183 0.081 0.500 0.140 4.70 4.44 4.33 5.98 20.87 6.58 ettm1 CSDI 0.001 0.018 0.152 0.078 0.037 0.116 1.37 1.05 1.07 2.25 3.66 2.64 DeepAR 0.001 0.015 0.175 0.100 0.057 0.135 1.48 1.09 1.49 2.42 3.67 2.82 GP-Copula 0.001 0.020 0.152 0.061 0.047 0.091 1.04 2.51 2.22 3.26 4.51 3.64 NsDiff 0.000 0.004 0.110 0.016 0.004 0.074 4.41 4.71 4.35 4.85 7.80 5.30 TACTiS-2 0.000 0.000 0.150 0.062 0.040 0.119 1.31 1.31 1.30 2.26 3.78 2.82 TimeGrad 0.003 0.063 0.204 0.113 0.240 0.043 10.58 7.31 6.70 8.15 14.17 8.80 TMDM 0.000 0.019 0.188 0.074 0.022 0.129 1.87 1.29 1.65 2.45 3.68 3.01 TSFlow 0.001 0.033 0.145 0.067 0.050 0.107 2.18 1.91 2.11 2.89 4.18 3.27 ettm2 CSDI 0.001 0.002 0.010 0.027 0.049 0.069 5.15 4.03 2.89 5.45 6.89 5.98 DeepAR 0.002 0.008 0.062 0.002 0.023 0.039 6.07 5.62 4.46 6.61 7.16 7.03 GP-Copula 0.004 0.009 0.042 0.003 0.000 0.039 4.97 5.18 3.86 6.71 7.17 7.21 NsDiff 0.000 0.004 0.055 0.029 0.006 0.025 8.95 8.84 7.66 8.99 9.82 9.34 TACTiS-2 0.001 0.001 0.007 0.013 0.047 0.069 4.64 4.19 2.62 4.16 6.81 5.09 TimeGrad 0.003 0.020 0.120 0.062 0.201 0.081 12.10 12.54 12.20 12.41 14.85 12.55 TMDM 0.002 0.005 0.015 0.019 0.028 0.047 6.09 5.13 4.01 6.17 6.92 6.41 TSFlow 0.002 0.003 0.041 0.001 0.024 0.016 13.64 14.05 13.38 13.63 16.48 13.75 exchange CSDI 0.000 0.026 0.036 0.128 0.326 0.102 -41.58 -25.69 -39.72 -42.43 -31.77 -42.04 DeepAR 0.016 0.023 0.082 0.049 0.411 0.003 -33.27 -19.66 -32.02 -33.12 -23.03 -32.77 GP-Copula 0.010 0.010 0.122 0.003 0.345 0.023 -39.91 -23.05 -36.72 -38.24 -29.11 -37.83 NsDiff 0.010 0.026 0.013 0.086 0.434 0.013 -32.50 -23.40 -30.04 -32.25 -23.10 -31.88 TACTiS-2 0.016 0.046 0.118 0.086 0.464 0.056 -42.47 -25.88 -39.44 -42.62 -28.35 -42.14 TimeGrad 0.003 0.066 0.490 0.500 0.480 0.500 -12.22 -9.98 -14.18 -19.46 -6.49 -19.43 TMDM 0.023 0.023 0.046 0.230 0.204 0.227 -13.76 -10.77 -16.31 -19.98 -9.52 -19.94 TSFlow 0.007 0.049 0.283 0.266 0.214 0.161 -36.90 -22.61 -37.73 -39.89 -27.48 -38.94 traffic DeepAR 0.001 0.078 0.094 0.205 0.111 0.174 -2751.54 -2041.51 -2780.99 -1496.57 -2275.99 -1413.46 GP-Copula 0.001 0.090 0.094 0.226 0.209 0.188 -2763.05 -1871.26 -2697.83 -1431.22 -2037.54 -1349.32 NsDiff 0.001 0.053 0.117 0.359 0.117 0.103 -2252.30 -1306.73 -2347.17 -1488.17 -2020.75 -1090.52 TACTiS-2 0.003 0.090 0.110 0.232 0.115 0.187 -2633.07 -2013.00 -2766.84 -1456.24 -2197.27 -1373.58 TimeGrad 0.006 0.093 0.127 0.369 0.323 0.356 -1977.18 -1547.09 -2157.70 -1253.83 -1439.78 -1193.78 TMDM 0.000 0.001 0.075 0.315 0.321 0.036 -2036.92 -905.44 -2123.10 -1467.98 -1649.90 -1086.10 TSFlow 0.011 0.067 0.103 0.271 0.125 0.238 -2693.86 -1980.92 -2595.17 -1519.05 -2249.78 -1434.20 weather CSDI 0.004 0.021 0.178 0.144 0.074 0.044 24.26 -2.85 -10.86 -3.75 53.21 5.73 DeepAR 0.004 0.014 0.086 0.179 0.045 0.114 37.54 30.15 30.30 32.18 67.44 35.49 GP-Copula 0.004 0.004 0.114 0.139 0.074 0.052 26.74 13.03 7.31 -34.41 58.39 0.05 NsDiff 0.004 0.009 0.039 0.081 0.367 0.015 49.90 49.44 48.36 50.91 104.03 52.31 TACTiS-2 0.002 0.001 0.155 0.158 0.070 0.074 23.38 -0.54 -8.53 0.08 53.40 6.62 TimeGrad 0.004 0.019 0.016 0.143 0.035 0.103 41.73 37.48 35.29 35.94 81.61 37.35 TMDM 0.002 0.040 0.014 0.133 0.072 0.057 34.83 30.26 30.21 32.24 77.09 34.27 TSFlow 0.001 0.090 0.031 0.284 0.055 0.225 47.89 66.87 76.82 25.99 71.28 29.01 Table 7: Detailed conformal wrapper comparison at the 95%95\% target. The table reports coverage gap and log-volume for each dataset–backbone pair across all conformal wrappers. Lower is better for both metrics. Coverage gap ↓ LogVol ↓ Dataset Backbone SPACE MultiDim SPCI Local Ell. Emp. Cop. CopulaCPTS Gauss. Cop. SPACE MultiDim SPCI Local Ell. Emp. Cop. CopulaCPTS Gauss. Cop. electricity DeepAR 0.004 0.111 0.068 0.171 0.042 0.052 1823.07 1759.72 1790.24 2072.63 2949.63 2166.29 GP-Copula 0.000 0.116 0.098 0.304 0.040 0.080 1854.07 1866.89 1864.68 2082.68 3051.98 2186.00 NsDiff 0.000 0.057 0.369 0.635 0.050 0.174 3683.00 1928.18 1895.89 2046.48 3144.25 2201.14 TACTiS-2 0.003 0.121 0.110 0.130 0.042 0.078 1903.58 1805.75 1821.25 2134.99 2957.35 2187.38 TimeGrad 0.003 0.067 0.402 0.630 0.044 0.349 4202.56 2073.73 2041.09 2126.94 3188.55 2236.49 TMDM 0.002 0.065 0.376 0.527 0.037 0.159 3757.88 1997.86 1963.54 2102.96 3220.00 2221.45 TSFlow 0.001 0.110 0.085 0.192 0.040 0.052 2967.57 1777.25 1781.26 2051.49 2915.85 2191.19 etth1 CSDI 0.003 0.000 0.034 0.011 0.049 0.021 12.96 8.25 9.33 12.86 20.79 13.36 DeepAR 0.000 0.008 0.024 0.011 0.030 0.024 11.97 8.60 9.25 12.45 17.49 13.10 GP-Copula 0.005 0.000 0.030 0.026 0.024 0.031 12.31 9.45 10.60 13.08 18.69 13.46 NsDiff 0.000 0.010 0.027 0.014 0.030 0.024 13.53 9.95 10.56 12.99 18.42 13.68 TACTiS-2 0.000 0.003 0.033 0.018 0.028 0.026 12.75 8.58 9.64 12.75 17.76 13.28 TimeGrad 0.003 0.010 0.109 0.038 0.010 0.029 16.81 12.42 13.86 14.95 20.49 15.13 TMDM 0.002 0.001 0.028 0.013 0.028 0.018 13.72 8.43 9.43 12.84 18.38 13.45 TSFlow 0.000 0.000 0.030 0.020 0.050 0.024 13.29 8.77 9.61 12.57 20.87 12.97 ettm1 CSDI 0.000 0.002 0.027 0.013 0.007 0.024 9.32 5.91 6.57 10.49 14.63 11.21 DeepAR 0.000 0.002 0.028 0.009 0.002 0.022 9.18 5.92 7.19 10.18 14.06 10.97 GP-Copula 0.001 0.003 0.028 0.024 0.005 0.024 7.55 7.00 8.26 11.77 14.89 11.77 NsDiff 0.007 0.003 0.023 0.009 0.006 0.013 11.44 8.33 8.43 11.27 15.32 11.47 TACTiS-2 0.000 0.006 0.026 0.009 0.001 0.019 9.45 6.18 7.07 10.58 14.39 11.32 TimeGrad 0.019 0.004 0.067 0.026 0.038 0.016 14.26 9.89 10.95 13.45 19.58 13.76 TMDM 0.000 0.003 0.025 0.008 0.002 0.023 9.51 5.73 6.72 10.31 14.23 11.10 TSFlow 0.000 0.004 0.032 0.013 0.000 0.025 9.55 6.53 7.45 10.49 14.21 11.31 ettm2 CSDI 0.001 0.005 0.016 0.007 0.003 0.008 11.95 7.65 7.36 10.87 13.51 11.34 DeepAR 0.000 0.001 0.027 0.011 0.003 0.003 12.64 8.16 8.62 11.33 13.79 11.75 GP-Copula 0.000 0.002 0.029 0.016 0.011 0.006 11.36 8.83 8.88 11.97 13.67 12.87 NsDiff 0.002 0.005 0.005 0.002 0.011 0.004 13.81 11.13 11.17 13.53 14.49 13.59 TACTiS-2 0.001 0.009 0.020 0.004 0.011 0.009 12.34 7.39 7.26 10.87 13.79 11.27 TimeGrad 0.002 0.001 0.024 0.013 0.037 0.015 14.95 14.38 14.84 16.01 19.11 16.09 TMDM 0.000 0.002 0.011 0.009 0.002 0.004 12.25 7.96 7.72 10.89 13.42 11.29 TSFlow 0.001 0.002 0.002 0.008 0.004 0.005 20.47 16.18 16.64 17.26 20.82 17.31 exchange CSDI 0.004 0.026 0.004 0.003 0.027 0.003 -34.50 -22.56 -33.74 -35.03 -26.83 -35.03 DeepAR 0.012 0.065 0.006 0.114 0.009 0.114 -28.47 -17.98 -27.50 -29.07 -21.76 -29.07 GP-Copula 0.004 0.042 0.016 0.088 0.026 0.088 -32.85 -20.47 -30.81 -32.31 -26.16 -32.31 NsDiff 0.004 0.042 0.011 0.128 0.011 0.091 -28.63 -21.37 -25.40 -28.80 -21.38 -28.37 TACTiS-2 0.004 0.036 0.020 0.006 0.024 0.006 -34.39 -22.68 -33.41 -34.39 -27.13 -34.39 TimeGrad 0.009 0.026 0.473 0.845 0.838 0.845 -10.88 -8.83 -12.23 -17.64 -6.17 -17.64 TMDM 0.009 0.026 0.034 0.374 0.246 0.364 -12.82 -9.70 -14.77 -18.17 -9.52 -18.01 TSFlow 0.003 0.039 0.101 0.118 0.095 0.078 -31.84 -20.47 -32.35 -34.26 -26.32 -33.38 traffic DeepAR 0.003 0.067 0.027 0.417 0.049 0.417 -1857.94 -1533.90 -1786.11 -1110.90 -1279.82 -1110.90 GP-Copula 0.002 0.057 0.002 0.467 0.047 0.467 -1920.81 -1373.97 -1677.17 -1047.90 -1245.06 -1047.90 NsDiff 0.006 0.028 0.034 0.553 0.026 0.553 -1537.94 -968.88 -1653.99 -1090.52 -1389.24 -1090.52 TACTiS-2 0.004 0.054 0.018 0.469 0.040 0.469 -1766.47 -1485.68 -1695.04 -1075.69 -1308.90 -1075.69 TimeGrad 0.001 0.050 0.033 0.709 0.031 0.709 -1292.86 -1205.61 -1273.84 -955.67 -929.23 -955.67 TMDM 0.003 0.043 0.027 0.486 0.026 0.486 -1414.25 -629.93 -1600.10 -1086.10 -1365.39 -1086.10 TSFlow 0.002 0.035 0.031 0.471 0.039 0.471 -1866.32 -1488.75 -1755.73 -1111.23 -1445.21 -1111.23 weather CSDI 0.001 0.015 0.022 0.065 0.007 0.035 49.05 8.05 8.75 31.59 111.77 33.00 DeepAR 0.008 0.007 0.020 0.137 0.029 0.078 60.59 37.36 42.74 51.66 121.18 54.18 GP-Copula 0.000 0.001 0.007 0.071 0.005 0.040 49.41 25.51 27.73 46.84 112.24 49.05 NsDiff 0.001 0.007 0.002 0.094 0.002 0.068 66.01 55.47 57.71 61.58 113.15 62.93 TACTiS-2 0.001 0.003 0.011 0.060 0.011 0.036 50.75 10.75 9.63 31.75 112.39 33.28 TimeGrad 0.007 0.004 0.105 0.273 0.023 0.241 61.43 42.84 46.47 48.32 114.37 50.19 TMDM 0.001 0.004 0.003 0.107 0.005 0.057 55.34 37.29 38.28 45.55 109.12 47.81 TSFlow 0.001 0.178 0.021 0.211 0.006 0.165 64.55 68.04 80.36 45.53 114.72 46.73 Appendix D Non-stationarity and Performance Gains Analysis In this appendix, we examine whether the performance gains of Space are associated with the severity of time-varying multivariate dependence structure. Our goal is to test whether datasets with stronger and more rapidly changing joint covariance geometry are precisely the settings in which modeling multivariate uncertainty structure is most beneficial. Local covariance drift index. For each dataset, let Yt∈ℝdY_t ^d denote the standardized multivariate target at time t. We quantify local covariance instability by comparing adjacent rolling windows. For a window size w, define WtL=Yt−w,…,Yt−1,WtR=Yt,…,Yt+w−1,W_t^L=\Y_t-w,…,Y_t-1\, W_t^R=\Y_t,…,Y_t+w-1\, (31) with empirical covariance matrices ΣtL _t^L and ΣtR _t^R. Let Σfull _full denote the covariance matrix of the full standardized series, and define a whitening operator A=(Σfull+λI)−1/2,A=( _full+λ I)^-1/2, (32) where λ>0λ>0 is a small ridge term used only for numerical stability. We then compute the local covariance-drift score δt=‖A⁡(ΣtL−ΣtR)A‖Fd. _t= \|A( _t^L- _t^R)A \|_F d. (33) This quantity is dimension-normalized and measures how much the local covariance structure changes from one temporal regime to the next, relative to the dataset’s own global covariance geometry. To obtain a single dataset-level summary, we aggregate the largest drift episodes using the mean of the top decile of δtt\ _t\_t: LCI=Top10Mean⁡(δtt),LCI=Top10Mean (\ _t\_t ), (34) where larger values indicate stronger local covariance drift. In the plots below, datasets are ordered from left to right by decreasing LCILCI. Interpretation. The local covariance drift index is intended to capture a specific form of nonstationarity: instability in the joint second-order geometry of the target. This is directly relevant for Space, whose construction relies on adaptive multivariate covariance structure when forming calibrated uncertainty sets. If joint covariance geometry is both meaningful and time-varying, then methods that explicitly account for such geometry should benefit most. Results and ranking. Using the local covariance drift index defined above, the datasets are ranked as follows: traffic (4.9394.939), weather (3.3263.326), electricity (3.1083.108), etth1 (1.7381.738), exchange (1.4871.487), ettm2 (1.0751.075), and ettm1 (0.8670.867). Thus, the left-hand side of Figure 4 corresponds to datasets with substantially stronger local covariance deformation over time, while the right-hand side contains the smaller and more stationary regimes. In particular, ettm1 has the smallest local covariance drift index, consistent with its comparatively stable dependence structure. Association with coverage-gap improvement. This ranking is informative for understanding where Space helps most. At the pair level, the local covariance drift index is positively associated with coverage-gap improvement, with Pearson correlation r=0.357r=0.357 (p=0.0080p=0.0080) and Spearman correlation ρ=0.281ρ=0.281 (p=0.0396p=0.0396). At the dataset level, where the sample size is necessarily small (n=7n=7), the association remains positive but is less statistically stable, with Pearson r=0.596r=0.596 (p=0.158p=0.158) and Spearman ρ=0.679ρ=0.679 (p=0.0938p=0.0938). Looking directly at the dataset means, the strongest gains appear in more difficult regimes such as electricity, whose local covariance drift index is 3.1083.108 and whose mean coverage-gap improvement is 0.03440.0344, and traffic, whose local covariance drift index is 4.9394.939 with mean improvement 0.02170.0217. By contrast, smaller and more stationary datasets exhibit much smaller gains; for example, ettm1 has the lowest local covariance drift index (0.8670.867) and also one of the smallest mean improvements (0.00350.0035). Overall, the figure supports the interpretation that Space is most beneficial when the target exhibits stronger time-varying multivariate dependence geometry, while performance becomes much closer to the baselines in lower-drift settings. SPACE versus all baselines. To make this trend visible at a dataset level, Figure 4 compares Space against the mean performance of all baseline methods, with datasets ordered by decreasing local covariance drift. The top panel reports mean coverage gap, and the bottom panel reports mean rolling gap. In both panels, points correspond to dataset means across dataset–forecaster pairs, and error bars show mean ±2± 2 standard errors. The raw points in the background show the individual dataset–forecaster pairs. The same qualitative pattern appears in both panels: Space offers the clearest advantage in the high-drift regime, while performance differences narrow substantially in the low-drift regime. Figure 4: Performance of Space versus the mean over all baselines, with datasets ordered from left to right by decreasing local covariance drift index. Top: mean coverage gap. Bottom: mean rolling gap. Error bars are constructed based on two standard errors of the mean across the dataset–forecaster pairs. The left side of the figure corresponds to higher-dimensional and more nonstationary datasets with stronger local covariance drift, where the advantage of Space is largest; the gap narrows for smaller and more stationary datasets such as ettm1. Overall, this appendix analysis is consistent with the view that Space is especially effective in datasets whose multivariate dependence geometry is both strong and temporally unstable. In contrast, when the target is lower-dimensional and closer to stationary covariance structure, Space remains competitive but its advantage over simpler baselines naturally becomes smaller. Appendix E Dynamic Window-Length Diagnostics In this appendix, we provide a visual diagnostic for the dynamic window length selected by Space, denoted by LtL_t. The goal is to illustrate the intended mechanism behind the dynamic rule. When the underlying target process changes regime, the selected memory length is designed to contract so that stale calibration history is downweighted or discarded. When the process becomes more stable again, the selected memory length can recover, allowing the method to reuse longer and more stable calibration windows. To avoid circularity, we detect regime changes directly from the original multivariate target stream yty_t, rather than from the selected LtL_t path or from Space diagnostics. Specifically, for each dataset–forecaster pair at the 90%90\% target level, we align the valid test times, standardize the multivariate target stream, reduce it by PCA to a low-dimensional representation, and then apply the energy-divisive nonparametric multivariate change-point method of 22. This provides a target-side notion of regime change that is external to the dynamic window-selection rule itself. Figure 5 is constructed to visualize three aspects of the dynamic rule: 1. Contraction: whether the selected window length becomes smaller around a detected change point; 2. Recovery: whether the selected window length increases again after the immediate post-change region; 3. Long-term behavior: whether the selected window length tends to grow as the detected regime becomes older and more stable. Each panel in Figure 5 contains three aligned subplots. The top subplot shows the selected dynamic window length LtL_t over valid test times. The middle subplot shows Target PC1, the first principal component of the standardized multivariate target stream, which provides a one-dimensional summary of the dominant variation in the target process. The bottom subplot shows target-side change points detected by the energy-divisive procedure. The intended visual pattern is as follows. Around a vertical change-point marker, a local downward reset in LtL_t indicates that the dynamic rule shortens its calibration memory when the target process changes. Over a more stable segment, LtL_t may increase again, indicating that the method can use a longer historical window when the local regime persists. Taken together, these illustrative plots are consistent with the intended interpretation of the dynamic window-selection mechanism. The selected memory length is not static; rather, it changes over time, contracts around detected target-side disturbances, and in the clearer examples rebuilds during more stable segments. This behavior reflects the design motivation for dynamic LtL_t: Space uses shorter memory when the recent past becomes unreliable and longer memory when the local regime appears stable enough to benefit from a broader calibration history. (a) ettm2/gpcopula (b) weather/tactis2 Figure 5: Illustrative dynamic-L overlays at the 90%90\% target level using energy-divisive target-side change-point detection [22]. In each plot, the top panel shows the selected window length LtL_t, the middle panel shows the first principal component of the standardized valid-test target stream (Target PC1), and the bottom panel shows detected target-side change points. Across these examples, the selected dynamic memory length exhibits active interior adaptation and contraction-and-rebuilding behavior around detected target-side regime changes. Appendix F Implementation details of Space F.1 Metric definitions For a fixed target miscoverage level α, empirical joint coverage is defined as Cov^=1N∑t=1Nyt∈t(α). Cov= 1N _t=1^NI\y_t _t(α)\. (35) The global coverage gap is the absolute deviation from the target coverage level, Gap=|Cov^−(1−α)|.Gap= | Cov-(1-α) |. (36) For rolling diagnostics, let w denote the sliding-window length and let T denote the set of valid rolling-window endpoints. The rolling joint coverage at endpoint t∈t is RC^t=1w∑s=t−w+1tys∈s(α). RC_t= 1w _s=t-w+1^tI\y_s _s(α)\. (37) The corresponding rolling coverage gap is gt=|RC^t−(1−α)|.g_t= | RC_t-(1-α) |. (38) We summarize these local deviations using three rolling metrics. The mean rolling gap is MeanRollingGap=1||∑t∈gt=1||∑t∈|RC^t−(1−α)|.Mean\ Rolling\ Gap= 1|T| _t g_t= 1|T| _t | RC_t-(1-α) |. (39) The P90 rolling gap is the 9090th percentile of the rolling coverage gaps, P90RollingGap=Q0.9(gt:t∈).P90\ Rolling\ Gap=Q_0.9\! (\g_t:t \ ). (40) The fraction of bad windows is the proportion of rolling windows whose gap exceeds a tolerance threshold δ, Frac.BadWindows=1||∑t∈gt>δ.Frac.\ Bad\ Windows= 1|T| _t I\g_t>δ\. (41) In all main experiments, we use w=30w=30 and δ=0.10δ=0.10. For efficiency, we report the mean log-volume of the prediction regions, LogVol=1N∑t=1Nlog⁡Vol⁡(t(α)). LogVol= 1N _t=1^N \! (C_t(α) ). (42) Lower values are better for coverage gap, rolling diagnostics, and LogVol. The size metric should be interpreted together with calibration, since small sets are meaningful only when coverage is close to the target. F.2 Compute resources All wrapper experiments were run on a single local workstation equipped with a 13th Gen Intel(R) Core(TM) i9-13980HX CPU, an NVIDIA GeForce RTX 4080 Laptop GPU and 32GB system memory. The reported benchmark evaluates all dataset–forecaster pairs and all target levels included in the full result table in Section C.2. Since Space is a post-hoc wrapper, the reported time refers to wrapper execution and evaluation on saved predictive sample streams, excluding the original training of the forecasting backbones. The complete Space run took 1h19m38s. MultiDimSPCI was the most computationally intensive among the baselines, requiring 14h58m33s under the same benchmark scope. The remaining conformal baselines were substantially lighter and completed in comparatively small additional time. Additional exploratory runs and preliminary tuning experiments were performed during development, but the above times describe the compute required to reproduce the reported final benchmark from the saved sample streams. F.3 Experimental configuration Unless stated otherwise, all forecasters use context window 9696, and one-step-ahead prediction (pred_len=1=1) to simulate one-step-ahead forecasting. The saved predictive sample streams have M=100M=100 forecast samples per time point. Note that, in principle, increasing the number of forecast samples reduces sample-geometry error. In the main benchmark, Space is evaluated under a single fixed configuration. Specifically, the candidate calibration lengths are L∈20,30,…,300,L∈\20,30,…,300\, (43) so that Lmin=20L_ =20, Lmax=300L_ =300, and the backward extension step size is h=10h=10. The probe-block length is set to p=20p=20, and the adaptive conformal inference update uses step size η=0.01η=0.01. For numerical robustness, the benchmark implementation uses a shrinkage-stabilized version of the sample-derived covariance estimator, with identity shrinkage parameter λshrinkage=0.30 _shrinkage=0.30. This shrinkage choice is an empirical implementation device introduced to improve stability of the local covariance estimate in practice, especially when the effective sample size is limited or the sample cloud is ill-conditioned. It is not a separate conceptual component of Space, and it is not required by the theoretical results in Section 4. Our theory is stated at the level of the sample-derived local geometry and its estimation error, so the guarantees do not rely on this particular shrinkage form or on the specific value of ρ. For dynamic same-regime calibration window selection, Space uses a DKW-style KS acceptance threshold to instantiate the same-regime test described in Section 3.2. At time t, the threshold is τt(cτ)=cτlog⁡(2KtJδτ)2p, _t(c_τ)=c_τ \! ( 2K_tJ _τ )2p, (44) where cτ=2.0c_τ=2.0 is a multiplicative tuning constant, KtK_t is the number of feasible candidate window lengths at time t, J=3J=3 is the number of diagnostics in the test bundle, δτ=0.05 _τ=0.05, and p=20p=20 is the probe-block length. The factor 2KtJ2K_tJ accounts for the two block comparisons, the feasible candidate lengths, and the diagnostic bundle. This threshold is the empirical implementation of the KS-style same-regime test; the theoretical analysis abstracts its behavior through Assumption (A2), which controls false rejection on clean blocks and false acceptance on contaminated extensions. Appendix G Dataset Details We use the following seven publicly available datasets for academic research and evaluation purposes. The details of the datasets are listed in Table 8. Table 8: Dataset summary for the saved predictive sample streams used in wrapper evaluation. Here d is the target dimension, T is the saved stream length, and Test denotes the Space/baseline wrapper test split. Dataset Freq. d T Test electricity H 321 5261 1053 ETTh1 H 7 3484 697 ETTm1,2 15 min 7 13936 2788 exchange D 8 1518 304 traffic H 862 3509 702 weather 10 min 21 10539 2108 1. Electricity. The dataset records electricity consumption for 321 customers and is used as a standard multivariate forecasting benchmark. URL: https://archive.ics.uci.edu/dataset/321/electricityloaddiagrams20112014. 2. ETTh1, ETTm1, and ETTm2. The datasets contain electricity transformer variables, including load and oil temperature, recorded at hourly or 15-minute resolution [45]. URL: https://github.com/zhouhaoyi/ETDataset. 3. Exchange. The dataset records daily exchange rates of eight countries and is commonly used in multivariate time-series forecasting benchmarks [16]. URL: https://github.com/laiguokun/multivariate-time-series-data. 4. Traffic. The dataset contains hourly road occupancy rates measured by sensors on San Francisco Bay Area freeways and is commonly distributed through multivariate time-series benchmark collections [16]. Raw data are associated with the Caltrans Performance Measurement System (PeMS), which requires an account for access. URL: https://pems.dot.ca.gov/; benchmark mirror: https://github.com/laiguokun/multivariate-time-series-data. 5. Weather. The dataset is commonly distributed through the Autoformer/Time-Series-Library benchmark collection [38]. URL: https://github.com/thuml/Autoformer.

Instant research discovery