Model-agnostic predictability quantification of multivariable time series

# Model-agnostic predictability quantification of multivariable time series
Source: [https://arxiv.org/html/2610.33209](https://arxiv.org/html/2610.33209)
Yuxuan ShuVasileios LamposCentre for Artificial IntelligenceDepartment of Computer ScienceUniversity College London\{yuxuan\.shu\.22,v\.lampos\}@ucl\.ac\.uk

###### Abstract

Multivariable time series forecasting plays an important role in various domains, where predictions depend on complex interactions among covariates\. Forecasting models are commonly evaluated across numerous benchmark datasets\. However, forecasting accuracy reflects both model capability and task predictability, making it difficult to characterise the difficulty of a task independently of the model used\. Existing predictability measures often fail to capture nonlinear dependencies or are limited to univariate series\. In this work, we proposeMultivariablePredictabilityScore \(MAPS\), a model\-agnostic metric for measuring the intrinsic predictability of multivariable time series forecasting tasks\.MAPSuses covariate information to assess predictability across different lookback windows and forecasting horizons, and provides a quantitative assessment of the predictive information available for forecasting\. We evaluateMAPSacross synthetic and real\-world forecasting datasets against66baseline predictability measures\. The former confirms that our method closely aligns with the expected reference predictability \(r=0\.963r=0\.963\)\. Empirical results further corroborate thatMAPSyields more reliable projections in the majority of experimental configurations\.

Source code:[github\.com/ClaudiaShu/MAPS](https://github.com/ClaudiaShu/MAPS)

## 1Introduction

Time\-series forecasting plays an important role across a wide range of domains and supports decision\-making across multidisciplinary applications\. Most forecasting problems are inherently multivariable\. A single target variable is predicted from a set of covariates such as sensor information\([Gneiting et al\., 2006](https://arxiv.org/html/2610.33209#bib.bib30);[Inman et al\., 2013](https://arxiv.org/html/2610.33209#bib.bib31)\), economic indicators\([Sims, 1980](https://arxiv.org/html/2610.33209#bib.bib25);[Stock and Watson, 2002b](https://arxiv.org/html/2610.33209#bib.bib26);[Stock and Watson, 2002a](https://arxiv.org/html/2610.33209#bib.bib27)\), or clinical measurements\([Lehman et al\., 2014](https://arxiv.org/html/2610.33209#bib.bib24);[Henderson et al\., 2000](https://arxiv.org/html/2610.33209#bib.bib28);[Wulfsohn and Tsiatis, 1997](https://arxiv.org/html/2610.33209#bib.bib29)\)\. However, forecasting tasks differ considerably in the predictive information they contain\. Some exhibit strong temporal regularities and stable relationships between covariates and targets, whereas others are dominated by noise or weak predictive signals\. This motivates a quantitative assessment of task predictability prior to modelling\. Such an assessment can indicate whether a dataset supports reliable forecasting at a given horizon and provide insight into the performance that may be achievable\.

Prior studies have attempted to characterise predictability in time series\. An early approach\([Sugihara et al\., 1994](https://arxiv.org/html/2610.33209#bib.bib20)\)uses a forecasting method that predicts target values from similar historical patterns and relies on the forecasting accuracy to distinguish structured dynamics from stochastic noise\.[Zadeh et al\. \(2025\)](https://arxiv.org/html/2610.33209#bib.bib1)use mutual information to measure the relationship between inputs and targets and how much of this relationship remains in a model’s prediction errors\.[Beck et al\. \(2025\)](https://arxiv.org/html/2610.33209#bib.bib3)assess predictability by testing whether forecasting models consistently outperform a naïve persistence model, considering significant reductions in forecast error as evidence of predictable structure in the data\. However, these approaches are model\-dependent, meaning that the resulting predictability estimate depends on both the data and the forecasting model rather than reflecting the intrinsic aleatoric property of the prediction task\. Early approaches from information theory deploy entropy to assess predictability\. Approximate entropy\([Pincus, 1991](https://arxiv.org/html/2610.33209#bib.bib4)\)measures the extent to which recurring patterns in a time series remain similar, where a lower value indicates that the series is more predictable\. Permutation entropy\([Bandt and Pompe, 2002](https://arxiv.org/html/2610.33209#bib.bib17)\)uses Shannon entropy to measure the irregularity of changes in a time series, with lower values indicating higher predictability\. However, neither method is formulated for multivariable forecasting tasks or considers how predictability changes across different forecasting horizons\. Alternative approaches utilised the frequency domain\.[Powell and Percival \(1979\)](https://arxiv.org/html/2610.33209#bib.bib18)introduced spectral entropy, which measures how consistently a time series follows repeating cycles, with lower values indicating clearer patterns and greater predictability\. Building on this, ForeCA\([Goerg, 2013](https://arxiv.org/html/2610.33209#bib.bib21)\)extends to multivariate settings by forming weighted sums of variables and identifying those with lower spectral entropy as more predictable\. More recently,[Wang et al\. \(2025a\)](https://arxiv.org/html/2610.33209#bib.bib13)transformed time series into the frequency domain and used distances between their frequency components as an indicator for forecasting model selection\. Despite their computational efficiency, these frequency\-domain approaches either consider a univariate setup\([Powell and Percival, 1979](https://arxiv.org/html/2610.33209#bib.bib18);[Wang et al\., 2025a](https://arxiv.org/html/2610.33209#bib.bib13)\), or do not account for nonlinear dependencies\([Goerg, 2013](https://arxiv.org/html/2610.33209#bib.bib21)\)\.

Some prior work accounts for the forecast horizon when assessing predictability\.[Wang et al\. \(2025b\)](https://arxiv.org/html/2610.33209#bib.bib2)measures how much frequency patterns vary across time\-series segments, including both the observed history and the target to be predicted, with greater variation indicating lower predictability\. Similarly,[Catt \(2026\)](https://arxiv.org/html/2610.33209#bib.bib19)assesses predictability at each forecast horizon using the mutual information between the observed past and the target value, with lower mutual information indicating less predictive information at that horizon\. However, these approaches are designed for univariate time series and do not account for covariate dependencies\. Furthermore, mutual information estimation in high\-dimensional spaces is known to be computationally prohibitive\([Kraskov et al\., 2004](https://arxiv.org/html/2610.33209#bib.bib22)\), with complexity scaling poorly as the number of variables increases\([Gao et al\., 2015](https://arxiv.org/html/2610.33209#bib.bib23)\)\(Figure[S6](https://arxiv.org/html/2610.33209#A3.F6)\)\.

Overall, existing predictability measures either depend on forecasting models, ignore forecasting horizons, or are limited to univariate settings\. This gap is particularly relevant to the evaluation of forecasting models across benchmarks\. Common long\-term forecasting evaluations repeatedly use a small collection of datasets, including ETT, Electricity, Traffic, Weather, Exchange\-Rate, and ILI\([Zeng et al\., 2023](https://arxiv.org/html/2610.33209#bib.bib34)\), while broader archives have mostly focused on collections of univariate series\([Godahewa et al\., 2021](https://arxiv.org/html/2610.33209#bib.bib45)\)\. Recent benchmarks such as ProbTS\([Zhang et al\., 2024](https://arxiv.org/html/2610.33209#bib.bib46)\), GIFT\-Eval\([Aksu et al\., 2024](https://arxiv.org/html/2610.33209#bib.bib47)\), and fev\-bench\([Shchur et al\., 2025](https://arxiv.org/html/2610.33209#bib.bib15)\)have expanded evaluation across domains and multivariate settings\. Nevertheless, benchmark tasks can still differ widely in the predictive information they contain\. Some exhibit strong recurring patterns that can be captured from the target history, whereas others rely on cross\-variable relationships or contain only weak predictive signals\. A model\-independent measure that considers the available covariates, lookback window, and forecast horizon can help interpret performance beyond conventional forecasting metrics\.

Motivated by the aforementioned observations, our main contributions are:

1. 1\.We proposeMultivariablePredictabilityScore \(MAPS\), a model\-agnostic measure that quantifies the predictability of a multivariable forecasting task by assessing whether similar historical samples lead to consistent future target vectors\.MAPSexplicitly accounts for the available covariates, lookback window, and forecasting horizon without requiring a forecasting model\.
2. 2\.We evaluateMAPSon synthetic time series with diverse dynamics and known predictability\.MAPSclosely follows the estimated predictability, achieving an average linear correlation ofr=0\.963r=0\.963\.
3. 3\.We compareMAPSagainst55existing predictability measures on multivariable Weather and Influenza\-Like Illness forecasting tasks across multiple lookback windows and forecasting horizons\.MAPSshows the strongest relationship with empirical forecasting performance, obtaining correlations ofr=0\.961r=0\.961on Weather andr=0\.827r=0\.827on ILI\.111We note that these correlations are measuring different relationships which stem from the different experimental setups for these two forecasting tasks \(see Section[4\.4](https://arxiv.org/html/2610.33209#S4.SS4)\)\.

![Refer to caption](https://arxiv.org/html/2610.33209v1/figure1.png)Figure 1:Overview of the proposed Multivariable Predictability Score \(MAPS\)\.
## 2Predictability of Multivariable Time Series – Task Definition

We focus on quantifying the predictability of datasets for multivariable forecasting in a model\-agnostic fashion\. Specifically, we consider the multivariable forecasting task over a data set𝐐∈ℝV×C\\mathbf\{Q\}\\in\\mathbb\{R\}^\{V\\times C\}whereVVdenotes the number of observations andCCis number of covariates, including11target andC−1C\-1exogenous variables\. From𝐐\\mathbf\{Q\}we extractNNsamples with a fixed stride ofMM\. This forms an input tensor𝐗∈ℝN×L×C\\mathbf\{X\}\\in\\mathbb\{R\}^\{N\\times L\\times C\}, whereLLis the lookback window length\. The corresponding future target values are captured in𝐘∈ℝN×H\\mathbf\{Y\}\\in\\mathbb\{R\}^\{N\\times H\}, whereHHis the forecasting horizon\. Each sample thus consists of an input𝐗t∈ℝL×C\\mathbf\{X\}\_\{t\}\\\!\\in\\\!\\mathbb\{R\}^\{L\\\!\\times\\\!C\}covering the past time steps\{t−L\+1,…,t\}\\\{t\-L\+1,\\dots,t\\\}, and the corresponding target vector𝐲t∈ℝH\\mathbf\{y\}\_\{t\}\\\!\\in\\\!\\mathbb\{R\}^\{H\}covering the target time steps\{t\+1,…,t\+H\}\\\{t\+1,\\dots,t\+H\\\}\. The task is to define a predictability function𝒫⁡\(𝐗→𝐘\)\\mathcal\{P\}\(\\mathbf\{X\}\\rightarrow\\mathbf\{Y\}\)that quantifies the intrinsic predictability of the future target values from the corresponding input windows, independently of any specific model or supervised learning procedure\.

## 3Multivariable Predictability Score \(MAPS\)

An overview of theMAPSmethod is shown in Figure[1](https://arxiv.org/html/2610.33209#S1.F1)\. We first standardise each variable \(or column\) of𝐐\\mathbf\{Q\}over its full time series\. We then partition the standardised time series of each variable intoNNoverlapping windows as described in section[2](https://arxiv.org/html/2610.33209#S2)\.𝐗t\\mathbf\{X\}\_\{t\}and𝐲t\\mathbf\{y\}\_\{t\}are additionally standardised to maintain local patterns within each sample\.

Dimensionality Reduction with Principal Component Analysis \(PCA\)\.Following standardisation, each input window𝐗t\\mathbf\{X\}\_\{t\}is reshaped into vector𝐳t∈ℝD\\mathbf\{z\}\_\{t\}\\in\\mathbb\{R\}^\{D\}, whereD=L×CD=L\\times C\. This forms the input state matrix held in𝐙∈ℝN×D\\mathbf\{Z\}\\in\\mathbb\{R\}^\{N\\times D\}\. We note that the subsequent predictability computation operates on thisDD\-dimensional representation and does not explicitly consider temporal adjacency\. For multivariable datasets with many covariates or long lookback windows,DDcan be rather large\. To reduce the dimensionality of this representation, we optionally apply PCA\([Jolliffe, 1986](https://arxiv.org/html/2610.33209#bib.bib16)\)to project𝐙\\mathbf\{Z\}ontod=min⁡\(d,N−1,D\)d=\\min\(d,N\-1,D\)leading principal components, whereddis an adjustable target dimensionality\. For simplicity, we use𝐙∈ℝN×d\\mathbf\{Z\}\\in\\mathbb\{R\}^\{N\\times d\}and𝐳t∈ℝd\\mathbf\{z\}\_\{t\}\\in\\mathbb\{R\}^\{d\}to denote the matrix and state vectors after this optional projection, respectively, withd=Dd=Dwhen PCA is not applied\.

Identify Similar States with Nearest Neighbours\.We usekk\-nearest neighbours \(kk\-NN\) to identify similar samples among theNNsliced windows𝐙\\mathbf\{Z\}\. For each input sample𝐳t\\mathbf\{z\}\_\{t\}, we identify itskknearest neighbours𝒩t\\mathcal\{N\}\_\{t\}using the Euclidean distance, excluding𝐳t\\mathbf\{z\}\_\{t\}itself\. The neighbourhood sizekkdetermines a local component for the proposed predictability score\. It should be large enough to provide statistically stable variance estimates, while small enough to preserve the local nature of the analysis\. Following established practices in nearest\-neighbour\([Stone, 1977](https://arxiv.org/html/2610.33209#bib.bib10);[Devroye et al\., 1994](https://arxiv.org/html/2610.33209#bib.bib11);[Levada et al\., 2024](https://arxiv.org/html/2610.33209#bib.bib9)\), we setk=⌊log2⁡\(N\)⌋k=\\lfloor\\log\_\{2\}\(N\)\\rfloor\.

### 3\.1Quantifying Dataset Predictability

Our method is motivated by local regularity in the relationship between observed states and their future evolution\([Stone, 1977](https://arxiv.org/html/2610.33209#bib.bib10)\)\. In deterministic dynamical systems, the current state determines the future trajectory\([Lorenz, 1963](https://arxiv.org/html/2610.33209#bib.bib6)\), while state\-space reconstruction provides a basis for studying how nearby observations evolve\([Takens, 1981](https://arxiv.org/html/2610.33209#bib.bib8);[Eckmann and Ruelle, 1985](https://arxiv.org/html/2610.33209#bib.bib7)\)\. Accordingly, when the input contains predictive information about the target, nearby samples should correspond to target vectors with similar directions\. We thus assess how consistently neighbouring samples align in the vector direction as an empirical measure of the predictive structure captured by the input samples\.

Given the standardised target vectors𝐲1,…,𝐲N∈ℝH\{\\mathbf\{y\}\_\{1\},\\ldots,\\mathbf\{y\}\_\{N\}\}\\in\\mathbb\{R\}^\{H\}, we first obtain the corresponding unit direction vectors as𝐲^t=𝐲t∥𝐲t∥\+ϵ\\hat\{\\mathbf\{y\}\}\_\{t\}=\\frac\{\\mathbf\{y\}\_\{t\}\}\{\\lVert\\mathbf\{y\}\_\{t\}\\rVert\+\\epsilon\}, whereϵ=10−8\\epsilon=10^\{\-8\}is a small constant that prevents division by zero\. The resulting target vectors thus have unit length\. The neighbourhood coherence between the current windowttand a neighbouring windowjjis then obtained by computing the cosine similarity⟨𝐲^t,𝐲^j⟩∈\[−1,1\]\\langle\\hat\{\\mathbf\{y\}\}\_\{t\},\\hat\{\\mathbf\{y\}\}\_\{j\}\\rangle\\in\[\-1,1\]\. Averaging over𝒩t\\mathcal\{N\}\_\{t\}, we have:

ct=1\|𝒩t\|​∑j∈𝒩t⟨𝐲^t,𝐲^j⟩,c\_\{t\}=\\frac\{1\}\{\|\\mathcal\{N\}\_\{t\}\|\}\\sum\_\{j\\in\\mathcal\{N\}\_\{t\}\}\\left\\langle\\hat\{\\mathbf\{y\}\}\_\{t\},\\hat\{\\mathbf\{y\}\}\_\{j\}\\right\\rangle\\,,\(1\)wherect∈\[−1,1\]c\_\{t\}\\in\[\-1,1\]measures the average directional agreement between windowttand the futures of windows in𝒩t\\mathcal\{N\}\_\{t\}, whose observed histories are most similar\.

The predictability of a time series at windowttis then given byst=\|ct\|s\_\{t\}=\\lvert c\_\{t\}\\rvert, where, by construction,st∈\[0,1\]s\_\{t\}\\in\[0,1\]\. Importantly, we take the absolute value ofctc\_\{t\}, after averaging the directional similarities\. This is because both positive \(ct→1c\_\{t\}\\rightarrow 1\) and negative \(ct→−1c\_\{t\}\\rightarrow\-1\) directional agreement of the target vectors are equally indicative of predictability\. In this case, large\|ct\|\\lvert c\_\{t\}\\rvertvalues indicate that thekk\-nearest neighbours of windowttfollow coherent future directions, reflecting a more predictable mapping between similar historical samples and their corresponding target vectors\. By contrast, neighbourhoods containing orthogonal future directions yieldct≈0c\_\{t\}\\approx 0, indicating that windows with similar observed histories do not evolve consistently\. In such cases, the input representation provides limited information about future dynamics\.

The overall predictability score of a time series dataset w\.r\.t\.the specific look\-back window and forecasting horizon is thus given by:

𝒫⁡\(𝐗→𝐲\)=1N​∑t=1Nst∈\[0,1\]\.\\mathcal\{P\}\(\\mathbf\{X\}\\rightarrow\\mathbf\{y\}\)=\\frac\{1\}\{N\}\\sum\_\{t=1\}^\{N\}s\_\{t\}\\in\[0,1\]\\,\\,\\,\.\(2\)Intuitively,𝒫⁡\(𝐗→𝐲\)≈1\\mathcal\{P\}\(\\mathbf\{X\}\\rightarrow\\mathbf\{y\}\)\\approx 1indicates high predictability, where samples that are close in the input space have futures that also evolve similarly, indicating a more deterministic input\-to\-target mapping\. Contrarily,𝒫⁡\(𝐗→𝐲\)≈0\\mathcal\{P\}\(\\mathbf\{X\}\\rightarrow\\mathbf\{y\}\)\\approx 0indicates low predictability, i\.e\.even nearly identical samples result in diverged directions, suggesting chaotic dynamics that are potentially less predictable\.

Computation cost\.Performing PCA costO⁡\(N​D​min⁡\{N,D\}\)O\(ND\\min\\\{N,D\\\}\), an exhaustive Euclidean nearest\-neighbour search overNNwindows costsO⁡\(N2​d\)O\(N^\{2\}d\), and conducting cosine similarity over the neighbouring target vectors costsO⁡\(N​k​H\)O\(NkH\)\. GivenddandHHare fixed andk=⌊log2⁡N⌋k=\\lfloor\\log\_\{2\}N\\rfloor, this complexity simplifies toO⁡\(N​D​min⁡\{N,D\}\+N2\)O\(ND\\min\\\{N,D\\\}\+N^\{2\}\)\. We additionally provide runtime comparison in Appendix[C\.6](https://arxiv.org/html/2610.33209#A3.SS6)\.

Constant target variable\.A special case emerges when the target variable is a constant\. Standardisation is applied over the entire sequence, so a constant target will become equal to zero in every window \(𝐲t=0\\mathbf\{y\}\_\{t\}=0\) and leave the direction undefined\. Given that a constant target can be forecast exactly by its own value, we setst=1s\_\{t\}=1regardless\. A near\-constant target is otherwise scored normally as it does not collapse to zero vector\. We further discuss this in Appendix[C\.1](https://arxiv.org/html/2610.33209#A3.SS1)\.

## 4Results

We evaluateMAPSon both synthetic and real\-world datasets and compare it with66baseline methods\. We further examine its behaviour under distribution shift and robustness to different configurations\.

### 4\.1Multivariable Time Series Datasets

Synthetic time series\.To evaluate our method under controlled conditions, we construct synthetic datasets with controllable predictability\. We begin by combining signals with known deterministic and stochastic components, including sinusoidal, linear trends, exponential and logarithmic growth, and Gaussian noise\. By varying the proportion and magnitude of each component, we control the level of predictability\. In this way, we form66basic synthetic data sets\. We further consider66complex dynamic systems to enable a more comprehensive evaluation: a linear Gaussian state\-space model, a smooth Gaussian process, a fractional Gaussian noise process, a two\-state hidden Markov model, a noisy Lorenz system, and a quasi\-periodic process\. These cover a range of different predictability decay behaviours w\.r\.t\.the forecasting horizon, enabling a horizon\-dependent evaluation of our predictability metric\([Dong et al\., 2025](https://arxiv.org/html/2610.33209#bib.bib49)\)\.

Real\-world tasks\.We form the benchmark multivariable time series datasets from the Weather \(WEA\) and Influenza\-Like Illness \(ILI\) forecasting tasks introduced by[Shu and Lampos \(2026\)](https://arxiv.org/html/2610.33209#bib.bib14)\. ForWEA, the forecasting target is variable ‘t850’ at a target location \(for2020locations in total\), and the time step is set to66hours\. We consider both a task that uses55covariates at the target location, and also expand this to a task with4545covariates across a3×33\\times 3grid of locations \(target and 8 neighbouring ones\)\.ILIcontains daily ILI rates and corresponding Google search frequencies for related queries for USA HHS Regions 2, 9, and 10, and England, UK\.

Further details for both synthetic and real\-world datasets are provided in Appendix[A](https://arxiv.org/html/2610.33209#A1)\.

### 4\.2Experiment configurations

Baseline methods\.We compareMAPSagainst55model\-agnostic baseline predictability methods\. These are a persistence baseline, ForeCA\([Goerg, 2013](https://arxiv.org/html/2610.33209#bib.bib21)\), Spectral Entropy\([Wang et al\., 2025a](https://arxiv.org/html/2610.33209#bib.bib13)\), Pattern Complexity \(PC\)\([Wang et al\., 2025b](https://arxiv.org/html/2610.33209#bib.bib2)\), and Mutual Information \(MI\)\. To adapt baselines that are originally designed for univariate forecasting to multivariable settings, we compute the scores for each input variable w\.r\.t\.the target and average across variables\. We additionally include Epiplexity\([Finzi et al\., 2026](https://arxiv.org/html/2610.33209#bib.bib12)\), a model\-dependent measurement, where we use ridge regression as the predictive model in our experiments\. Further details are provided in Appendix[B](https://arxiv.org/html/2610.33209#A2)\.

Experiment setup\.We conduct experiments with varying lookback window lengths \(LL\) and forecasting horizons \(HH\)\. Specifically, we considerL∈\{28,56,84\}L\\\!\\in\\\!\\\{28,56,84\\\},H∈\{4,12,28\}H\\\!\\in\\\!\\\{4,12,28\\\}forWEA, andL∈\{28,56\}L\\\!\\in\\\!\\\{28,56\\\},H∈\{7,14,21,28\}H\\\!\\in\\\!\\\{7,14,21,28\\\}forILI\. We obtaind=50d=50principal components and set the neighbourhood size tok=⌊log2⁡N⌋k=\\lfloor\\log\_\{2\}N\\rfloor\. Unless otherwise stated, we sample rolling windows using a stride ofM=10M=10\.

Evaluation metrics\.We use the linear correlation coefficientrrbetween estimated and reference predictability scores as the primary metric in our evaluation\. For synthetic data, the reference value is the designed predictability of each data\-generating process\. For real\-world data, we use the performance of a competitive machine learning model as the reference\. Machine learning performance is mainly determined by the linear correlation coefficientρ\\rhobetween the model predictions and observed targets, a metric that is invariant to both scaling and translation\. We additionally consider a scaled absolute forecasting error as a complementary measure\. Specifically, leteedenote the MAE of the model on the test set, andσy\\sigma\_\{y\}be the standard deviation of the target series over the test period\. We first compute the standardised MAEe/σye/\\sigma\_\{y\}\([Gasteiger et al\., 2020](https://arxiv.org/html/2610.33209#bib.bib48)\), which expresses absolute forecasting error in units of target standard deviation\. We further transform this into the scaled MAE,s=σy/\(σy\+e\)∈\(0,1\]s=\\sigma\_\{y\}/\(\\sigma\_\{y\}\+e\)\\in\(0,1\], with higher values indicating better forecasting performance\.

Forecasting configuration\.For the controlled synthetic systems, we use ridge regression as an interpretable forecasting baseline, with the penalty coefficient set to11\. For the real\-world experiments, we use the performance of multivariable forecasting models as an empirical reference to evaluate the predictability scores\. Predictability is computed on the training set, and forecasting performance is evaluated on the test set\. Training settings and the Weather and ILI data splits are detailed in Appendix[B](https://arxiv.org/html/2610.33209#A2)\.

### 4\.3Predictability under Controlled Conditions with Synthetic Data Sets

![Refer to caption](https://arxiv.org/html/2610.33209v1/simple_synthetic_data.png)Figure 2:Basic synthetic time series used to validateMAPS\. The exogenous variables are shown in the top row, the corresponding target series in the middle row, and the predictability w\.r\.t\.HHin the bottom row, each is a combination of the exogenous inputs\. Predictability scores are computed using a rolling window where the input includes all covariates, i\.e\.both exogenous variables and the target\.We commence by conducting experiments using66basic synthetic data sets over800800time steps with known predictability and consider a forecasting setting ofL=50L=50\. The first two are deterministic, and their expected predictability score is 1\. The3−53\-5series adds i\.i\.d\.Gaussian noise to33deterministic functions, making the series less predictable\. The sixth is only Gaussian noise and is therefore not predictable\. Outcomes are depicted in Figure[2](https://arxiv.org/html/2610.33209#S4.F2)\.MAPSassigns a score of11to both deterministic series across differentHH, about0\.770\.77to the series with added noise, and close to00to pure Gaussian noise, corroborating the expected predictability\. Contrarily, persistence reflects the pattern of the series rather than its predictability\.

Figure 3:Reference predictability compared toMAPS, ForeCA, Spectral, Persistence, and the empirical performance \(using linear correlation\) of a ridge regression model for complex synthetic time series tasks across increasing forecasting horizons \(HH\) for a fixed lookback window \(L=16L=16\)\.We further evaluateMAPSon66synthetic univariate tasks that exhibit more complex and nonlinear characteristics and hence have different degrees of predictability\. This univariate setting \(C=1C=1\) serves as a special case of the forecasting task considered byMAPS\. For each task, we generate sequences of4,0004\{,\}000time steps and compare a reference \(expected\) predictability score withMAPSw\.r\.t\.different forecasting horizons \(HH\), under a fixed look\-back window ofL=16L=16\. We define the reference predictability as

P\(L,H\)=1−𝔼\[∑h=1HVar\(yt\+h∣yt−L\+1:t\)\]∑h=1HVar⁡\(yt\+h\),P\(L,H\)=1\-\\frac\{\\mathbb\{E\}\\\!\\left\[\\sum\_\{h=1\}^\{H\}\\operatorname\{Var\}\(y\_\{t\+h\}\\mid y\_\{t\-L\+1:t\}\)\\right\]\}\{\\sum\_\{h=1\}^\{H\}\\operatorname\{Var\}\(y\_\{t\+h\}\)\}\\,\\,,\(3\)whereyty\_\{t\}is the generated future target value at timett\. We use the known data\-generating process to obtain possible future sequences given the observations in the lookback window and calculate how much uncertainty remains across these futures\.P⁡\(L,H\)∈\[0,1\]P\(L,H\)\\in\[0,1\]quantifies uncertainty after observing past values within the look\-back windowLL\. A value of11indicates that the past fully determines the target value, while00dictates that the past provides no predictive information \(see also Appendix[A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1)\)\.

Figure[3](https://arxiv.org/html/2610.33209#S4.F3)compares the reference predictability withMAPS,22frequency\-based baselines \(ForeCA and Spectral\), persistence baseline, alongside the average correlation between predictions and ground truth for a ridge regression model as a baseline\. It can be observed thatMAPSstrongly correlates withP⁡\(L,H\)P\(L,H\)\(r=\.963r=\.963\)\. Their scores actually overlap, with the exception of the quasi\-periodic task where MAPS slightly overestimates predictability\. Ridge regression performance estimates exhibit slightly weaker agreement with the target predictability than MAPS, but still a significant one \(r=\.924r=\.924\)\. This alignment further motivates the use of forecasting model performance as the reference predictability in the follow\-up experiments\. Persistence achieves a lower correlation \(r=\.852r=\.852\) and oftentimes underestimates predictability, particularly at longer forecasting horizons\. Given this mediocre performance and the fact that its score reflects the error of the persistence predictor rather than data predictability, we do not report Persistence in the following experiments\. Additionally, the two spectral baselines perform poorly on datasets without strong seasonality and fail to capture how predictability varies across differentL/HL/Hsettings\.

![Refer to caption](https://arxiv.org/html/2610.33209v1/combined_weather.png)Figure 4:Predictability scores ofMAPSand55baselines for theWEAtask compared against the forecasting performance of Sonnet across2020locations \(lookbackL=56L\{=\}56, horizonH=12H\{=\}12\) for two configurations \(55or4545covariates\)\. Location order reflects Sonnet’s correlation performance\.
### 4\.4Evaluating Predictability Against Empirical Forecasting Performance

We further evaluate predictability measurements on real\-world data sets \(WEAandILI\)\. ForWEA, we first average predictability and performance scores across the test periods and then compute their correlation over2020locations for each forecasting horizon\. We perform this for two configurations, considering just the target location \(55covariates\) or including a grid of88neighbouring locations \(4545covariates\)\. Results are shown in Figure[4](https://arxiv.org/html/2610.33209#S4.F4)and Table[1](https://arxiv.org/html/2610.33209#S4.T1)\. Across the compared methods,MAPSshows the strongest correlation with the performance \(ρ\\rho\) of Sonnet\. Among the baselines, Epiplexity and ForeCA are the most competitive\. However, Epiplexity is model\-dependent, meaning that its effectiveness is closely related to the model used\.

Additionally, including climate indicators from neighbouring regions \(45 covariates\) is expected to improve the reference performanceρ\\rhoover using only local indicators \(5 covariates\) by providing additional predictive information\. Importantly,MAPSreflects this improvement\. The configuration with more covariates yields higherMAPSscores\. The baseline measures largely fail to capture this distinction, except for the model\-based method Epiplexity and MI\.

Table 1:Average linear correlation \(across all test periods\) between predictability scores and Sonnet’s forecasting performance, which is measured by correlation \(ρ\\rho\) or scaled MAE \(ss\)\. Correlation breakdowns for lookback windows \(LL\), forecasting horizons \(HH\), and covariate number \(CC\) are provided for theWEAtask, over the2020considered locations\. The right\-most two columns enumerate the average correlation among all combinations of\(L,H\)\(L,H\)across regions forWEAandILI\.Results across all forecasting configurations ofWEAare enumerated in Table[1](https://arxiv.org/html/2610.33209#S4.T1)\. When we compare the predictability score with target correlationρ\\rho,MAPSyields the strongest correlationrrin 9 out of 18 configurations\. It performs best in five of six configurations withH=12H=12and in four out of six configurations withH=28H=28\. Epiplexity is the most competitive baseline when referred toρ\\rho, obtaining stronger results for the shortest forecasting horizon\(H=4\)\(H=4\)\. Turning to the scaled\-MAE scoress,MAPSyields the strongest correlation in 13 out of 18 configurations\. ForeCA is the most competitive baseline forss, achieving the highest correlation in six configurations, while Epiplexity leads in one\. Notably,MAPSis more competitive at longer forecasting horizons\. It leads or ties in five of six settings atH=12H=12and four of six atH=28H=28for target correlationρ\\rho, and in five and six settings, respectively, for scaled MAE\.MAPS’s alignment withssis also significantly stronger than Epiplexity’s \(pairedtt\-test,p<0\.001p<0\.001\)\. The final two columns provide a complementary within\-location analysis across differentL,HL,Hsettings\. Baselines that do not produce task\-specific scores \(ForeCA and Spectral\) are marked “–” because their scores remain constant\.MAPSachieves the highest within\-location correlation inWEA/ILI, reaching0\.961/0\.8270\.961/0\.827withρ\\rho, and0\.946/0\.8170\.946/0\.817withss\. These results show that the alignment betweenMAPSand forecasting performance is consistent across datasets and\(L,H\)\(L,H\)settings, capturing changes in forecasting difficulty across tasks\.

### 4\.5Predictability Under Distribution Shift

Figure 5:Correlation between predictability score and downstream forecasting performance across weather regions under distribution shift\. \(A\) Predictability computed on the source \(training\) data\. \(B\) Predictability computed on the combined training and test data\.Table 2:Linear correlationrrbetween each predictability measure and the downstream forecasting performance across the same cross\-location weather tasks\.

Table 3:Component ablation of MAPS\. Correlations are averaged over all combinations of\(L,H\)\(L,H\)across locations forWEAandILI\. Best results arebolded\.
Table 4:MAPS predictability scores for England ILI under different PCA variance thresholds for each\(L,H\)\(L,H\)task\. “No PCA” uses all input variables, while the other rows retain the fewest principal components needed to explain70%70\\%,80%80\\%, or90%90\\%of the variance\.

We evaluate whether dataset predictability provides a useful prior for forecasting performance under real\-world distribution shift\. Specifically, we train a forecasting model on one region’s weather data and evaluate it on another, constructing a diverse set of cross\-location transfer tasks to approximate the distribution shift\. To ensure a comprehensive evaluation, we train the forecasting model on each of the 20 locations and evaluate it on five other locations, selecting the two geographically nearest and three geographically farthest locations for each case\.

Figure[5](https://arxiv.org/html/2610.33209#S4.F5)\.A shows that the predictability score of the training data alone shows a weak relationship with downstream performance \(r=0\.088r=0\.088\)\. This is expected as the source dataset provides limited information about the predictability of the target series if there is a distribution shift\. In contrast, when predictability is estimated on the combined training and deployment data \(Figure[5](https://arxiv.org/html/2610.33209#S4.F5)\.B\), a slightly more positive correlation emerges \(r=0\.684r=0\.684\)\. Table[2](https://arxiv.org/html/2610.33209#S4.T2)also shows that other baselines present a similar change\. It can be concluded that the current predictability scores measure the intrinsic predictability of the time series on which they are computed, and they do so accurately once that distribution is observed\. However, it should be noted that in deployment settings subject to shift, a score computed at training time would overstate or understate the difficulty of the task\.

### 4\.6Ablation study forMAPS

In this section, we perform an ablation study by replacing one component ofMAPSat a time to justify our design choice\. ‘w/ variance’ replaces the direction\-based similarity measure with a variance\-ratio objective using the mean variance of the target and neighbouring features\. Further detail in[B\.2](https://arxiv.org/html/2610.33209#A2.SS2)\. ‘w/ signed directions’ is instead of taking the absolute value, we map\[−1,1\]\[\-1,1\]directly onto\[0,1\]\[0,1\], so a negative alignment is mapped below0\.50\.5\. ‘w/o standardisation’ removes the z\-score standardisation applied before window slicing\. ‘w/o window standardisation’ removes the window\-wise standardisation for each sliced window\. Finally, ‘w/o PCA’ omits the PCA projection before the nearest\-neighbour search\. Table[4](https://arxiv.org/html/2610.33209#S4.T4)shows how each design choice performs under controlled predictability scenarios and contributes toMAPS\.MAPSscore consistently achieves the strongest alignment with downstream forecasting performance on the real\-world datasets, obtaining the highest correlation on bothWEA\(r=0\.961r=0\.961\) andILI\(r=0\.827r=0\.827\)\. Removing the standardisation has little effect, while removing PCA substantially degrades performance alignment on theILIdataset \(the dataset that contains more covariates\), suggesting that dimensionality reduction improves the reliability of neighbourhood estimation in higher\-dimensional input spaces\. Replacing the neighbourhood direction alignment with a variance measurement or using signed directions also reduces agreement with forecasting performance\.

### 4\.7Sensitivity analysis

Figure 6:Effect of window stride onMAPSacross nine\(L,H\)\(L,H\)settings across2020locations \(WEA\)\. \(A\) Averaged predictability score; \(B\) Correlation across2020locations\. The axis breaks at stride 100, which separates the equally spaced stride values from a linear scale over 100–500\.Effect of PCA dimensionality\.InMAPS, PCA is applied to remove redundancy among all available variables, thereby reducing the computational cost of the followingkk\-nearest\-neighbour search\. This is particularly important for theILIdataset, where each dataset contains2​k2\\text\{k\}–5​k5\\text\{k\}covariates\. Table[4](https://arxiv.org/html/2610.33209#S4.T4)reports theMAPSscore after retaining different proportions of the variance\. The score is consistently higher after PCA than the no PCA version\. We note that these results present as a sensitivity analysis ofMAPSto the input\. They do not, by themselves, show that a higher score is more accurate\.

Effect of window stride\.As shown in Figure[6](https://arxiv.org/html/2610.33209#S4.F6)\.A, increasing the stride \(MM\) lowers the absolute predictability score because adjacent windows overlap less, making neighbouring samples less similar\. Despite this change in score magnitude, its correlation with downstream forecasting performance remains stable across a wide range of strides \(Figure[6](https://arxiv.org/html/2610.33209#S4.F6)\.B\)\. Stride 1 performs slightly worse, as heavily overlapping windows may be selected as neighbours rather than revealing recurring patterns\. Correlation begins to drop at larger strides \(250 and 500\), where too few windows remain to accurately capture the local structure of the dataset\. Arguably, other factors such as dataset size and seasonality may also contribute to this decline, though we leave a closer analysis to future work\. Overall,MAPSis robust to the choice of stride\.

## 5Conclusion

We proposeMAPS, a model\-agnostic measure for quantifying the predictability of a multivariable time\-series forecasting task given its covariates, lookback window, and forecasting horizon\. Across synthetic systems with known predictability,MAPSclosely follows the reference predictability, achieving an average linear correlation ofr=0\.963r=0\.963\. Furthermore, across Weather and Influenza\-Like Illness forecasting tasks,MAPSshows the strongest relationship with empirical forecasting performance among the evaluated predictability measures, obtaining correlations ofr=0\.961r=0\.961andr=0\.827r=0\.827, respectively\. Overall,MAPSprovides a complementary measure for characterising forecasting tasks and interpreting model performance across heterogeneous datasets\.

### Acknowledgments

V\. Lampos would like to thank all levels of support from grant EP/X031276/1 \(EPSRC\) and the SOFAIR Lab \(UKRI\)\.

## References

- T\. Aksu, G\. Woo, J\. Liu, X\. Liu, C\. Liu, S\. Savarese, C\. Xiong, and D\. SahooGift\-eval: A benchmark for general time series forecasting model evaluation\.arXiv preprint arXiv:2410\.10393\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p4.1)\.
- Bandt and Pompe \(2002\)C\. Bandt and B\. PompePermutation Entropy: A Natural Complexity Measure for Time Series\.Phys\. Rev\. Lett\.88,pp\. 174102\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p2.1)\.
- Becket al\.\(2025\)N\. Beck, J\. Dovern, and S\. VoglMind the naive forecast\! a rigorous evaluation of forecasting models for time series with low predictability\.Applied Intelligence55\(6\),pp\. 395\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p2.1)\.
- Catt \(2026\)P\. M\. CattOn the Limits of Prediction: Forecastability Profiles and Information Decay in Time Series\.arXiv preprint arXiv:2603\.20546\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p3.1)\.
- Changet al\.\(2004\)P\. Chang, R\. Saravanan, T\. DelSole, and F\. WangPredictability of Linear Coupled Systems\. Part I: Theoretical Analyses\.Journal of Climate17\(7\),pp\. 1474–1486\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p1.1)\.
- DelSole \(2004\)T\. DelSolePredictability and Information Theory\. Part I: Measures of Predictability\.Journal of the Atmospheric Sciences61\(20\),pp\. 2425–2440\.Cited by:[§B\.3](https://arxiv.org/html/2610.33209#A2.SS3.p6.1)\.
- Devroyeet al\.\(1994\)L\. Devroye, L\. Gyorfi, A\. Krzyzak, and G\. LugosiOn the Strong Universal Consistency of Nearest Neighbor Regression Function Estimates\.The Annals of Statistics22\(3\),pp\. 1371 – 1385\.Cited by:[§3](https://arxiv.org/html/2610.33209#S3.p3.1)\.
- Donget al\.\(2025\)C\. Dong, D\. Faranda, A\. Gualandi, V\. Lucarini, and G\. MengaldoTime\-lagged recurrence: A data\-driven method to estimate the predictability of dynamical systems\.Proceedings of the National Academy of Sciences122\(20\),pp\. e2420252122\.Cited by:[§4\.1](https://arxiv.org/html/2610.33209#S4.SS1.p1.1)\.
- Eckmann and Ruelle \(1985\)J\. \-P\. Eckmann and D\. RuelleErgodic theory of chaos and strange attractors\.Rev\. Mod\. Phys\.57,pp\. 617–656\.Cited by:[§3\.1](https://arxiv.org/html/2610.33209#S3.SS1.p1.1)\.
- Finziet al\.\(2026\)M\. Finzi, S\. Qiu, Y\. Jiang, P\. Izmailov, J\. Z\. Kolter, and A\. G\. WilsonFrom Entropy to Epiplexity: Rethinking Information for Computationally Bounded Intelligence\.arXiv preprint arXiv:2601\.03220\.Cited by:[§B\.3](https://arxiv.org/html/2610.33209#A2.SS3.p7.1),[§4\.2](https://arxiv.org/html/2610.33209#S4.SS2.p1.1)\.
- Gaoet al\.\(2015\)S\. Gao, G\. Ver Steeg, and A\. GalstyanEfficient Estimation of Mutual Information for Strongly Dependent Variables\.InProceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics,Proceedings of Machine Learning Research, Vol\.38,pp\. 277–286\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p3.1)\.
- Gasteigeret al\.\(2020\)J\. Gasteiger, J\. Groß, and S\. GünnemannDirectional Message Passing for Molecular Graphs\.InInternational Conference on Learning Representations,Cited by:[§4\.2](https://arxiv.org/html/2610.33209#S4.SS2.p3.1)\.
- Gneitinget al\.\(2006\)T\. Gneiting, K\. Larson, K\. Westrick, M\. G\. Genton, and E\. AldrichCalibrated Probabilistic Forecasting at the Stateline Wind Energy Center\.Journal of the American Statistical Association101\(475\),pp\. 968–979\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Godahewaet al\.\(2021\)R\. Godahewa, C\. Bergmeir, G\. I\. Webb, R\. J\. Hyndman, and P\. Montero\-MansoMonash time series forecasting archive\.arXiv preprint arXiv:2105\.06643\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p4.1)\.
- Goerg \(2013\)G\. GoergForecastable Component Analysis\.InProceedings of the 30th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.28,pp\. 64–72\.Cited by:[§B\.3](https://arxiv.org/html/2610.33209#A2.SS3.p3.1),[§1](https://arxiv.org/html/2610.33209#S1.p2.1),[§4\.2](https://arxiv.org/html/2610.33209#S4.SS2.p1.1)\.
- Gordonet al\.\(1993\)N\.J\. Gordon, D\.J\. Salmond, and A\.F\.M\. SmithNovel approach to nonlinear/non\-Gaussian Bayesian state estimation\.IEE Proceedings F \(Radar and Signal Processing\)140,pp\. 107–113\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p3.1)\.
- Hendersonet al\.\(2000\)R\. Henderson, P\. Diggle, and A\. DobsonJoint modelling of longitudinal measurements and event time data\.Biostatistics1\(4\),pp\. 465–480\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Inmanet al\.\(2013\)R\. H\. Inman, H\. T\.C\. Pedro, and C\. F\.M\. CoimbraSolar forecasting methods for renewable energy integration\.Progress in Energy and Combustion Science39\(6\),pp\. 535–576\.External Links:ISSN 0360\-1285Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Jolliffe \(1986\)I\. T\. JolliffePrincipal Components in Regression Analysis\.InPrincipal Component Analysis,pp\. 129–155\.Cited by:[§3](https://arxiv.org/html/2610.33209#S3.p2.1)\.
- Kitagawa \(1996\)G\. KitagawaMonte Carlo Filter and Smoother for Non\-Gaussian Nonlinear State Space Models\.Journal of Computational and Graphical Statistics5\(1\),pp\. 1–25\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p3.1)\.
- Kraskovet al\.\(2004\)A\. Kraskov, H\. Stögbauer, and P\. GrassbergerEstimating mutual information\.Phys\. Rev\. E69,pp\. 066138\.Cited by:[§B\.3](https://arxiv.org/html/2610.33209#A2.SS3.p6.1),[§1](https://arxiv.org/html/2610.33209#S1.p3.1)\.
- Lehmanet al\.\(2014\)L\. H\. Lehman, R\. P\. Adams, L\. Mayaud, G\. B\. Moody, A\. Malhotra, R\. G\. Mark, and S\. NematiA physiological time series dynamics\-based approach to patient monitoring and outcome prediction\.IEEE J Biomed Health Inform19\(3\),pp\. 1068–1076\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Levadaet al\.\(2024\)A\. L\. M\. Levada, F\. Nielsen, and M\. F\. C\. HaddadAdaptivekk\-nearest neighbor classifier based on the local estimation of the shape operator\.arXiv preprint arXiv:2409\.05084\.Cited by:[§3](https://arxiv.org/html/2610.33209#S3.p3.1)\.
- Lorenz \(1963\)E\. N\. LorenzDeterministic Nonperiodic Flow\.Journal of Atmospheric Sciences20\(2\),pp\. 130–141\.Cited by:[§3\.1](https://arxiv.org/html/2610.33209#S3.SS1.p1.1)\.
- Parket al\.\(2011\)C\. Park, J\. Z\. Huang, and Y\. DingDomain Decomposition Approach for Fast Gaussian Process Regression of Large Spatial Data Sets\.Journal of Machine Learning Research12\(47\),pp\. 1697–1728\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p1.2)\.
- Pincus \(1991\)S\. M\. PincusApproximate entropy as a measure of system complexity\.\.Proceedings of the National Academy of Sciences88\(6\),pp\. 2297–2301\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p2.1)\.
- Powell and Percival \(1979\)G\. E\. Powell and I\. C\. PercivalA spectral entropy method for distinguishing regular and irregular motion of Hamiltonian systems\.Journal of Physics A: Mathematical and General12\(11\),pp\. 2053\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p2.1)\.
- Rabiner \(1989\)L\.R\. RabinerA tutorial on hidden Markov models and selected applications in speech recognition\.Proceedings of the IEEE77\(2\),pp\. 257–286\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p2.1)\.
- Ranjanet al\.\(2011\)P\. Ranjan, R\. Haynes, and R\. KarstenA Computationally Stable Approach to Gaussian Process Interpolation of Deterministic Computer Simulation Data\.Technometrics53\(4\),pp\. 366–378\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p1.4)\.
- Raspet al\.\(2020\)S\. Rasp, P\. D\. Dueben, S\. Scher, J\. A\. Weyn, S\. Mouatadid, and N\. ThuereyWeatherBench: a benchmark data set for data\-driven weather forecasting\.Journal of Advances in Modeling Earth Systems12\(11\),pp\. e2020MS002203\.Cited by:[§A\.4](https://arxiv.org/html/2610.33209#A1.SS4.p1.1)\.
- Shchuret al\.\(2025\)O\. Shchur, A\. F\. Ansari, C\. Turkmen, L\. Stella, N\. Erickson, P\. Guerron, M\. Bohlke\-Schneider, and Y\. Wangfev\-bench: A realistic benchmark for time series forecasting\.arXiv preprint arXiv:2509\.26468\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p4.1)\.
- Shu and Lampos \(2025\)Y\. Shu and V\. LamposDeformTime: capturing variable dependencies with deformable attention for time series forecasting\.Transactions on Machine Learning Research\.Cited by:[§A\.3](https://arxiv.org/html/2610.33209#A1.SS3.p1.1),[§B\.1](https://arxiv.org/html/2610.33209#A2.SS1.p1.1)\.
- Shu and Lampos \(2026\)Y\. Shu and V\. LamposSonnet: Spectral Operator Neural Network for Multivariable Time Series Forecasting\.Proceedings of the AAAI Conference on Artificial Intelligence\.Cited by:[§A\.1](https://arxiv.org/html/2610.33209#A1.SS1.p1.1),[§A\.4](https://arxiv.org/html/2610.33209#A1.SS4.p1.1),[§B\.1](https://arxiv.org/html/2610.33209#A2.SS1.p1.1),[§4\.1](https://arxiv.org/html/2610.33209#S4.SS1.p2.1)\.
- Sims \(1980\)C\. A\. SimsMacroeconomics and Reality\.Econometrica48\(1\),pp\. 1–48\.External Links:ISSN 00129682, 14680262Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Stock and Watson \(2002a\)J\. H\. Stock and M\. W\. WatsonMacroeconomic Forecasting Using Diffusion Indexes\.Journal of Business & Economic Statistics20\(2\),pp\. 147–162\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Stock and Watson \(2002b\)J\. H\. Stock and M\. W\. WatsonForecasting Using Principal Components from a Large Number of Predictors\.Journal of the American Statistical Association97\(460\),pp\. 1167–1179\.External Links:ISSN 01621459Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Stone \(1977\)C\. J\. StoneConsistent Nonparametric Regression\.The Annals of Statistics5\(4\),pp\. 595–620\.Cited by:[§3\.1](https://arxiv.org/html/2610.33209#S3.SS1.p1.1),[§3](https://arxiv.org/html/2610.33209#S3.p3.1)\.
- Sugiharaet al\.\(1994\)G\. Sugihara, B\. T\. Grenfell, R\. M\. May, and H\. TongNonlinear forecasting for the classification of natural time series\.Philosophical Transactions of the Royal Society of London, Series A: Physical and Engineering Sciences348\(1688\),pp\. 477–495\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p2.1)\.
- Takens \(1981\)F\. TakensDetecting strange attractors in turbulence\.InDynamical Systems and Turbulence, Warwick 1980,pp\. 366–381\.External Links:ISBN 978\-3\-540\-38945\-3Cited by:[§3\.1](https://arxiv.org/html/2610.33209#S3.SS1.p1.1)\.
- Tereninet al\.\(2024\)A\. Terenin, D\. R\. Burt, A\. Artemev, S\. Flaxman, M\. van der Wilk, C\. E\. Rasmussen, and H\. GeNumerically Stable Sparse Gaussian Processes via Minimum Separation using Cover Trees\.Journal of Machine Learning Research25\(26\),pp\. 1–36\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p1.4)\.
- Wanget al\.\(2025a\)O\. Wang, P\. Quan, K\. Yang, and M\. SrivastavaSpectral Predictability as a Fast Reliability Indicator for Time Series Forecasting Model Selection\.arXiv preprint arXiv:2511\.08884\.Cited by:[§A\.1](https://arxiv.org/html/2610.33209#A1.SS1.p1.1),[§B\.3](https://arxiv.org/html/2610.33209#A2.SS3.p4.1),[§1](https://arxiv.org/html/2610.33209#S1.p2.1),[§4\.2](https://arxiv.org/html/2610.33209#S4.SS2.p1.1)\.
- Wanget al\.\(2025b\)Y\. Wang, H\. Wu, Y\. Ma, Y\. Fang, Z\. Zhang, Y\. Liu, S\. Wang, Z\. Ye, Y\. Xiang, J\. Wang,et al\.Accuracy law for the future of deep time series forecasting\.arXiv preprint arXiv:2510\.02729\.Cited by:[§B\.3](https://arxiv.org/html/2610.33209#A2.SS3.p5.1),[§1](https://arxiv.org/html/2610.33209#S1.p3.1),[§4\.2](https://arxiv.org/html/2610.33209#S4.SS2.p1.1)\.
- Wenet al\.\(2019\)G\. Wen, N\. I\. Fox, and P\. S\. MarketA Gaussian mixture method for specific differential phase retrieval at X\-band frequency\.Atmospheric Measurement Techniques12,pp\. 5613–5637\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p2.1)\.
- Williams and Rasmussen \(1995\)C\. Williams and C\. RasmussenGaussian Processes for Regression\.InAdvances in Neural Information Processing Systems,Vol\.8,pp\.\.Cited by:[§A\.2\.1](https://arxiv.org/html/2610.33209#A1.SS2.SSS1.p1.2)\.
- Wulfsohn and Tsiatis \(1997\)M\. S\. Wulfsohn and A\. A\. TsiatisA Joint Model for Survival and Longitudinal Data Measured with Error\.Biometrics53\(1\),pp\. 330–339\.External Links:ISSN 0006341X, 15410420Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p1.1)\.
- Zadehet al\.\(2025\)S\. G\. Zadeh, V\. Shaj, P\. Jahnke, G\. Neumann, and T\. BreitenbachTowards Measuring Predictability: To which extent data\-driven approaches can extract deterministic relations from data exemplified with time series prediction and classification\.Transactions on Machine Learning Research\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p2.1)\.
- Zenget al\.\(2023\)A\. Zeng, M\. Chen, L\. Zhang, and Q\. XuAre Transformers Effective for Time Series Forecasting?\.InProceedings of the AAAI Conference on Artificial Intelligence,pp\. 11121–11128\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p4.1)\.
- Zhanget al\.\(2024\)J\. Zhang, X\. Wen, Z\. Zhang, S\. Zheng, J\. Li, and J\. BianProbTS: Benchmarking Point and Distributional Forecasting across Diverse Prediction Horizons\.InAdvances in Neural Information Processing Systems,Vol\.37,pp\. 48045–48082\.Cited by:[§1](https://arxiv.org/html/2610.33209#S1.p4.1)\.
- Zhang and Yan \(2023\)Y\. Zhang and J\. YanCrossformer: Transformer Utilizing Cross\-Dimension Dependency for Multivariate Time Series Forecasting\.InInternational Conference on Learning Representations,Cited by:[§C\.4](https://arxiv.org/html/2610.33209#A3.SS4.p1.1)\.

## Appendix

## Appendix ADatasets and Metrics

### A\.1Choice of model performance metrics

Figure S1:Relationship between SMAPE and linear correlation across weather forecasting tasks for Sonnet against target correlationρ\\rho\(left\) and scaled MAE scoress\(right\)\. Each point represents one task for one of 20 locations, with colours indicating different locations\.The proposedMAPSscore is designed to quantify the intrinsic predictability of a time series under a rolling\-window setting in a model\-agnostic way\. As described in Section[4\.2](https://arxiv.org/html/2610.33209#S4.SS2), for real\-world datasets, where reference predictability is difficult to quantify directly, we assess the score by comparing it with the performance of a forecasting model on the same task under the metrics of linear correlation \(ρ\\rho\) and scaled MAE \(ss\)\. Notably, the SMAPE score of ML models has recently been used to benchmark predictability scores\[[Wang et al\., 2025a](https://arxiv.org/html/2610.33209#bib.bib13)\]\. In Figure[S1](https://arxiv.org/html/2610.33209#A1.F1), we compare SMAPE with linear correlation \(ρ\\rho\) and scaled MAE \(ss\) for Sonnet\[[Shu and Lampos, 2026](https://arxiv.org/html/2610.33209#bib.bib14)\], where each location is shown in a different colour\. Within each location, bothρ\\rhoandssshow a negative relationship with SMAPE as expected\. Models with higher correlation tend to have lower prediction errors under the same data distribution\. However, this relationship becomes weaker when results are evaluated across locations\. This suggests that these metrics capture different aspects of forecasting performance: a predictability score that correlates well withρ\\rhoandssneed not correlate well with SMAPE\. For completeness, we additionally report the relationship between the proposed predictability score and target SMAPE in Appendix[C\.3](https://arxiv.org/html/2610.33209#A3.SS3)\.

### A\.2Generation of Complex Synthetic Datasets with Known Predictability

Figure S2:Generated time series for the six synthetic processes with known predictability\. Each series containsT=4,000T=4\{,\}000time steps, of which the first600600are shown\.We construct six synthetic time series with known predictability to evaluate the proposed method under controlled conditions\. Figure[S2](https://arxiv.org/html/2610.33209#A1.F2)shows the generated series\. For the linear Gaussian state\-space model \(SSM\), smooth Gaussian Process \(GP\), fractional Gaussian noise, and quasi\-periodic process, we characterise temporal dependence through autocovariance functions\. Throughout this section,yt∈ℝy\_\{t\}\\in\\mathbb\{R\}denotes the observed univariate series at timett, andγ⁡\(h\)=Cov⁡\(yt,yt\+h\)∈ℝ\\gamma\(h\)=\\operatorname\{Cov\}\(y\_\{t\},y\_\{t\+h\}\)\\in\\mathbb\{R\}denotes its autocovariance function at laghh, which does not depend onttbecause the process is stationary\.

Linear Gaussian state\-space model\.We first generate synthetic time series using a22\-dimensional linear Gaussian state\-space model\. Letzt=\(zt,1,zt,2\)⊤∈ℝ2z\_\{t\}=\(z\_\{t,1\},z\_\{t,2\}\)^\{\\top\}\\in\\mathbb\{R\}^\{2\}be the latent state at timett\. The state transition is given by

zt\+1=A​zt\+wt,yt=zt,1\+vt,z\_\{t\+1\}=Az\_\{t\}\+w\_\{t\},\\qquad y\_\{t\}=z\_\{t,1\}\+v\_\{t\}\\,\\,,\(S1\)wherewt∼𝒩⁡\(0,σw2​I2\)w\_\{t\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{w\}^\{2\}I\_\{2\}\)andvt∼𝒩⁡\(0,σv2\)v\_\{t\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{v\}^\{2\}\)are independent Gaussian noise\. The transition matrix is defined as

A=κ​\(cos⁡\(ω\)−sin⁡\(ω\)sin⁡\(ω\)cos⁡\(ω\)\),ω=2​π​f,A=\\kappa\\begin\{pmatrix\}\\cos\(\\omega\)&\-\\sin\(\\omega\)\\\\ \\sin\(\\omega\)&\\phantom\{\-\}\\cos\(\\omega\)\\end\{pmatrix\},\\qquad\\omega=2\\pi f\\,\\,,\(S2\)whereffcontrols the frequency andκ<1\\kappa<1controls the rate at which the process loses temporal dependence\. This construction produces a stationary, stochastically perturbed oscillatory process with exponentially decaying temporal dependence, whileσv2\\sigma\_\{v\}^\{2\}controls the noise level\. In our experiments, we setf=0\.06f=0\.06,κ=0\.985\\kappa=0\.985,σw2=1\\sigma\_\{w\}^\{2\}=1, andσv2=0\.35\\sigma\_\{v\}^\{2\}=0\.35, and initialise withz0=\(0,0\)⊤z\_\{0\}=\(0,0\)^\{\\top\}\.

Smooth Gaussian Process\.The smooth GP is defined by the squared\-exponential covariance function, given by:

γ⁡\(h\)=\{1\+σϵ2h=0exp⁡\(−h22​ℓ2\)h≠0,\\gamma\(h\)=\\begin\{cases\}1\+\\sigma\_\{\\epsilon\}^\{2\}&h=0\\\\ \\displaystyle\\exp\\left\(\-\\frac\{h^\{2\}\}\{2\\ell^\{2\}\}\\right\)&h\\neq 0\\,\\,,\\end\{cases\}\(S3\)wherehhdenotes the temporal lag between two time points\.ℓ\\ellcontrols the rate at which temporal correlations decay with increasing time lag, whileσϵ2\\sigma\_\{\\epsilon\}^\{2\}controls the variance of the noise, which comes into play only whenh=0h=0\. In our experiments, we setℓ=12\\ell=12andσϵ2=0\.05\\sigma\_\{\\epsilon\}^\{2\}=0\.05\.

Fractional Gaussian noise process\.The next synthetic time series is generated using fractional Gaussian noise, a zero\-mean stationary GP with covariance function

γ⁡\(h\)=\{1\+σϵ2h=012​\[\(h\+1\)2​η−2​h2​η\+\|h−1\|2​η\]h≠0,\\gamma\(h\)=\\begin\{cases\}1\+\\sigma\_\{\\epsilon\}^\{2\}&h=0\\\\ \\displaystyle\\frac\{1\}\{2\}\\left\[\(h\+1\)^\{2\\eta\}\-2h^\{2\\eta\}\+\|h\-1\|^\{2\\eta\}\\right\]&h\\neq 0\\,\\,,\\end\{cases\}\(S4\)wherehhdenotes the temporal lag andη∈\(0,1\)\\eta\\in\(0,1\)controls how much a time series depends on the past values\. Forη\>1/2\\eta\>1/2, the process exhibits long\-range dependence, with temporal correlations reducing according to a power law\. Larger values ofη\\etatherefore produce stronger and more persistent dependence across long time lags\. The termσϵ2\\sigma\_\{\\epsilon\}^\{2\}is the noise variance\. We setη=0\.95\\eta=0\.95andσϵ2=0\.01\\sigma\_\{\\epsilon\}^\{2\}=0\.01\.

Two\-state hidden Markov model\.The next series is generated by a two\-state hidden Markov model with Gaussian emissions\. At each time pointtt, the process is assumed to be in one of two possible states, denoted byzt∈\{1,2\}z\_\{t\}\\in\\\{1,2\\\}\. The transition probabilities are given by

Pr⁡\(zt=j∣zt−1=i\)=Ai​j,A=\(λ1−λ1−λλ\),\\Pr\(z\_\{t\}=j\\mid z\_\{t\-1\}=i\)=A\_\{ij\},\\qquad A=\\begin\{pmatrix\}\\lambda&1\-\\lambda\\\\ 1\-\\lambda&\\lambda\\end\{pmatrix\}\\,\\,,\(S5\)wherePr⁡\(⋅\)\\Pr\(\\cdot\)denotes probability,Ai​jA\_\{ij\}is the probability of transitioning from stateiiat timet−1t\-1to statejjat timett, andλ\\lambdais the probability of remaining in the same state\. Given the current state, the time series is generated from a state\-specific Gaussian distribution:

yt\|zt=i∼𝒩⁡\(μi,σi2\),y\_\{t\}\\mid z\_\{t\}=i\\sim\\mathcal\{N\}\(\\mu\_\{i\},\\sigma\_\{i\}^\{2\}\)\\,\\,,\(S6\)whereμi\\mu\_\{i\}is the mean of stateii, andσi\\sigma\_\{i\}is its standard deviation\. The two hidden states have their own mean and standard deviation, allowing the series to switch between two regimes\. We setλ=0\.99\\lambda=0\.99,\(μ1,μ2\)=\(−2,2\)\(\\mu\_\{1\},\\mu\_\{2\}\)=\(\-2,2\), and\(σ1,σ2\)=\(1\.4,1\.4\)\(\\sigma\_\{1\},\\sigma\_\{2\}\)=\(1\.4,1\.4\)\.

Noisy Lorenz system\.This series is generated through a stochastic Lorenz state\-space model\. Letzt∈ℝ3z\_\{t\}\\in\\mathbb\{R\}^\{3\}denote the latent state, which evolves according to

zt\+1=ΦRK4​\(zt\)\+wt,yt=zt,1\+vt,z\_\{t\+1\}=\\Phi\_\{\\mathrm\{RK4\}\}\(z\_\{t\}\)\+w\_\{t\},\\qquad y\_\{t\}=z\_\{t,1\}\+v\_\{t\}\\,\\,,\(S7\)wherezt\+1z\_\{t\+1\}is the state at the next time point,ΦRK4\\Phi\_\{\\mathrm\{RK4\}\}is the function that takes the current stateztz\_\{t\}and computes its next state according to the Lorenz differential equations using the fourth\-order Runge–Kutta \(RK4\) method,wt∼𝒩⁡\(0,σw2​I3\)w\_\{t\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{w\}^\{2\}I\_\{3\}\)is process noise that introduces randomness into the state evolution,zt,1z\_\{t,1\}is the first component of the latent state, andvt∼𝒩⁡\(0,σv2\)v\_\{t\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{v\}^\{2\}\)is the added noise\. We use the standard Lorenz parameters\(σ,ρ,β\)=\(10,28,8/3\)\(\\sigma,\\rho,\\beta\)=\(10,28,8/3\), an integration step ofΔ​t=0\.05\\Delta t=0\.05, and noise varianceσw2=0\.02\\sigma\_\{w\}^\{2\}=0\.02andσv2=4\\sigma\_\{v\}^\{2\}=4\. We initialise withz0=\(1,1,1\)⊤z\_\{0\}=\(1,1,1\)^\{\\top\}\.

Quasi\-periodic process\.We next combine several sinusoidal components with frequencies that do not have a common period\. This yields

yt=∑j=1Jaj​cos⁡\(2​π​fj​t\+ϕj\)\+ϵt,y\_\{t\}=\\sum\_\{j=1\}^\{J\}a\_\{j\}\\cos\(2\\pi f\_\{j\}t\+\\phi\_\{j\}\)\+\\epsilon\_\{t\}\\,\\,,\(S8\)whereJJis the number of sinusoidal components,aja\_\{j\}is the amplitude of componentjj,fjf\_\{j\}is its frequency,ϕj\\phi\_\{j\}is its fixed phase, andϵt∼𝒩⁡\(0,σϵ2\)\\epsilon\_\{t\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{\\epsilon\}^\{2\}\)is independent noise\. Because the frequencies do not share a common period, the combined signal does not exactly repeat over time\. The underlying signal is fully deterministic because the amplitudes, frequencies, and phases are fixed, while the added noise can make the observed series appear irregular\. We setJ=2J=2,\(a1,a2\)=\(1,0\.7\)\(a\_\{1\},a\_\{2\}\)=\(1,0\.7\),f1=0\.017f\_\{1\}=0\.017,f2=\(5−1\)​f1f\_\{2\}=\(\\sqrt\{5\}\-1\)f\_\{1\}, andσϵ2=0\.01\\sigma\_\{\\epsilon\}^\{2\}=0\.01\.

#### A\.2\.1Reference Predictability

Following[Chang et al\. \[2004\]](https://arxiv.org/html/2610.33209#bib.bib35), we define reference predictability as one minus the variance that remains across possible future sequences given the lookback window, normalised by the total variance in the future target values\. For a lookback window ofLLtime steps and a forecasting horizonHH, the reference predictability is given by

P\(L,H\)=1−𝔼\[∑h=1HVar\(yt\+h∣yt−L\+1:t\)\]∑h=1HVar⁡\(yt\+h\)\.P\(L,H\)=1\-\\frac\{\\mathbb\{E\}\\\!\\left\[\\sum\_\{h=1\}^\{H\}\\operatorname\{Var\}\(y\_\{t\+h\}\\mid y\_\{t\-L\+1:t\}\)\\right\]\}\{\\sum\_\{h=1\}^\{H\}\\operatorname\{Var\}\(y\_\{t\+h\}\)\}\\,\\,\.\(S9\)For each past input window, we compute how much the possible future target values can vary under the known data\-generating process\. The variance is computed directly for the44GP\-based and the hidden Markov models, and approximated with generated futures for the noisy Lorenz system\. The reference predictability score is then computed by11minus the normalised variance averaged across the sliced windows\.P⁡\(L,H\)=1P\(L,H\)=1thus indicates a completely predictable target, whereasP⁡\(L,H\)=0P\(L,H\)=0means that observing the past provides no information about the target value\. For each synthetic system, we first compute its conditional predictive variance\. For the linear Gaussian SSM, smooth GP, fractional Gaussian noise, and quasi\-periodic Gaussian process, the conditional variance is available analytically from the autocovariance functionγ⁡\(h\)\\gamma\(h\)\[[Williams and Rasmussen, 1995](https://arxiv.org/html/2610.33209#bib.bib37),[Park et al\., 2011](https://arxiv.org/html/2610.33209#bib.bib36)\]\. Let

Σp​q=γ⁡\(\|p−q\|\),cp=γ⁡\(H\+L−1−p\),\\Sigma\_\{pq\}=\\gamma\(\|p\-q\|\),\\qquad c\_\{p\}=\\gamma\(H\+L\-1\-p\)\\,\\,,\(S10\)wherep,q=0,…,L−1p,q=0,\\ldots,L\-1index positions in the lookback window\.Σ\\Sigmais the covariance matrix of theLLpast input\.𝐜=\(c0,…,cL−1\)⊤\\mathbf\{c\}=\(c\_\{0\},\\ldots,c\_\{L\-1\}\)^\{\\top\}is a vector that describes how strongly each of theLLpast input is related to the target target valueyt\+Hy\_\{t\+H\}\. Because the past input and the target are jointly Gaussian, the variance remaining after observing the past is given by:

Var\(yt\+H∣yt−L\+1:t\)=γ\(0\)−𝐜⊤Σ−1𝐜\.\\operatorname\{Var\}\(y\_\{t\+H\}\\mid y\_\{t\-L\+1:t\}\)=\\gamma\(0\)\-\\mathbf\{c\}^\{\\top\}\\Sigma^\{\-1\}\\mathbf\{c\}\.\(S11\)Following[Ranjan et al\. \[2011\]](https://arxiv.org/html/2610.33209#bib.bib39)and[Terenin et al\. \[2024\]](https://arxiv.org/html/2610.33209#bib.bib38), we add10−10​γ​\(0\)10^\{\-10\}\\gamma\(0\)to diagonal ofΣ\\Sigmafor stability\.

For each lookback window ofLL, the forward algorithm\[[Rabiner, 1989](https://arxiv.org/html/2610.33209#bib.bib40)\]estimates the posterior probabilities of the two hidden states at the end of the window\. The posterior is propagatedHHsteps through the transition matrix, giving future\-state probabilitiesπt\+H,i\\pi\_\{t\+H,i\}in stateiiat timet\+Ht\+H\. Using the Gaussian\-mixture variance formula of[Wen et al\. \[2019\]](https://arxiv.org/html/2610.33209#bib.bib41), the predictive variance is given by

Var\(yt\+H∣yt−L\+1:t\)=∑i=12πt\+H,i\(μi2\+σi2\)−\(∑i=12πt\+H,iμi\)2\.\\operatorname\{Var\}\(y\_\{t\+H\}\\mid y\_\{t\-L\+1:t\}\)=\\sum^\{2\}\_\{i=1\}\\pi\_\{t\+H,i\}\(\\mu\_\{i\}^\{2\}\+\\sigma\_\{i\}^\{2\}\)\-\\left\(\\sum^\{2\}\_\{i=1\}\\pi\_\{t\+H,i\}\\mu\_\{i\}\\right\)^\{2\}\.\(S12\)We average the predictive variance over400400lookback windows and use the result to calculate the reference predictability scoreP⁡\(L,H\)P\(L,H\)\.

For the noisy Lorenz system, we generate500500possible hidden Lorenz states for each lookback window from its target values using a bootstrap particle filter\[[Gordon et al\., 1993](https://arxiv.org/html/2610.33209#bib.bib42),[Kitagawa, 1996](https://arxiv.org/html/2610.33209#bib.bib43)\]\. We then evolve these states under the known Lorenz dynamics and measure how much their future target values differ\. Since the target is the first state variable plus Gaussian noise, the predictive variance is estimated as

Var\(yt\+H∣yt−L\+1:t\)≈Var\(zt\+H,1∣yt−L\+1:t\)\+σv2\.\\operatorname\{Var\}\(y\_\{t\+H\}\\mid y\_\{t\-L\+1:t\}\)\\approx\\operatorname\{Var\}\(z\_\{t\+H,1\}\\mid y\_\{t\-L\+1:t\}\)\+\\sigma\_\{v\}^\{2\}\.\(S13\)Here,Var\\operatorname\{Var\}is the variance estimated from the propagated particles\. We average this variance over400400lookback windows for obtainingP⁡\(L,H\)P\(L,H\)\.

### A\.3Construction of theILIdatasets

Following[Shu and Lampos \[2025\]](https://arxiv.org/html/2610.33209#bib.bib5), we construct daily influenza\-like illness \(ILI\) forecasting datasets for England and US Health and Human Services Regions 2, 9, and additionally included region 10\. The forecasting target is the ILI rate for each region\. ILI rates are obtained from the Royal College of General Practitioners \(RCGP\) for England and from the Centers for Disease Control and Prevention \(CDC\) for the US regions\. Since the original ILI rates are reported weekly, we linearly interpolate them to a daily resolution\. We use daily frequencies of health\-related Google search trends as exogenous inputs\. For the US regions, search frequencies are collected at the state level and combined to obtain a regional time series\. Additionally, we consider the lag in reported ILI rates, specifically77days for England and1414days for the US regions, following the setup in DeformTime\[[Shu and Lampos, 2025](https://arxiv.org/html/2610.33209#bib.bib5)\]\.

### A\.4Selection of prediction regions from the WeatherBench dataset

Figure S3:Distribution of the 20 locations used in the weather forecasting tasks\.We form the weather forecasting \(WEA\) dataset using climate information extracted from the WeatherBench repository\[[Rasp et al\., 2020](https://arxiv.org/html/2610.33209#bib.bib32)\], covering the period from 1979 to 2018, following the data collection process of Sonnet\[[Shu and Lampos, 2026](https://arxiv.org/html/2610.33209#bib.bib14)\]\. Data is sampled from a uniform spatial grid over the globe with a resolution of5\.625∘5\.625^\{\\circ\}in both latitude and longitude, where each grid point represents a geographical location\. To capture a broader range of climatic conditions, we select2020spatially distant locations distributed across different continents, hemispheres, climate zones, and elevations, as labelled in Figure[S3](https://arxiv.org/html/2610.33209#A1.F3)\. This allows the dataset to cover diverse weather patterns\. For each location, climate data is sampled from the grid point closest to its geographical coordinates\. The forecasting target is the 850 hPa temperature in Kelvin \(T850\), a widely used indicator of atmospheric conditions\. Following previous work\[[Shu and Lampos, 2026](https://arxiv.org/html/2610.33209#bib.bib14)\], we include44additional climate variables as exogenous inputs, namely the 500 hPa geopotential \(Z500\), 2\-metre air temperature \(t2m\), 10\-metre zonal wind \(u10\), and 10\-metre meridional wind \(v10\)\. The data is resampled to a temporal resolution of 6 hours\. To investigate the contribution of spatial context, we construct two variants for each location\. The first contains 45 covariates\. It includes all55above\-mentioned variables from the88neighbouring grid points surrounding the target location\. This forms a3×33\\times 3spatial window centred on the selected grid point, resulting in 45 input variables\. The second contains55local climate variables at the target grid point, without accounting for the neighbouring grid points\.

## Appendix BSupplementary Experiment Setup

### B\.1Forecasting setup and data splits

Models are trained using Adam with a learning rate of2×10−42\\times 10^\{\-4\}\. We set the maximum training epochs to100100with early stopping based on validation loss, with a patience of55epochs\. Following Sonnet\[[Shu and Lampos, 2026](https://arxiv.org/html/2610.33209#bib.bib14)\]forWEA, years 2016, 2017, and 2018 are used for testing, with the preceding year used for validation and all earlier data \(since 1980\) used for training\. Following DeformTime\[[Shu and Lampos, 2025](https://arxiv.org/html/2610.33209#bib.bib5)\]forILI, we use44test seasons from 2015/16 to 2018/19\. The preceding99seasons are used for training, and three6060\-day periods from the most recent training seasons are selected as a held\-out set for validation\. All inputs are standardised using training\-set statistics\. For each task, the predictability score is computed on the training set, while the forecasting performance is evaluated on the test set\. All experiments are conducted using a server running Ubuntu 24\.04\.4, with 2 AMD EPYC 9354 CPUs, 768GB of DDR5 RAM, and 3 NVIDIA L40S GPUs with 48GB VRAM each\.

### B\.2Details of the ablation components with variance

We describe the variance\-based component evaluated in the ablation study in Section[4\.6](https://arxiv.org/html/2610.33209#S4.SS6), which replaces the direction\-based similarity measure used byMAPS\. Specifically, given the neighbourhood𝒩t\\mathcal\{N\}\_\{t\}defined in Section[3\.1](https://arxiv.org/html/2610.33209#S3.SS1), we compute:

vt=1H​∑h=1H1k−1​∑j∈𝒩t\(yj,h−y¯t,h\)2,y¯t,h=1k​∑j∈𝒩tyj,h\.v\_\{t\}=\\frac\{1\}\{H\}\\sum\_\{h=1\}^\{H\}\\frac\{1\}\{k\-1\}\\sum\_\{j\\in\\mathcal\{N\}\_\{t\}\}\\left\(y\_\{j,h\}\-\\bar\{y\}\_\{t,h\}\\right\)^\{2\},\\qquad\\bar\{y\}\_\{t,h\}=\\frac\{1\}\{k\}\\sum\_\{j\\in\\mathcal\{N\}\_\{t\}\}y\_\{j,h\}\.\(S14\)Here,yj,hy\_\{j,h\}denotes thehh\-th element of the standardised target vector𝐲j\\mathbf\{y\}\_\{j\}, andy¯t,h\\bar\{y\}\_\{t,h\}is its mean across the neighbours of windowtt\. Thus,vtv\_\{t\}measures the variance of the neighbouring targets averaged across theHHforecast steps\. We define the window\-level predictability as

ut=1−min⁡\(vt,1\),u\_\{t\}=1\-\\min\(v\_\{t\},1\)\\,\\,,\(S15\)and the overall variance\-based score as

𝒫var=1N​∑t=1Nut\.\\mathcal\{P\}\_\{\\mathrm\{var\}\}=\\frac\{1\}\{N\}\\sum\_\{t=1\}^\{N\}u\_\{t\}\.\(S16\)A smallervtv\_\{t\}indicates that neighbouring input windows have more clustered future targets and therefore yielding a higher predictability score\.

### B\.3Summary of baseline predictability methods

We compare against66baseline methods for dataset predictability\.

Persistence \(Naïve\)\.The persistence baseline predicts target values with the last available observation\. We thereby define persistence predictability as the reduction in absolute error relative to predicting the mean target value:

P=1−∑t=1N∑h=1H\|yt\+h−yt\|∑t=1N∑h=1H\|yt\+h−y¯\|,y¯=1N​H​∑t=1N∑h=1Hyt\+h\.P=1\-\\frac\{\\sum\_\{t=1\}^\{N\}\\sum\_\{h=1\}^\{H\}\|y\_\{t\+h\}\-y\_\{t\}\|\}\{\\sum\_\{t=1\}^\{N\}\\sum\_\{h=1\}^\{H\}\|y\_\{t\+h\}\-\\bar\{y\}\|\},\\qquad\\bar\{y\}=\\frac\{1\}\{NH\}\\sum\_\{t=1\}^\{N\}\\sum\_\{h=1\}^\{H\}y\_\{t\+h\}\.
Forecastable Component Analysis \(ForeCA\)\.Proposed by[Goerg \[2013\]](https://arxiv.org/html/2610.33209#bib.bib21), ForeCA defines the predictability of a univariate time series as the complement of the normalised Shannon entropy of its power spectral density\. We apply this measure to the complete target series for each training split and report the average across splits\.222Our implementation uses the[ForeCA package](https://cran.r-project.org/web/packages/ForeCA/index.html)\.

Spectral predictability\.Proposed by[Wang et al\. \[2025a\]](https://arxiv.org/html/2610.33209#bib.bib13), this method quantifies predictability as one minus the normalised spectral entropy of the power spectrum of a univariate time series\. For each training split, we compute the score over the complete series of each covariate and the target separately and average the resulting scores across variables\. We then report the average across training splits\.

Pattern Complexity \(PC\)\.Proposed by[Wang et al\. \[2025b\]](https://arxiv.org/html/2610.33209#bib.bib2), PC measures the variation in frequency\-domain patterns across forecasting windows, with larger values indicating greater complexity\. For each input variable, including the target history, we concatenate itsLL\-step lookback window with the correspondingHH\-step future target, compute PC across the resulting windows, and average the scores across variables\. We then report the average across training splits\.

Mutual Information \(MI\)\.Motivated by the predictability measurement proposed in[DelSole \[2004\]](https://arxiv.org/html/2610.33209#bib.bib44), we adapt mutual information to multivariable forecasting tasks by constructing paired samples from rolling windows, with anLL\-step look\-back window as input and anHH\-step target horizon as output\. We estimate MI for each input feature using[Kraskov et al\. \[2004\]](https://arxiv.org/html/2610.33209#bib.bib22)and report the average score\.

Epiplexity\.Following[Finzi et al\. \[2026\]](https://arxiv.org/html/2610.33209#bib.bib12), epiplexity measures the learnable information in a dataset, estimated as the area between the training loss curve and the final loss\. We approximate epiplexity using ridge regression as a baseline measure of dataset predictability, with higher values indicating more learnable structure\. Notably, since epiplexity is model\-dependent, it may yield different values under different architectures\.

## Appendix CSupplementary Results

### C\.1Results with synthetic data

![Refer to caption](https://arxiv.org/html/2610.33209v1/simple_functions_data.png)Figure S4:Additional synthetic data and correspondingMAPSscores\. Exogenous variables are shown in the first and third rows, with the corresponding target series directly below\.Table S1:Predictability scores of bounded baseline methods on additional synthetic time series\. All scores lie in\[0,1\]\[0,1\], with higher values indicating greater estimated predictability\.To validate whetherMAPScaptures predictability across series with broader patterns, we construct additional synthetic data as shown in Figure[S4](https://arxiv.org/html/2610.33209#A3.F4), with detailed results reported in Table[S1](https://arxiv.org/html/2610.33209#A3.T1)\. We show rows of exogenous\-variable plots, each followed by the corresponding target\-variable plots\. Each dataset contains one or more exogenous variables, with the target sequence generated as a known function of the inputs\. These data include exponential, logarithmic, and periodic patterns such as the chirp, sawtooth, etc\. It can be observed thatMAPSassigns scores close to 1 across these deterministic systems\. In comparison, Spectral produces consistently high scores of approximately0\.920\.92, while ForeCA remains near0\.450\.45\. MI correctly identifies the monotonic deterministic functions but assigns substantially lower scores to the damped sine and chirp signals \(0\.270\.27and0\.560\.56, respectively\)\. Additionally, we include a constant series with sparse outlier spikes\. This series receives a highMAPSscore of0\.970\.97\. This result indicates that infrequent irregularities have a limited effect on the overall predictability of the series\. Since the underlying series remains largely constant, the outliers introduce limited uncertainty relative to the overall signal\.

Table S2:Subtraction\-correctedMAPSpredictability score analysis w\.r\.t\.different neighbourhood sizekk, averaged over2020locations\.
Table S3:Linear correlation betweenMAPSscore and Sonnet’s forecasting performance, taken across the2020regions within each setting\.

### C\.2Further analysis on the number of nearest neighbours \(kk\) in MAPS

We use theWEAdataset to assess how the neighbourhood sizekkaffects the predictability score\. Results are enumerated in Table[S3](https://arxiv.org/html/2610.33209#A3.T3)\. For each of the2020weather regions and99forecasting tasks, we evaluate the effect for values ofkkranging from11to all available neighbours and average the results across regions\. To compare results across different values ofkk, we remove the similarity that occurs by chance\. Specifically, we randomly shuffle the target vectors𝐲t\\mathbf\{y\}\_\{t\}across the input samples while keeping the input vectors𝐳t\\mathbf\{z\}\_\{t\}and their neighbourhoods𝒩t\\mathcal\{N\}\_\{t\}unchanged\. We then recompute the predictability score using the shuffled targets and report𝒫subtraction=𝒫⁡\(𝐗→𝐲\)−𝒫⁡\(𝐗→𝐲shuffled\)\\mathcal\{P\}\_\{\\mathrm\{subtraction\}\}=\\mathcal\{P\}\(\\mathbf\{X\}\\rightarrow\\mathbf\{y\}\)\-\\mathcal\{P\}\(\\mathbf\{X\}\\rightarrow\\mathbf\{y\}\_\{\\mathrm\{shuffled\}\}\)\. The shuffled score measures the neighbourhood similarity expected when the inputs and their futures are unrelated\. Subtracting it removes the chance similarity associated with each value ofkk\.

This subtraction is particularly important for small values ofkk\. Whenk=1k=1, each input window is evaluated only on the target vector of its nearest input neighbour\. Increasingkkprovides a more stable estimate of local future behaviour by combining similarity scores from multiple nearby inputs\. Atk=8k=8, the originalMAPSscore exceeds the permutation baseline by at least0\.350\.35across all forecasting tasks, showing that it captures substantially more structure than𝒫subtraction\\mathcal\{P\}\_\{\\mathrm\{subtraction\}\}\. Overall, intermediate values ofkkprovide the clearest predictability signal: overly small neighbourhoods offer limited evidence beyond chance, whereas very large neighbourhoods dilute the local relationship between inputs and their futures\.

We also examine how wellMAPSreflects forecasting performance across the same range ofkkin Table[S3](https://arxiv.org/html/2610.33209#A3.T3)\. Specifically, we compute the linear correlationrrbetweenMAPSand Sonnet’s forecasting performanceρ\\rho\. The correlation remains high across a broad range ofkk, including our default choicek=⌊log2⁡\(N\)⌋=12k=\\lfloor\\log\_\{2\}\(N\)\\rfloor=12\. Interestingly, the threeH=4H=4tasks reach the highest correlation atk=1k=1, likely because the target prediction of each input sample is compared only with that of its closest neighbour, preserving small differences in short\-term dynamics\. On the other hand, short\-horizon targets are also easier to predict, reflected in the forecasting performance\. Meanwhile, most longer\-horizon settings peak at moderate values ofkk\. This indicates that the default setting ofkkprovides a robust choice without requiringkkto be tuned separately for each task\.

### C\.3Results with additional metrics

![Refer to caption](https://arxiv.org/html/2610.33209v1/combined_weather_smape_sonnet.png)Figure S5:Relationship between different predictability measures and the target SMAPE across all evaluated locations under the Sonnet forecasting model\.Figure[S5](https://arxiv.org/html/2610.33209#A3.F5)compares the predictability scores with Sonnet’s SMAPE\. For the reasons discussed in Section[A\.1](https://arxiv.org/html/2610.33209#A1.SS1), we report these results for completeness rather than use SMAPE as the primary reference\. It can be observed that all methods correlate poorly with SMAPE, and stronger alignment with SMAPE does not necessarily indicate a more accurate predictability measure\.

### C\.4Results with additional machine learning model

To verify that the analysis ofMAPSis not limited to a single forecasting model, we further provide results with Crossformer[Zhang and Yan \[2023\]](https://arxiv.org/html/2610.33209#bib.bib33)in place of Sonnet\. Table[S4](https://arxiv.org/html/2610.33209#A3.T4)reports the linear correlationrracross2020locations between each predictability score and Crossformer’s performance, measured byρ\\rhoandss\. When comparing the predictability scores with Crossformer’s target correlationρ\\rho,MAPSyields the strongest correlationrrin88of the1818configurations\. It performs best in five of the six configurations withH=12H=12and three of the six configurations withH=28H=28\. Epiplexity is the most competitive baseline, leading all six configurations withH=4H=4and three configurations withH=28H=28, while ForeCA leads the remainingH=12H=12configuration\. Turning to the scaled MAE scoress,MAPSyields the strongest correlation in1313of the1818configurations, including five of the six configurations at bothH=12H=12andH=28H=28\. ForeCA leads the remaining five configurations\. These results show that the strong alignment betweenMAPSand forecasting performance is not limited to Sonnet\.

Table S4:Linear correlationrrbetween predictability scores and Crossformer’s forecasting performance, measured by target correlationρ\\rhoor scaled MAEss, across the2020WEAlocations\. Predictability scores and model performance are averaged over the 2016–2018 test periods before computingrr\. Results are reported across lookback windows \(LL\), forecasting horizons \(HH\), and covariate numbers \(CC\)\. The best results are bolded\.
### C\.5Additional analysis on predictability with univariate setup

Table S5:Linear correlationrrbetween target\-only predictability scores and Sonnet’s forecasting performance\. Arrows indicate whetherrrincreases or decreases relative to the corresponding result in Table[1](https://arxiv.org/html/2610.33209#S4.T1)and==indicates no change\.We further evaluate all methods in a univariate setting, where only the target series is used as input\. The resulting scores are compared with Sonnet’s performance usingC=45C=45andC=5C=5inputs\. This analysis therefore evaluates how much information about forecasting difficulty is lost when the predictability measures exclude the external covariates\. The results are presented in Table[S5](https://arxiv.org/html/2610.33209#A3.T5)\. Compared with the multivariable results in Table[1](https://arxiv.org/html/2610.33209#S4.T1), the correlation betweenMAPSandρ\\rhodecreases in all1818configurations\. This suggests that the exogenous covariates provide useful information about forecasting difficulty\. Spectral improves in every configuration, MI generally improves at longer horizons, PC performs worse, and ForeCA changes only slightly\.

Figure S6:Runtime analysis of predictability estimators\. We measure the median per\-window runtime \(ms\) versus the number of variables, lookback window length, and forecast horizon\.
### C\.6Runtime Comparison

Figure[S6](https://arxiv.org/html/2610.33209#A3.F6)compares the median per\-window runtime\. We change one factor at a time overC∈\{4,16,64,256,1024\}C\\in\\\{4,16,64,256,1024\\\},L∈\{7,14,28,56,112,224\}L\\in\\\{7,14,28,56,112,224\\\}, andH∈\{1,4,7,14,28,56,112\}H\\in\\\{1,4,7,14,28,56,112\\\}, holding the others atL=28L=28,H=7H=7, andC=64C=64\. Across the evaluated settings, the runtime ofMAPSdoes not change significantly because its dominant nearest\-neighbour search operates in the fixeddd\-dimensional PCA space\. By contrast, all baselines become more computationally expensive as the number of variables increases\. Additionally, MI also requires more runtime for longer lookback windows\.

## Appendix DLimitations

The main limitation is thatMAPSestimates predictability relative to the observed data distribution and should not be interpreted as a theoretical upper bound on forecasting performance\. As demonstrated under distribution shift, a score computed on training data may not characterise an unseen deployment distribution\. Furthermore,MAPSevaluates directional consistency and does not capture differences in forecast magnitude\. Finally, although we evaluate synthetic systems and real\-world weather and ILI tasks, further experiments across a broader range of domains are required to establish its general applicability\.