Paper deep dive
Beyond MSE: Rethinking the Evaluation Metric and Benchmarking for Irregular Time Series Forecasting
Rongwen Li, Haixin Xie, Xiao Wang, Changjian Chen
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/19/2026, 4:42:18 AM
Summary
The paper critiques the use of Mean Squared Error (MSE) for evaluating irregular time series forecasting, arguing that MSE is biased by sample-specific timestamp sampling distributions. It proposes Continuous-time Squared Error (CSE), an importance-weighted metric that estimates continuous-time risk more accurately. The authors provide theoretical proofs that CSE's asymptotic error is no greater than MSE's and introduce a systematic benchmark with synthetic, semi-synthetic, and real-world datasets to validate CSE's effectiveness.
Entities (7)
Relation Signals (6)
CSE → eliminatesinfluenceof → timestamp sampling distributions
confidence 95% · CSE, which employs importance weighting to eliminate the influence of the timestamp sampling distributions.
CSE → estimates → Continuous-time Risk
confidence 95% · CSE's asymptotic estimation error with respect to continuous-time risk is no greater than that of MSE.
MSE → isbiasedby → timestamp sampling distributions
confidence 95% · MSE is determined not only by the model prediction but also by the sample-specific timestamp sampling distributions, leading to a biased assessment
CSE → haslowerasymptoticerrorthan → MSE
confidence 90% · CSE’s asymptotic estimation error with respect to continuous-time risk is no greater than that of MSE.
Benchmark → validates → CSE
confidence 90% · we construct a systematic benchmark ... to validate the effectiveness of CSE
GSE → decomposesdiscrepancywith → MSE and CSE
confidence 85% · we further decompose the discrepancy between MSE and CSE into a sampling-dependence gap and a temporal-distribution gap... G_samp = MSE - GSE, G_time = GSE - CSE
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Existing research on irregular time-series forecasting has primarily focused on model design, while evaluation metrics remain insufficiently studied. Existing benchmarks typically use mean squared error (MSE) as the evaluation metric. We show that, in irregular forecasting, MSE is determined not only by the model prediction but also by the sample-specific timestamp sampling distributions, leading to a biased assessment of the models' continuous-time predictive performance. To address this issue, we propose the Continuous-time Squared Error (CSE), which employs importance weighting to eliminate the influence of the timestamp sampling distributions. We further theoretically prove that CSE's asymptotic estimation error with respect to continuous-time risk is no greater than that of MSE. Finally, we construct a systematic benchmark covering synthetic, semi-synthetic, and eight real-world datasets to validate the effectiveness of CSE and systematically evaluate models' continuous-time predictive performance. Experiments show that CSE can recover continuous-time risk more accurately than MSE, while relying solely on MSE may not fully reflect models' continuous-time predictive performance in real-world scenarios. Our code can be obtained at this https URL.
Tags
Links
- Source: https://arxiv.org/abs/2608.17293v1
- Canonical: https://arxiv.org/abs/2608.17293v1
Trouble viewing inline? Open PDF directly →
Full Text
100,767 characters extracted from source content.
Expand or collapse full text
Beyond MSE: Rethinking the Evaluation Metric and Benchmarking for Irregular Time Series Forecasting Rongwen Li Affiliation: College of Computer Science and Electronic Engineering, Hunan University, Changsha, Hunan, China Haixin Xie Affiliation: College of Computer Science and Electronic Engineering, Hunan University, Changsha, Hunan, China Xiao Wang Affiliation: College of Computer Science and Electronic Engineering, Hunan University, Changsha, Hunan, China Changjian Chen Affiliation: College of Computer Science and Electronic Engineering, Hunan University, Changsha, Hunan, China Correspondence to: changjianchen@hnu.edu.cn Abstract Existing research on irregular time-series forecasting has primarily focused on model design, while evaluation metrics remain insufficiently studied. Existing benchmarks typically use mean squared error (MSE) as the evaluation metric. We show that, in irregular forecasting, MSE is determined not only by the model prediction but also by the sample-specific timestamp sampling distributions, leading to a biased assessment of the models’ continuous-time predictive performance. To address this issue, we propose the Continuous-time Squared Error (CSE), which employs importance weighting to eliminate the influence of the timestamp sampling distributions. We further theoretically prove that CSE’s asymptotic estimation error with respect to continuous-time risk is no greater than that of MSE. Finally, we construct a systematic benchmark covering synthetic, semi-synthetic, and eight real-world datasets to validate the effectiveness of CSE and systematically evaluate models’ continuous-time predictive performance. Experiments show that CSE can recover continuous-time risk more accurately than MSE, while relying solely on MSE may not fully reflect models’ continuous-time predictive performance in real-world scenarios. Our code can be obtained at https://github.com/hnu-vis/ITS-Bench. Keywords: Machine Learning, ICML 1 Introduction Irregular time series are ubiquitous in real-world applications such as healthcare, environmental monitoring, and human activity analysis (7; 6; 1). Unlike regularly sampled time series, their observations occur at non-uniform timestamps, and different samples may exhibit different sampling patterns. In recent years, considerable research effort has been devoted to modeling such non-uniform temporal processes, leading to the development of continuous-time (2), attention-based (20), graph-based (23), and patch-based (15) models. However, compared with the rapid progress in model design, evaluation metrics for irregular time-series forecasting remain insufficiently studied. Figure 1: Illustration of evaluation bias under non-uniform sampling. Although Model B has lower continuous-time error, Model A obtains a lower MSE because of over-dense sampling within the interval [0.0,0.4][0.0,0.4]. For time-series forecasting, an ideal evaluation metric is the continuous-time discrepancy between the predicted trajectory and the true trajectory X over the entire future interval. Given X over a time span T, let ℓ(t,X) (t,X) denote the prediction error of a model at time t. This continuous-time discrepancy is defined as the continuous-time risk: Rct=X[1||∫ℓ(t,X)t].R_ct=E_X [ 1|T| _T (t,X)\,dt ]. (1) This risk directly characterizes the model’s continuous-time predictive performance over the future interval. Existing studies typically adopt mean squared error (MSE) and mean absolute error (MAE), which are widely used in regularly sampled forecasting. We argue that in irregular settings, MSE and MAE can be affected by the sample-specific timestamp sampling distributions and may therefore fail to accurately reflect a model’s performance. As illustrated in Fig. 1, although Model B has a lower continuous-time error than Model A (i.e., the orange line of Model B is closer to the black ground truth line over the whole time span), Model A obtains a lower MSE because of over-dense sampling within the interval [0.0,0.4][0.0,0.4]. By theoretical analysis, we find that MSE provides a biased estimate of the continuous-time risk RctR_ct, as it actually measures the observation-time risk Robs=X[∫ℓ(t,X)p(t∣X)t]R_obs=E_X[ _T (t,X)p(t X)\,dt], where p(t∣X)p(t X) is the conditional sampling distribution for the timestamps. This distribution is sample-specific and thus causes model performance evaluations to be confounded by sample-specific sampling patterns. To eliminate the influence of the timestamp sampling distributions on evaluation results, we begin directly from RctR_ct and construct its empirical estimator, termed the Continuous-time Squared Error (CSE), through importance sampling. CSE reweights errors according to the observation-time density p(t∣X)p(t X), reducing the contribution of densely sampled regions while increasing that of sparsely sampled regions, thereby mitigating evaluation bias caused by timestamp sampling distributions. Returning to Fig. 1, CSE correctly identifies Model B as having the lower continuous-time error. We further prove that the asymptotic estimation error of CSE with respect to the continuous-time risk is no greater than that of MSE. Finally, we further decompose the discrepancy between MSE and CSE into a sampling-dependence gap and a temporal-distribution gap to investigate their respective effects in real-world settings. Similarly, we can derive the Continuous-time Absolute Error (CAE) as the debiased counterpart to MAE. To evaluate the effectiveness of CSE, we require access to the true trajectory X of each dataset. To this end, we construct a systematic benchmark for irregular time-series forecasting covering synthetic, semi-synthetic, and real-world datasets. Both the synthetic and semi-synthetic datasets have the true continuous or regular trajectory. Based on a comprehensive evaluation of this benchmark, we have the following conclusions: (1) on synthetic and semi-synthetic datasets, CSE can recover continuous-time risk more accurately than MSE under non-uniform sampling; (2) on real-world datasets, relying solely on MSE may not fully reflect models’ continuous-time predictive performance; (3) the discrepancies between MSE and CSE are jointly influenced by the sampling-dependence and temporal-distribution gaps. The main contributions are summarized as follows: • Revealing the biased evaluation of MSE for model’s continuous-time predictive performance. We show that MSE estimates the observation-distribution risk rather than the continuous-time risk, thereby coupling model evaluation with sampling mechanisms. • A continuous-time evaluation metric with theoretical guarantees. We propose Continuous-Time Squared Error (CSE), which estimates the continuous-time risk from irregular timestamps through importance weighting, and prove that its asymptotic estimation error is no greater than MSE. • A systematic benchmark for continuous-time evaluation. We construct a benchmark covering synthetic, semi-synthetic, and real-world datasets to validate the effectiveness of CSE, systematically evaluate models’ continuous-time predictive performance, and analyze the sources of discrepancies between MSE and CSE. 2 Related Work 2.1 Irregular Time Series Forecasting Irregular time series are typically characterized by sparse, asynchronous, and nonuniform sampling. Existing studies have primarily focused on modeling continuous-time dynamics from limited observations. One line of work explicitly describes the continuous evolution of latent states. For example, Latent ODE, Neural CDE, and Continuous Recurrent Unit handle arbitrary time intervals based on ordinary differential equations, controlled differential equations, and continuous state transitions, respectively (18; 10; 19). Another line of work directly models irregular observations together with their temporal information. SeFT represents a time series as a set of observations, mTAN employs multi-time attention to aggregate asynchronous observations, and Raindrop uses graph structures to capture dynamic dependencies among variables (8; 20; 26). More recent methods further integrate multiscale representations with graph structures. For instance, Warpformer models temporal patterns at different scales through time warping, whereas GraFITi formulates irregular forecasting as an information propagation process over a graph (24; 23). Despite these advances, existing studies primarily focus on model design and typically evaluate predictions using MSE or MAE at observed future timestamps. In contrast, we study the evaluation protocol itself and examine how irregular sampling affects model comparison and selection. 2.2 Benchmarks and Evaluation for Irregular Time Series With the growing interest in irregular time series, prior studies have developed related benchmarks from the perspectives of data resources, software frameworks, and standardized experimental settings. The MIMIC-IV benchmark organizes unified tasks for sparse and irregular clinical time series and compares a range of representative models (3). PYRREGULAR provides unified data structures and analytical tools, together with a standardized benchmark covering multiple datasets and classification methods (21). Time-IMM further introduces new datasets and evaluation benchmarks for irregular multimodal multivariate time series (4). For forecasting tasks, Physiome-ODE constructs controlled datasets from biological ODEs to provide more challenging and discriminative environments for comparing irregular forecasting models (11). Existing benchmarks mainly improve datasets, tasks, and experimental standardization rather than the evaluation objective itself. Physiome-ODE is the most closely related work: it provides controlled ODE-generated trajectories to improve model discrimination, but still evaluates predictions using MSE at sampled query times (11). Consequently, it does not examine whether model conclusions depend on the query-time distribution or recover risk over the entire future interval. In contrast, we take continuous-time risk as the evaluation target, propose CSE to estimate it, and decompose the resulting evaluation discrepancy into sampling-dependence and temporal-distribution gaps. 3 Preliminaries 3.1 Irregular Time Series Forecasting The dataset contains S irregular time series. For the s-th sequence, its historical observations are denoted by s=(qs,i,s,i)i=1HsO_s=\(q_s,i,x_s,i)\_i=1^H_s, where qs,i∈ℝq_s,i and s,i∈ℝDx_s,i ^D denote the irregular timestamp and the corresponding multivariate observation, respectively. Given the historical observations sO_s, a forecasting model fθf_θ aims to estimate the latent process over a future time interval T. For any query time t∈t , the model output is given by ^s(t)=fθ(s,t) X_s(t)=f_θ(O_s,t). Let s(t)X_s(t) denote the ground-truth latent process, and define the prediction error at time t as ℓ(t,Xs)=ℒ(s(t),^s(t)) (t,X_s)=L(X_s(t), X_s(t)). In practical test datasets, the ground-truth values are available only at a finite set of future timestamps ts,jj=1Ls\t_s,j\_j=1^L_s, which are generated from the conditional sampling distribution p(t∣Xs)p(t X_s). 3.2 MSE and MAE under Irregular Sampling Existing studies on irregular time series forecasting typically compute MSE and MAE at the actually observed future timestamps. In the following, we use MSE as the example to illustrate the idea, i.e., ℓ(t,Xs)=‖s(t)−^s(t)‖2 (t,X_s)= \|X_s(t)- X_s(t) \|^2. Accordingly, the MSE on the test set can be written as MSE=1S∑s=1S1Ls∑j=1Lsℓ(ts,j,Xs).MSE= 1S _s=1^S 1L_s _j=1^L_s (t_s,j,X_s). (2) This metric first averages the prediction errors over the future observations of each trajectory and then assigns equal weight to different trajectories. When the future timestamps are regarded as conditional samples drawn from p(t∣X)p(t X), the above expression is essentially a Monte Carlo estimator of the following observation-distribution risk: Robs=X[∫ℓ(t,X)p(t∣X)t]=X[T|X[ℓ(T,X)]].R_obs=E_X [ _T (t,X)p(t X)\,dt ]=E_X [E_T X[ (T,X)] ]. (3) Therefore, MSE aggregates prediction errors according to the empirical frequency of future timestamps: densely observed temporal regions receive higher weights in the evaluation, whereas sparsely observed regions receive lower weights. In other words, the overall evaluation objective of MSE is determined not only by the model’s prediction errors but also by the conditional sampling distribution p(t∣X)p(t X). 4 Continuous-Time Risk for Irregular Forecasting In this section, we first define the continuous-time risk that directly corresponds to the continuous-time forecasting objective. We then construct its empirical estimator, CSE, using importance sampling and analyze the asymptotic properties of CSE relative to MSE. Finally, we introduce the global-time risk and decompose the evaluation discrepancy into two components associated with sample-dependent sampling and population-level temporal nonuniformity. 4.1 Continuous-Time Risk As discussed above, the observation-distribution risk corresponding to MSE is affected by the conditional sampling distribution p(t∣X)p(t X), thereby coupling model evaluation with sample-specific sampling mechanisms. A natural way to reduce this dependence is to adopt a common target-time measure across all trajectories. For continuous-time forecasting, an ideal evaluation objective should also measure the discrepancy between the predicted trajectory and the ground-truth process over the entire future interval. We therefore adopt the uniform temporal measure and define Rct=X[1||∫ℓ(t,X)t].R_ct=E_X [ 1|T| _T (t,X)\,dt ]. (4) RctR_ct assigns equal weight to each unit of time within the future interval and is independent of the sampling distribution of the actual observation times. It therefore characterizes the forecasting capability of a model over the entire continuous-time process. 4.2 Estimating Continuous-Time Risk The test set provides ground-truth observations only at a finite number of future timestamps, making it impossible to directly compute the continuous-time integral over the entire interval T. To estimate RctR_ct from the available irregular observations, we use the actual conditional sampling distribution p(t∣X)p(t X) as the proposal distribution and rewrite the continuous-time risk through importance sampling as Rct R_ct =X[∫ℓ(t,X)1||p(t∣X)p(t∣X)dt] =E_X [ _T (t,X) 1|T|p(t X)p(t X)\,dt ] (5) =X[T|X[ℓ(T,X)1||p(T∣X)]]. =E_X [E_T X [ (T,X) 1|T|p(T X) ] ]. Here, 1/(||p(t∣X))1/(|T|p(t X)) corrects the actual conditional sampling distribution toward the uniform temporal distribution. It reduces the contribution of each observation in densely sampled regions while increasing the contribution of each observation in sparsely sampled regions, allowing each temporal location to contribute to the evaluation according to the length of the time interval it represents. In practical datasets, the conditional sampling density p(t∣Xs)p(t X_s) is typically unknown. Let p^(ts,i∣Xs) p(t_s,i X_s) denote the conditional temporal density estimated from the future timestamps of the s-th trajectory. Since standard importance sampling may exhibit high finite-sample variance due to extreme weights, we adopt self-normalized importance sampling and define the resulting empirical metric under squared error as the Continuous-Time Squared Error (CSE): CSE=1S∑s=1S∑i=1Lsℓ(ts,i,Xs)/p^(ts,i∣Xs)∑i=1Ls1/p^(ts,i∣Xs).CSE= 1S _s=1^S _i=1^L_s (t_s,i,X_s)/ p(t_s,i X_s) _i=1^L_s1/ p(t_s,i X_s). (6) The constant 1/||1/|T| in the uniform temporal density cancels between the numerator and denominator. For each trajectory, a continuous-time error estimate is first obtained through self-normalized weighting, after which these estimates are averaged equally across all test trajectories. The same construction naturally yields the Continuous-time Absolute Error (CAE) under absolute error. In the following, we focus on CSE under squared error. To reduce the influence of an evaluation timestamp on its own density estimate, we employ leave-one-out Gaussian kernel density estimation: p^(ts,i∣Xs)=1(Ls−1)hs∑j=1j≠iLsK(ts,i−ts,jhs), p(t_s,i X_s)= 1(L_s-1)h_s _ subarraycj=1\\ j≠ i subarray^L_sK ( t_s,i-t_s,jh_s ), (7) where K(⋅)K(·) is the Gaussian kernel function and hsh_s is the bandwidth associated with the s-th trajectory. Excluding the current evaluation point prevents it from introducing additional kernel mass at its own location, thereby yielding a more robust estimate of the local temporal density. To further demonstrate the advantage of CSE over MSE, we theoretically compare their asymptotic estimation errors with respect to the continuous-time risk and introduce the following assumption and theorem. Assumption 4.1. Let T be a compact interval, and assume that different test trajectories are mutually independent. Conditional on XsX_s, the future timestamps are independently and identically distributed according to p(t∣Xs)p(t X_s). Suppose that the loss function is bounded, the conditional sampling density satisfies p(t∣X)≥ρ>0p(t X)≥ρ>0, and the density estimator p^(t∣X) p(t X) is strongly consistent over T. Theorem 4.2 (Asymptotic Non-Inferiority of CSE). Let L¯=min1≤s≤SLs L= _1≤ s≤ SL_s. Under Assumption 4.1, as both the number of test trajectories S and the minimum number of future observations per trajectory L¯ L tend to infinity, the following inequality holds almost surely: limS,L¯→∞|CSE−Rct|≤limS,L¯→∞|MSE−Rct|. _S, L→∞ |CSE-R_ct |≤ _S, L→∞ |MSE-R_ct |. (8) In particular, equality holds when p(t∣X)=1/||p(t X)=1/|T| almost everywhere. The complete proof is provided in the Appendix A. Theorem 4.2 shows that, under the stated conditions, the asymptotic estimation error of CSE with respect to the continuous-time risk is no greater than that of MSE. When future observation times are uniformly distributed, the two metrics are asymptotically equivalent. When the observation-distribution risk differs from the continuous-time risk, CSE achieves a strictly smaller asymptotic error. 4.3 Sources of Evaluation Discrepancy Because MSE and CSE converge to the observation-distribution risk and the continuous-time risk, respectively, they may yield different model evaluation results. To further analyze this discrepancy, we introduce the global-time risk and decompose the overall discrepancy into two components arising from sample-dependent sampling and population-level temporal nonuniformity. Marginalizing the conditional sampling distribution over trajectories yields pT(t)=X[p(t∣X)]p_T(t)=E_X[p(t X)], which describes the distribution of timestamps at the dataset level. Based on this distribution, we define the global-time risk as Rglobal=X[∫ℓ(t,X)pT(t)t].R_global=E_X [ _T (t,X)p_T(t)\,dt ]. (9) RglobalR_global removes the influence of sample-specific sampling mechanisms while retaining the population-level nonuniform temporal coverage of the dataset. In this sense, it preserves the global temporal distribution induced by the data collection process, but eliminates the instance-level irregularity that may bias model evaluation. It therefore serves as an intermediate reference between the observed risk RobsR_obs and the continuous-time risk RctR_ct. Based on this risk, the overall evaluation discrepancy can be decomposed as Robs−Rct R_obs-R_ct =Robs−Rglobal⏟Gsamp+Rglobal−Rct⏟Gtime, = R_obs-R_global_G_samp+ R_global-R_ct_G_time, (10) whereGsamp G_samp =∫CovX(ℓ(t,X),p(t∣X))dt, = _TCov_X ( (t,X),p(t X) )\,dt, Gtime G_time =CovU∼Unif()(ℓ¯(U),||pT(U)), =Cov_U (T) ( (U),|T|p_T(U) ), where ℓ¯(t)=X[ℓ(t,X)] (t)=E_X[ (t,X)]. The detailed derivation is provided in the Appendix B. GsampG_samp is determined by the covariance, across trajectories, between the model error ℓ(t,X) (t,X) and the sample-specific sampling density p(t∣X)p(t X) at a fixed temporal location. If trajectories with larger errors are more likely to be observed at that time, this term is positive; otherwise, it is negative. When the sampling process is independent of the trajectory, i.e., p(t∣X)=pT(t)p(t X)=p_T(t), we have Gsamp=0G_samp=0. Therefore, GsampG_samp characterizes the effect induced by sample-dependent sampling. GtimeG_time is determined by the covariance, across time, between the average model error ℓ¯(t) (t) and the population-level relative sampling density ||pT(t)|T|p_T(t). If temporal regions with larger average errors have higher population-level sampling densities, this term is positive; otherwise, it is negative. When pT(t)=1/||p_T(t)=1/|T|, we have Gtime=0G_time=0. Therefore, GtimeG_time characterizes the actual effect of population-level temporal nonuniformity on the evaluation of the current model. To estimate the above discrepancies in practice, we first estimate the marginal temporal density using p^T(t)=1S∑r=1Sp^(t∣Xr) p_T(t)= 1S _r=1^S p(t X_r). We then use p^T(t) p_T(t) as the target distribution to construct a self-normalized importance-sampling estimator of RglobalR_global, termed the Global-Time Squared Error (GSE). Its detailed form is provided in the Appendix C. The two discrepancy components can then be estimated as G^samp=MSE−GSE,G^time=GSE−CSE. G_samp=MSE-GSE, G_time=GSE-CSE. (11) This decomposition enables us to separately quantify the effects of sample-dependent sampling and population-level temporal nonuniformity on model scores and rankings in practice, providing a basis for analyzing evaluation discrepancies in the subsequent experiments. 5 Benchmark Protocol This section presents the benchmark setup, including the synthetic, semi-synthetic, and real-world datasets, representative baseline models, and unified training and evaluation protocols adopted to ensure fair and consistent comparisons. 5.1 Datasets Synthetic and Semi-Synthetic Data. To obtain directly computable continuous-time risks, we construct two fully synthetic datasets and two semi-synthetic datasets. Synthetic-Regime gradually transitions from a smooth, low-amplitude regime to a high-frequency periodic regime, simulating a process whose forecasting difficulty varies substantially over time. Synthetic-Multiscale simultaneously contains trends, multiscale periodic patterns, and local transient variations, covering more complex temporal patterns. Both fully synthetic datasets provide analytical continuous ground truth that can be queried at arbitrary future times. The semi-synthetic datasets are constructed from the regularly sampled ETTm1 and Weather datasets (22). Specifically, we sample a finite number of historical and future timestamps from complete regular windows to form irregular samples, and use the equally weighted error over the complete future grid as the reference continuous-time risk. Training target times uniformly cover the future interval, whereas test timestamps are sampled from pα(u)=αBeta(u,10,60)+(1−α)Uniform(0,1),p_α(u)= (u;10,60)+(1-α)Uniform(0,1), (12) where α∈0,0.3,0.5,0.9,0.99α∈\0,0.3,0.5,0.9,0.99\ controls the degree of temporal nonuniformity. All sampling conditions share the same training data and model predictions, while only the test-time distribution is varied, thereby isolating the effect of the evaluation protocol itself. Detailed data-generation procedures and window configurations are provided in the Appendix D.1. Real-World Datasets. We select eight real-world irregular time series datasets spanning different application domains and sampling patterns. USHCN contains long-term climate observations from weather stations across the United States (16). MIMIC-I consists of sparse clinical records from ICU patients (9). HumanActivity contains irregularly collected location-sensor measurements recorded during human activities (13). GDELT records multivariate event sequences triggered by international events (12). RepoHealth describes software repository development activities, whose sampling frequency varies with project activity (4). StudentLife contains mobile sensing data driven by user behavior and daily routines (17). FNSPID consists of stock prices and trading information and exhibits structured temporal gaps induced by trading sessions (4). CESNET contains traffic records generated by network devices, with timestamps affected by system scheduling delays and logging jitter (4). For all datasets, we use only numerical time series and follow the corresponding public settings for preprocessing and data splitting. Detailed dataset statistics, preprocessing procedures, and window configurations are provided in the Appendix D.2. 5.2 Baselines We select eleven representative irregular time series models covering different temporal modeling paradigms. These include recurrent or continuous-time models, namely GRU-D and NeuralFlow (5; 2); set-, attention-, or Transformer-based models, including SeFT, mTAN, and Warpformer (8; 20; 24); graph- or hypergraph-based models, including GraFITi, HyperIMTS, tPatchGNN, and ASTGI (23; 13; 25; 14); and pre-alignment and patch-aggregation methods, including KAFNet and APN (27; 15). Detailed descriptions and hyperparameter configurations for all models are provided in the Appendix D.3. 5.3 Training and Evaluation Protocol All models use the same data splits and input information, and their training objective is uniformly set to the MSE over observed timestamps. The remaining training and hyperparameter settings follow the corresponding public implementations, with model selection performed according to validation MSE. In each independent run, every model produces only one set of test predictions, on which MSE, GSE, and CSE are computed. This ensures that the differences among the metrics arise solely from the temporal weighting applied during test-time evaluation. Both the conditional and marginal temporal densities are estimated exclusively from the test timestamps and are shared across all models on the same dataset. We further compute G^samp G_samp and G^time G_time to analyze the sources of evaluation discrepancies. All experiments are independently repeated using five random seeds. The main text reports the mean results, while complete results with standard deviations are provided in the Appendix. Detailed training configurations, density-estimation parameters, and other implementation details are also provided in the Appendix D.4. 6 Experiments Using the proposed benchmark, we conduct experiments to investigate three questions: (1) Can CSE recover continuous-time risk and model rankings more accurately than MSE under non-uniform sampling? (2) Does observation-point MSE fully reflect models’ continuous-time predictive performance on real-world datasets? (3) How do the sampling-dependence and temporal-distribution gaps influence MSE–CSE discrepancies? 6.1 Can CSE Better Recover Continuous-Time Risk and Rankings than MSE? In this section, we evaluate the ability of CSE to recover continuous-time risks and model rankings in controlled environments with complete future ground truth. The experiments include two fully synthetic datasets and two semi-synthetic datasets, and consider five representative models: SeFT, GRU-D, GraFITi, tPatchGNN, and KAFNet. For the fully synthetic datasets, we compute RctR_ct through dense temporal integration. For the semi-synthetic datasets, we compute the reference risk over the complete regular future grid. The main text reports the results on Synthetic-Regime and Weather, while the remaining results are provided in the Appendix E.1. Table 1: Forecasting performance under MSE and CSE on eight real-world datasets. Superscripts denote ranks. The best and second-best results are highlighted in bold and underlined, respectively. CESNET ×10−1× 10^-1 FNSPID ×10−1× 10^-1 GDELT ×100× 10^0 HumanActivity ×10−2× 10^-2 MIMIC-I ×10−1× 10^-1 RepoHealth ×10−1× 10^-1 StudentLife ×10−1× 10^-1 USHCN ×10−1× 10^-1 Average Model MSE CSE MSE CSE MSE CSE MSE CSE MSE CSE MSE CSE MSE CSE MSE CSE Rank APN 8.76(5)8.76 (5) 8.58(5)8.58 (5) 2.38¯(2) 2.38 (2) 2.38¯(2) 2.38 (2) 1.02(1)1.02 (1) 0.99(1)0.99 (1) 5.68(6)5.68 (6) 5.65(6)5.65 (6) 6.66(6)6.66 (6) 6.72(6)6.72 (6) 3.36¯(2) 3.36 (2) 3.57(1)3.57 (1) 6.61¯(2) 6.61 (2) 6.68(3)6.68 (3) 6.04(4)6.04 (4) 6.71(5)6.71 (5) 3.6¯ 3.6 ASTGI 8.47(1)8.47 (1) 8.26(1)8.26 (1) 3.28(6)3.28 (6) 3.25(6)3.25 (6) 1.10(4)1.10 (4) 1.04(4)1.04 (4) 5.46(3)5.46 (3) 5.41(3)5.41 (3) 10.50(11)10.50 (11) 10.62(11)10.62 (11) 6.27(7)6.27 (7) 6.53(7)6.53 (7) 6.65(3)6.65 (3) 6.65¯(2) 6.65 (2) 5.14(3)5.14 (3) 5.66(3)5.66 (3) 4.74.7 GRU-D 10.22(10)10.22 (10) 9.91(10)9.91 (10) 5.24(8)5.24 (8) 5.20(8)5.20 (8) 1.22(10)1.22 (10) 1.15(10)1.15 (10) 8.56(9)8.56 (9) 8.56(9)8.56 (9) 6.67(7)6.67 (7) 6.75(7)6.75 (7) 11.68(11)11.68 (11) 11.95(11)11.95 (11) 7.21(8)7.21 (8) 7.23(8)7.23 (8) 7.78(11)7.78 (11) 8.59(11)8.59 (11) 9.39.3 GraFITi 8.50¯(2) 8.50 (2) 8.28¯(2) 8.28 (2) 3.21(4)3.21 (4) 3.22(4)3.22 (4) 1.09(3)1.09 (3) 1.01¯(2) 1.01 (2) 5.44¯(2) 5.44 (2) 5.40¯(2) 5.40 (2) 5.89(1)5.89 (1) 6.02¯(2) 6.02 (2) 3.32(1)3.32 (1) 3.61¯(2) 3.61 (2) 6.72(6)6.72 (6) 6.70(4)6.70 (4) 4.50(1)4.50 (1) 4.89(1)4.89 (1) 2.42.4 HyperIMTS 8.90(7)8.90 (7) 8.65(7)8.65 (7) 2.72(3)2.72 (3) 2.71(3)2.71 (3) 1.18(8)1.18 (8) 1.09(8)1.09 (8) 5.42(1)5.42 (1) 5.38(1)5.38 (1) 5.94¯(2) 5.94 (2) 5.92(1)5.92 (1) 3.72(3)3.72 (3) 3.97(3)3.97 (3) 6.54(1)6.54 (1) 6.56(1)6.56 (1) 6.67(8)6.67 (8) 7.30(8)7.30 (8) 4.14.1 KAFNet 8.73(4)8.73 (4) 8.54(4)8.54 (4) 2.23(1)2.23 (1) 2.22(1)2.22 (1) 1.08¯(2) 1.08 (2) 1.02(3)1.02 (3) 5.62(4)5.62 (4) 5.58(4)5.58 (4) 6.41(3)6.41 (3) 6.47(3)6.47 (3) 5.72(6)5.72 (6) 5.98(6)5.98 (6) 6.66(4)6.66 (4) 6.72(5)6.72 (5) 7.14(9)7.14 (9) 7.82(9)7.82 (9) 4.34.3 NeuralFlow 9.98(8)9.98 (8) 9.69(8)9.69 (8) 7.33(10)7.33 (10) 7.29(10)7.29 (10) 1.18(9)1.18 (9) 1.09(9)1.09 (9) 17.53(11)17.53 (11) 17.51(11)17.51 (11) 7.05(8)7.05 (8) 7.19(9)7.19 (9) 10.13(10)10.13 (10) 10.44(10)10.44 (10) 7.83(10)7.83 (10) 7.86(10)7.86 (10) 6.07(5)6.07 (5) 6.74(6)6.74 (6) 9.09.0 SeFT 10.11(9)10.11 (9) 9.79(9)9.79 (9) 8.70(11)8.70 (11) 8.66(11)8.66 (11) 1.13(6)1.13 (6) 1.05(5)1.05 (5) 9.38(10)9.38 (10) 9.36(10)9.36 (10) 7.10(9)7.10 (9) 7.15(8)7.15 (8) 9.13(8)9.13 (8) 9.44(8)9.44 (8) 8.94(11)8.94 (11) 8.95(11)8.95 (11) 7.64(10)7.64 (10) 8.41(10)8.41 (10) 9.19.1 Warpformer 8.58(3)8.58 (3) 8.37(3)8.37 (3) 3.66(7)3.66 (7) 3.63(7)3.63 (7) 1.30(11)1.30 (11) 1.20(11)1.20 (11) 5.65(5)5.65 (5) 5.61(5)5.61 (5) 6.42(4)6.42 (4) 6.50(5)6.50 (5) 5.04(5)5.04 (5) 5.33(5)5.33 (5) 6.78(7)6.78 (7) 6.80(7)6.80 (7) 6.13(6)6.13 (6) 6.67(4)6.67 (4) 5.95.9 mTAN 10.47(11)10.47 (11) 10.14(11)10.14 (11) 6.13(9)6.13 (9) 6.09(9)6.09 (9) 1.12(5)1.12 (5) 1.05(6)1.05 (6) 7.33(8)7.33 (8) 7.30(8)7.30 (8) 7.47(10)7.47 (10) 7.54(10)7.54 (10) 10.09(9)10.09 (9) 10.39(9)10.39 (9) 7.32(9)7.32 (9) 7.34(9)7.34 (9) 4.84¯(2) 4.84 (2) 5.23¯(2) 5.23 (2) 7.97.9 tPatchGNN 8.86(6)8.86 (6) 8.61(6)8.61 (6) 3.26(5)3.26 (5) 3.24(5)3.24 (5) 1.15(7)1.15 (7) 1.08(7)1.08 (7) 5.80(7)5.80 (7) 5.76(7)5.76 (7) 6.44(5)6.44 (5) 6.48(4)6.48 (4) 4.55(4)4.55 (4) 4.82(4)4.82 (4) 6.68(5)6.68 (5) 6.73(6)6.73 (6) 6.36(7)6.36 (7) 7.05(7)7.05 (7) 5.85.8 Fig. 4(a) and Fig. 4(b) show the estimation errors of MSE and CSE with respect to RctR_ct. Under uniform sampling, both metrics accurately recover the continuous-time risk. As α increases, the test-time distribution progressively deviates from the uniform distribution, causing the discrepancy between MSE and RctR_ct to increase consistently, whereas CSE maintains a substantially lower estimation error. For example, when α=0.9α=0.9, the relative errors of MSE and CSE on Synthetic-Regime are 80.6%80.6\% and 34.7%34.7\%, respectively. Under the more extreme setting of α=0.99α=0.99, although the error of CSE increases, it remains substantially lower than that of MSE. To further assess the recovery of model rankings, Table 2 compares the evaluation results of the five models under α=0.9α=0.9. On both Synthetic-Regime and Weather, CSE recovers exactly the same ranking of all five models as RctR_ct, whereas MSE produces a different ranking on both datasets. Overall, CSE effectively mitigates the influence of non-uniform test sampling on model evaluation. Across both synthetic and semi-synthetic datasets, CSE recovers continuous-time risk and model rankings more accurately than MSE, validating its effectiveness as a continuous-time evaluation metric. Table 2: Model evaluation under non-uniform test sampling with α=0.9α=0.9. Superscripts denote ranks; the best and second-best results are highlighted in bold and underlined, respectively. Synthetic-Regime ×10−2× 10^-2 Weather ×10−1× 10^-1 Model MSE CSE RctR_ct MSE CSE RctR_ct SeFT 1.09(3)1.09 (3) 8.70(1)8.70 (1) 6.26(1)6.26 (1) 4.40(5)4.40 (5) 4.96(5)4.96 (5) 4.85(5)4.85 (5) GRU-D 0.65(1)0.65 (1) 8.81¯(2) 8.81 (2) 6.48¯(2) 6.48 (2) 3.08(1)3.08 (1) 4.61(3)4.61 (3) 4.46(3)4.46 (3) GraFITi 2.57(5)2.57 (5) 9.88(5)9.88 (5) 7.66(5)7.66 (5) 3.89(4)3.89 (4) 4.92(4)4.92 (4) 4.76(4)4.76 (4) tPatchGNN 0.91¯(2) 0.91 (2) 9.38(4)9.38 (4) 6.91(4)6.91 (4) 3.88(3)3.88 (3) 3.69(1)3.69 (1) 3.53(1)3.53 (1) KAFNet 1.54(4)1.54 (4) 9.04(3)9.04 (3) 6.75(3)6.75 (3) 3.48¯(2) 3.48 (2) 3.78¯(2) 3.78 (2) 3.65¯(2) 3.65 (2) (a) Sampling-dependence gap G^samp G_samp (b) Temporal-distribution gap G^time G_time Figure 2: Normalized decomposition of the evaluation discrepancy into the sampling-dependence gap G^samp G_samp and temporal-distribution gap G^time G_time across models and datasets. 6.2 Does MSE Fully Reflect Continuous-Time Predictive Performance on Real-World Data? The previous section shows that CSE more accurately recovers continuous-time risk and model rankings when the underlying future trajectory is available. We now examine whether observation-point MSE provides a sufficient evaluation of models’ continuous-time predictive performance on real-world datasets, where the complete continuous trajectory is unavailable. Table 1 reports the results of 11 models on eight datasets, while Fig. 3 summarizes pairwise rank inversions and changes in the selection of the best-performing model. Comparisons between MAE and CAE, together with complete results including standard deviations, are provided in the Appendix E.2. The comparison shows that CSE is not a fixed-direction or fixed-ratio correction of MSE. The direction and magnitude of score changes vary across both datasets and models, indicating that the two metrics cannot be converted through simple rescaling. Moreover, five of the eight datasets exhibit pairwise rank inversions when replacing MSE with CSE. MIMIC-I and StudentLife each contain three inversions, GDELT and USHCN each contain two, and RepoHealth contains one. The best-performing model also changes on MIMIC-I and RepoHealth. The effects further differ across models. GraFITi ranks first on three datasets under MSE but on only one dataset under CSE. In contrast, the number of first-place rankings increases from two to three for HyperIMTS and from one to two for APN. Although GraFITi retains the best average rank, MSE and CSE do not provide fully consistent conclusions regarding the relative advantages of different models. Overall, the substantial discrepancies in model scores, rankings, and model selection suggest that relying solely on observation-point MSE may not fully reflect models’ continuous-time predictive performance on real-world datasets. (a) Pairwise rank inversions across datasets (b) Changes in Top-1 model selection Figure 3: Model ranking changes when replacing MSE with CSE on eight real-world datasets. Panel (a) reports the number of pairwise rank inversions, while Panel (b) compares the frequency with which each model is ranked first. 6.3 How Do the Sampling-Dependence and Temporal-Distribution Gaps Influence MSE–CSE Discrepancies? The previous section demonstrates that MSE and CSE can produce different model scores, rankings, and model selections on real-world datasets. We now use the risk decomposition introduced above to analyze the sources of these discrepancies. Fig. 2 reports the normalized sampling-dependence gap, G^samp=(MSE−GSE)/MSE G_samp=(MSE-GSE)/MSE, and the normalized temporal-distribution gap, G^time=(GSE−CSE)/GSE G_time=(GSE-CSE)/GSE. The relative contributions of the two gaps vary substantially across datasets. GDELT and RepoHealth are primarily affected by the sampling-dependence gap, whose cross-model averages are approximately 5.9%5.9\% and −5.8%-5.8\%, respectively. CESNET is mainly affected by the temporal-distribution gap, with an average value of approximately 3.3%3.3\%. On USHCN, the sampling-dependence and temporal-distribution gaps are approximately −7.8%-7.8\% and −2.2%-2.2\%, respectively, indicating that both components contribute to the discrepancy between MSE and CSE. The two gaps are comparatively small on the remaining datasets. Within the same dataset, models generally exhibit similar signs for each discrepancy component, suggesting that the overall direction is largely determined by the dataset’s sampling structure. However, the magnitude varies across models because of their different temporal error patterns. This model-dependent effect can alter model rankings when the original performance differences are small. For example, on MIMIC-I, the sampling-dependence gaps of GraFITi and HyperIMTS are −2.2%-2.2\% and +0.8%+0.8\%, respectively, changing the best-performing model from GraFITi under MSE to HyperIMTS under CSE. On RepoHealth, GraFITi exhibits a more negative sampling-dependence gap than APN, giving it a larger apparent advantage under MSE; this advantage disappears under CSE, making APN the best-performing model. (a) Synthetic-Regime (b) Weather Figure 4: Relative errors of MSE and CSE with respect to RctR_ct under increasing sampling non-uniformity. Overall, the discrepancies between MSE and CSE are jointly influenced by the sampling-dependence and temporal-distribution gaps, whose relative contributions vary across datasets and models. 7 Conclusion This work revisits the evaluation of irregular time series forecasting from a continuous-time perspective. We show that observation-point MSE and MAE effectively estimate the observation-distribution risk induced by the sampling distribution, rather than the continuous-time risk over the entire future interval. To mitigate this objective mismatch, we propose the Continuous-Time Squared Error (CSE), which estimates the continuous-time risk through self-normalized importance weighting, and prove that its asymptotic estimation error is no greater than that of MSE. We further decompose the discrepancy between the two risks into a sampling-dependence gap and a temporal-distribution gap, thereby characterizing how different sampling structures affect model evaluation. Finally, we construct a benchmark covering synthetic, semi-synthetic, and real-world datasets to validate CSE and systematically evaluate models’ continuous-time predictive performance. References Afrifa-Yamoah et al. (2020) E. Afrifa-Yamoah, U. A. Mueller, S. M. Taylor, and A. J. Fisher Missing data imputation of high-resolution temporal climate time series data. Meteorological Applications 27 (1), p. e1873. Cited by: §1. Bilos et al. (2021) M. Bilos, J. Sommer, S. S. Rangapuram, T. Januschowski, and S. Günnemann Neural flows: efficient alternative to neural ODEs. In Advances in Neural Information Processing Systems, p. 21325–21337. Cited by: §D.3, §1, §5.2. Bui et al. (2024) H. Bui, H. Warrier, and Y. Gupta Benchmarking with MIMIC-IV, an irregular, spare clinical time series dataset. arXiv preprint arXiv:2401.15290. External Links: Document Cited by: §2.2. Chang et al. (2025) C. Chang, J. Hwang, Y. Shi, H. Wang, W. Peng, T. Chen, and W. Wang Time-IMM: a dataset and benchmark for irregular multimodal multivariate time series. CoRR abs/2506.10412. Cited by: §D.2, §2.2, §5.1. Che et al. (2016) Z. Che, S. Purushotham, K. Cho, D. A. Sontag, and Y. Liu Recurrent neural networks for multivariate time series with missing values. CoRR abs/1606.01865. Cited by: §D.3, §5.2. Decorte et al. (2024) T. Decorte, S. Mortier, J. J. Lembrechts, F. J. R. Meysman, S. Latré, E. Mannens, and T. Verdonck Missing value imputation of wireless sensor data for environmental monitoring. Sensors 24 (8), p. 2416. Cited by: §1. Ghassemi et al. (2015) M. Ghassemi, M. A. F. Pimentel, T. Naumann, T. Brennan, D. A. Clifton, P. Szolovits, and M. Feng A multivariate timeseries modeling approach to severity of illness assessment and forecasting in ICU with sparse, heterogeneous clinical data. In Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence, January 25-30, 2015, Austin, Texas, USA, p. 446–453. Cited by: §1. Horn et al. (2020) M. Horn, M. Moor, C. Bock, B. Rieck, and K. M. Borgwardt Set functions for time series. In International Conference on Machine Learning, p. 4353–4363. Cited by: §D.3, §2.1, §5.2. Johnson et al. (2016) A. E. W. Johnson, T. J. Pollard, L. Shen, L. H. Lehman, M. Feng, M. Ghassemi, B. Moody, P. Szolovits, L. A. Celi, and R. G. Mark MIMIC-I, a freely accessible critical care database. Scientific Data 3, p. 160035. External Links: Document Cited by: §D.2, §5.1. Kidger et al. (2020) P. Kidger, J. Morrill, J. Foster, and T. J. Lyons Neural controlled differential equations for irregular time series. In Advances in Neural Information Processing Systems, Cited by: §2.1. Klötergens et al. (2025) C. Klötergens, V. K. Yalavarthi, R. Scholz, M. Stubbemann, S. Born, and L. Schmidt-Thieme Physiome-ODE: a benchmark for irregularly sampled multivariate time-series forecasting based on biological ODEs. In International Conference on Learning Representations, Cited by: §2.2, §2.2. Leetaru and Schrodt (2013) K. Leetaru and P. A. Schrodt GDELT: global data on events, location, and tone. ISA Annual Convention. Cited by: §5.1. Li et al. (2025) B. Li, Y. Luo, Z. Liu, J. Zheng, J. Lv, and Q. Ma HyperIMTS: hypergraph neural network for irregular multivariate time series forecasting. In International Conference on Machine Learning, p. 35502–35518. Cited by: §D.3, §5.1, §5.2. Liu et al. (2026a) X. Liu, X. Qiu, H. Cheng, X. Wu, C. Guo, B. Yang, and J. Hu ASTGI: adaptive spatio-temporal graph interactions for irregular multivariate time series forecasting. In International Conference on Learning Representations, Cited by: §D.3, §5.2. Liu et al. (2026b) X. Liu, X. Qiu, X. Wu, Z. Li, C. Guo, J. Hu, and B. Yang Rethinking irregular time series forecasting: a simple yet effective baseline. In AAAI Conference on Artificial Intelligence, p. 23873–23881. Cited by: §D.3, §1, §5.2. Menne et al. (2015) M. J. Menne, C. N. Jr. Williams, and R. S. Vose Long-term daily and monthly climate records from stations across the contiguous united states (u.s. historical climatology network). Note: Environmental System Science Data Infrastructure for a Virtual Ecosystem (ESS-DIVE) Cited by: §D.2, §5.1. Nepal et al. (2024) S. Nepal, W. Liu, A. Pillai, W. Wang, V. Vojdanovski, J. F. Huckins, C. Rogers, M. L. Meyer, and A. T. Campbell Capturing the college experience: a four-year mobile sensing study of mental health, resilience and behavior of college students during the pandemic. Proceedings of the ACM on interactive, mobile, wearable and ubiquitous technologies 8 (1), p. 1–37. Cited by: §5.1. Rubanova et al. (2019) Y. Rubanova, T. Q. Chen, and D. Duvenaud Latent ordinary differential equations for irregularly-sampled time series. In Advances in Neural Information Processing Systems, p. 5321–5331. Cited by: §2.1. Schirmer et al. (2022) M. Schirmer, M. Eltayeb, S. Lessmann, and M. Rudolph Modeling irregular time series with continuous recurrent units. In International Conference on Machine Learning, p. 19388–19405. Cited by: §D.3, §2.1. Shukla and Marlin (2021) S. N. Shukla and B. M. Marlin Multi-time attention networks for irregularly sampled time series. In International Conference on Learning Representations, Cited by: §D.3, §1, §2.1, §5.2. Spinnato and Landi (2025) F. Spinnato and C. Landi PYRREGULAR: a unified framework for irregular time series, with classification benchmarks. CoRR abs/2505.06047. Cited by: §2.2. Wu et al. (2021) H. Wu, J. Xu, J. Wang, and M. Long Autoformer: decomposition transformers with auto-correlation for long-term series forecasting. In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, p. 22419–22430. Cited by: §5.1. Yalavarthi et al. (2024) V. K. Yalavarthi, K. Madhusudhanan, R. Scholz, N. Ahmed, J. Burchert, S. Jawed, S. Born, and L. Schmidt-Thieme GraFITi: graphs for forecasting irregularly sampled time series. In AAAI Conference on Artificial Intelligence, p. 16255–16263. Cited by: §D.3, §1, §2.1, §5.2. Zhang et al. (2023) J. Zhang, S. Zheng, W. Cao, J. Bian, and J. Li Warpformer: a multi-scale modeling approach for irregular clinical time series. In ACM SIGKDD Conference on Knowledge Discovery and Data Mining, p. 3273–3285. Cited by: §D.3, §2.1, §5.2. Zhang et al. (2024) W. Zhang, C. Yin, H. Liu, X. Zhou, and H. Xiong Irregular multivariate time series forecasting: a transformable patching graph neural networks approach. In International Conference on Machine Learning, p. 60179–60196. Cited by: §D.3, §5.2. Zhang et al. (2022) X. Zhang, M. Zeman, T. Tsiligkaridis, and M. Zitnik Graph-guided network for irregularly sampled multivariate time series. In International Conference on Learning Representations, Cited by: §2.1. Zhou et al. (2026) Z. Zhou, Y. Huang, Y. Wang, Y. Wu, J. Kwok, and Y. Liang Revitalizing canonical pre-alignment for irregular multivariate time series forecasting. In AAAI Conference on Artificial Intelligence, p. 29115–29123. Cited by: §D.3, §5.2. Appendix A Proof of Theorem 1 For the s-th trajectory, denote ps(t)=p(t∣Xs),ℓs(t)=ℓ(t,Xs),p_s(t)=p(t X_s), _s(t)= (t,X_s), (13) and define rct(Xs) r_ct(X_s) =1||∫ℓs(t)t, = 1|T| _T _s(t)\,dt, (14) robs(Xs) r_obs(X_s) =∫ℓs(t)ps(t)t. = _T _s(t)p_s(t)\,dt. (15) Accordingly, Rct=X[rct(X)],Robs=X[robs(X)].R_ct=E_X[r_ct(X)], R_obs=E_X[r_obs(X)]. (16) We first introduce the oracle importance-sampling terms As=1Ls∑i=1Lsℓs(ts,i)ps(ts,i),Bs=1Ls∑i=1Ls1ps(ts,i),A_s= 1L_s _i=1^L_s _s(t_s,i)p_s(t_s,i), B_s= 1L_s _i=1^L_s 1p_s(t_s,i), (17) and their estimated counterparts A^s=1Ls∑i=1Lsℓs(ts,i)p^(ts,i∣Xs),B^s=1Ls∑i=1Ls1p^(ts,i∣Xs). A_s= 1L_s _i=1^L_s _s(t_s,i) p(t_s,i X_s), B_s= 1L_s _i=1^L_s 1 p(t_s,i X_s). (18) By Assumption 1, ps(t)≥ρ>0p_s(t)≥ρ>0 and p^(t∣Xs) p(t X_s) is uniformly consistent on T. Hence, almost surely, for sufficiently large LsL_s, inft∈p^(t∣Xs)≥ρ2, _t p(t X_s)≥ ρ2, (19) and therefore supt∈|1p^(t∣Xs)−1ps(t)| _t | 1 p(t X_s)- 1p_s(t) | ≤2ρ2supt∈|p^(t∣Xs)−ps(t)| ≤ 2ρ^2 _t | p(t X_s)-p_s(t) | (20) →a.s.0. a.s.0. Let M<∞M<∞ satisfy 0≤ℓ(t,X)≤M0≤ (t,X)≤ M. It follows that |A^s−As| | A_s-A_s | ≤Msupt∈|1p^(t∣Xs)−1ps(t)|→a.s.0, ≤ M _t | 1 p(t X_s)- 1p_s(t) | a.s.0, (21) |B^s−Bs| | B_s-B_s | ≤supt∈|1p^(t∣Xs)−1ps(t)|→a.s.0. ≤ _t | 1 p(t X_s)- 1p_s(t) | a.s.0. (22) Conditional on XsX_s, the future timestamps are independently and identically distributed according to ps(t)p_s(t). The strong law of large numbers therefore gives As A_s →a.s.T|Xs[ℓs(T)ps(T)]=∫ℓs(t)t, a.s.E_T X_s [ _s(T)p_s(T) ]= _T _s(t)\,dt, (23) Bs B_s →a.s.T|Xs[1ps(T)]=||. a.s.E_T X_s [ 1p_s(T) ]=|T|. (24) Combining the above results yields A^s→a.s.∫ℓs(t)t,B^s→a.s.||. A_s a.s. _T _s(t)\,dt, B_s a.s.|T|. (25) Define the trajectory-level CSE estimator as r^ct,s=∑i=1Lsℓs(ts,i)/p^(ts,i∣Xs)∑i=1Ls1/p^(ts,i∣Xs)=A^sB^s. r_ct,s= _i=1^L_s _s(t_s,i)/ p(t_s,i X_s) _i=1^L_s1/ p(t_s,i X_s)= A_s B_s. (26) Since ||>0|T|>0, the continuous mapping theorem gives r^ct,s→a.s.1||∫ℓs(t)t=rct(Xs). r_ct,s a.s. 1|T| _T _s(t)\,dt=r_ct(X_s). (27) For MSE, the conditional strong law directly gives r^obs,s=1Ls∑i=1Lsℓs(ts,i)→a.s.T|Xs[ℓs(T)]=robs(Xs). r_obs,s= 1L_s _i=1^L_s _s(t_s,i) a.s.E_T X_s[ _s(T)]=r_obs(X_s). (28) Let L¯=min1≤s≤SLs L= _1≤ s≤ SL_s. For the joint limit below, Assumption 1 is understood to hold uniformly over the test trajectories. Consequently, Eqs. (27) and (28) imply 1S∑s=1S|r^ct,s−rct(Xs)| 1S _s=1^S | r_ct,s-r_ct(X_s) | →a.s.0, a.s.0, (29) 1S∑s=1S|r^obs,s−robs(Xs)| 1S _s=1^S | r_obs,s-r_obs(X_s) | →a.s.0 a.s.0 (30) as L¯→∞ L→∞. Since the test trajectories are mutually independent and the trajectory-level risks are bounded, another application of the strong law yields 1S∑s=1Srct(Xs)→a.s.Rct,1S∑s=1Srobs(Xs)→a.s.Robs. 1S _s=1^Sr_ct(X_s) a.s.R_ct, 1S _s=1^Sr_obs(X_s) a.s.R_obs. (31) Table 3: Complete MSE and CSE results on CESNET, FNSPID, GDELT, and HumanActivity. Values are mean ± standard deviation over three random seeds; standard deviations are rounded to three decimal places. Model CESNET ×10−1× 10^-1 FNSPID ×10−1× 10^-1 GDELT ×100× 10^0 HumanActivity ×10−2× 10^-2 MSE CSE MSE CSE MSE CSE MSE CSE APN ±0.007(5)8.76\!±\!0.007^(5) ±0.008(5)8.58\!±\!0.008^(5) 2.38¯±0.143(2) 2.38\!±\!0.143^(2) 2.38¯±0.144(2) 2.38\!±\!0.144^(2) ±0.001(1)1.02\!±\!0.001^(1) ±0.001(1)0.99\!±\!0.001^(1) ±0.035(6)5.68\!±\!0.035^(6) ±0.039(6)5.65\!±\!0.039^(6) ASTGI ±0.137(1)8.47\!±\!0.137^(1) ±0.122(1)8.26\!±\!0.122^(1) ±0.213(6)3.28\!±\!0.213^(6) ±0.211(6)3.25\!±\!0.211^(6) ±0.003(4)1.10\!±\!0.003^(4) ±0.003(4)1.04\!±\!0.003^(4) ±0.167(3)5.46\!±\!0.167^(3) ±0.167(3)5.41\!±\!0.167^(3) GRU-D ±0.464(10)10.22\!±\!0.464^(10) ±0.466(10)9.91\!±\!0.466^(10) ±0.800(8)5.24\!±\!0.800^(8) ±0.791(8)5.20\!±\!0.791^(8) ±0.003(10)1.22\!±\!0.003^(10) ±0.003(10)1.15\!±\!0.003^(10) ±1.497(9)8.56\!±\!1.497^(9) ±1.515(9)8.56\!±\!1.515^(9) GraFITi 8.50¯±0.113(2) 8.50\!±\!0.113^(2) 8.28¯±0.102(2) 8.28\!±\!0.102^(2) ±0.755(4)3.21\!±\!0.755^(4) ±0.760(4)3.22\!±\!0.760^(4) ±0.002(3)1.09\!±\!0.002^(3) 1.01¯±0.002(2) 1.01\!±\!0.002^(2) 5.44¯±0.038(2) 5.44\!±\!0.038^(2) 5.40¯±0.034(2) 5.40\!±\!0.034^(2) HyperIMTS ±0.073(7)8.90\!±\!0.073^(7) ±0.060(7)8.65\!±\!0.060^(7) ±0.593(3)2.72\!±\!0.593^(3) ±0.592(3)2.71\!±\!0.592^(3) ±0.002(8)1.18\!±\!0.002^(8) ±0.002(8)1.09\!±\!0.002^(8) ±0.037(1)5.42\!±\!0.037^(1) ±0.039(1)5.38\!±\!0.039^(1) KAFNet ±0.012(4)8.73\!±\!0.012^(4) ±0.003(4)8.54\!±\!0.003^(4) ±0.047(1)2.23\!±\!0.047^(1) ±0.043(1)2.22\!±\!0.043^(1) 1.08¯±0.005(2) 1.08\!±\!0.005^(2) ±0.005(3)1.02\!±\!0.005^(3) ±0.021(4)5.62\!±\!0.021^(4) ±0.025(4)5.58\!±\!0.025^(4) NeuralFlow ±0.103(8)9.98\!±\!0.103^(8) ±0.118(8)9.69\!±\!0.118^(8) ±0.138(10)7.33\!±\!0.138^(10) ±0.140(10)7.29\!±\!0.140^(10) ±0.005(9)1.18\!±\!0.005^(9) ±0.005(9)1.09\!±\!0.005^(9) ±2.512(11)17.53\!±\!2.512^(11) ±2.523(11)17.51\!±\!2.523^(11) SeFT ±0.183(9)10.11\!±\!0.183^(9) ±0.177(9)9.79\!±\!0.177^(9) ±0.958(11)8.70\!±\!0.958^(11) ±0.957(11)8.66\!±\!0.957^(11) ±0.003(6)1.13\!±\!0.003^(6) ±0.003(5)1.05\!±\!0.003^(5) ±0.422(10)9.38\!±\!0.422^(10) ±0.406(10)9.36\!±\!0.406^(10) Warpformer ±0.067(3)8.58\!±\!0.067^(3) ±0.058(3)8.37\!±\!0.058^(3) ±0.216(7)3.66\!±\!0.216^(7) ±0.220(7)3.63\!±\!0.220^(7) ±0.016(11)1.30\!±\!0.016^(11) ±0.016(11)1.20\!±\!0.016^(11) ±0.143(5)5.65\!±\!0.143^(5) ±0.141(5)5.61\!±\!0.141^(5) mTAN ±0.153(11)10.47\!±\!0.153^(11) ±0.106(11)10.14\!±\!0.106^(11) ±0.743(9)6.13\!±\!0.743^(9) ±0.738(9)6.09\!±\!0.738^(9) ±0.003(5)1.12\!±\!0.003^(5) ±0.003(6)1.05\!±\!0.003^(6) ±0.108(8)7.33\!±\!0.108^(8) ±0.107(8)7.30\!±\!0.107^(8) tPatchGNN ±0.140(6)8.86\!±\!0.140^(6) ±0.114(6)8.61\!±\!0.114^(6) ±0.269(5)3.26\!±\!0.269^(5) ±0.269(5)3.24\!±\!0.269^(5) ±0.006(7)1.15\!±\!0.006^(7) ±0.006(7)1.08\!±\!0.006^(7) ±0.121(7)5.80\!±\!0.121^(7) ±0.132(7)5.76\!±\!0.132^(7) Therefore, as S,L¯→∞S, L→∞, CSE→a.s.Rct,MSE→a.s.Robs.CSE a.s.R_ct, a.s.R_obs. (32) It follows that limS,L¯→∞|CSE−Rct| _S, L→∞ |CSE-R_ct | =0, =0, (33) limS,L¯→∞|MSE−Rct| _S, L→∞ |MSE-R_ct | =|Robs−Rct| = |R_obs-R_ct | (34) almost surely. Hence, limS,L¯→∞|CSE−Rct|≤limS,L¯→∞|MSE−Rct|. _S, L→∞ |CSE-R_ct |≤ _S, L→∞ |MSE-R_ct |. (35) When p(t∣X)=1/||p(t X)=1/|T| almost everywhere, we have Robs=RctR_obs=R_ct, and both asymptotic estimation errors are zero. More generally, equality holds whenever Robs=RctR_obs=R_ct, whereas the inequality is strict when Robs≠RctR_obs≠ R_ct. The self-normalized estimator may be biased at a finite sample size; the result above concerns its consistency and asymptotic estimation error. Appendix B Derivation of the Evaluation Discrepancy Recall that pT(t)=X[p(t∣X)],ℓ¯(t)=X[ℓ(t,X)].p_T(t)=E_X[p(t X)], (t)=E_X[ (t,X)]. (36) For the sampling-dependence gap, Gsamp G_samp =Robs−Rglobal =R_obs-R_global =X[∫ℓ(t,X)p(t∣X)t]−X[∫ℓ(t,X)pT(t)t] =E_X [ _T (t,X)p(t X)\,dt ]-E_X [ _T (t,X)p_T(t)\,dt ] =∫X[ℓ(t,X)p(t∣X)]−X[ℓ(t,X)]X[p(t∣X)]t = _T \E_X[ (t,X)p(t X)]-E_X[ (t,X)]E_X[p(t X)] \dt =∫CovX(ℓ(t,X),p(t∣X))t. = _TCov_X ( (t,X),p(t X) )dt. (37) For the temporal-distribution gap, let U∼Unif()U (T). Since U[||pT(U)]=1E_U[|T|p_T(U)]=1, we have Gtime G_time =Rglobal−Rct =R_global-R_ct =∫ℓ¯(t)pT(t)t−1||∫ℓ¯(t)t = _T (t)p_T(t)\,dt- 1|T| _T (t)\,dt =U[ℓ¯(U)||pT(U)]−U[ℓ¯(U)]U[||pT(U)] =E_U [ (U)|T|p_T(U) ]-E_U[ (U)]E_U[|T|p_T(U)] =CovU∼Unif()(ℓ¯(U),||pT(U)). =Cov_U (T) ( (U),|T|p_T(U) ). (38) Therefore, Robs−Rct=Gsamp+Gtime.R_obs-R_ct=G_samp+G_time. (39) Here, GsampG_samp captures the effect of sample-dependent sampling, whereas GtimeG_time captures the effect of population-level temporal non-uniformity. Appendix C Global-Time Squared Error The marginal timestamp density is estimated by averaging the trajectory-specific conditional densities: p^T(t)=1S∑r=1Sp^(t∣Xr). p_T(t)= 1S _r=1^S p(t X_r). (40) Using p^(t∣Xs) p(t X_s) as the proposal density and p^T(t) p_T(t) as the target density, we define the Global-Time Squared Error as GSE=1S∑s=1S∑i=1Lsℓ(ts,i,Xs)p^T(ts,i)p^(ts,i∣Xs)∑i=1Lsp^T(ts,i)p^(ts,i∣Xs).GSE= 1S _s=1^S _i=1^L_s (t_s,i,X_s) p_T(t_s,i) p(t_s,i X_s) _i=1^L_s p_T(t_s,i) p(t_s,i X_s). (41) GSE evaluates each trajectory under the dataset-level marginal timestamp distribution, thereby removing sample-specific sampling effects while retaining population-level temporal non-uniformity. It serves as the empirical counterpart of RglobalR_global and an intermediate reference between MSE and CSE. Table 4: Complete MSE and CSE results on MIMIC-I, RepoHealth, StudentLife, and USHCN. Values are mean ± standard deviation over three random seeds; standard deviations are rounded to three decimal places. Model MIMIC-I ×10−1× 10^-1 RepoHealth ×10−1× 10^-1 StudentLife ×10−1× 10^-1 USHCN ×10−1× 10^-1 MSE CSE MSE CSE MSE CSE MSE CSE APN ±0.038(6)6.66\!±\!0.038^(6) ±0.038(6)6.72\!±\!0.038^(6) 3.36¯±0.069(2) 3.36\!±\!0.069^(2) ±0.080(1)3.57\!±\!0.080^(1) 6.61¯±0.009(2) 6.61\!±\!0.009^(2) ±0.006(3)6.68\!±\!0.006^(3) ±0.554(4)6.04\!±\!0.554^(4) ±0.686(5)6.71\!±\!0.686^(5) ASTGI ±0.000(11)10.50\!±\!0.000^(11) ±0.000(11)10.62\!±\!0.000^(11) ±2.310(7)6.27\!±\!2.310^(7) ±2.378(7)6.53\!±\!2.378^(7) ±0.036(3)6.65\!±\!0.036^(3) 6.65¯±0.038(2) 6.65\!±\!0.038^(2) ±0.169(3)5.14\!±\!0.169^(3) ±0.179(3)5.66\!±\!0.179^(3) GRU-D ±0.043(7)6.67\!±\!0.043^(7) ±0.046(7)6.75\!±\!0.046^(7) ±0.695(11)11.68\!±\!0.695^(11) ±0.672(11)11.95\!±\!0.672^(11) ±0.045(8)7.21\!±\!0.045^(8) ±0.045(8)7.23\!±\!0.045^(8) ±0.690(11)7.78\!±\!0.690^(11) ±0.756(11)8.59\!±\!0.756^(11) GraFITi ±0.264(1)5.89\!±\!0.264^(1) 6.02¯±0.278(2) 6.02\!±\!0.278^(2) ±0.341(1)3.32\!±\!0.341^(1) 3.61¯±0.342(2) 3.61\!±\!0.342^(2) ±0.164(6)6.72\!±\!0.164^(6) ±0.162(4)6.70\!±\!0.162^(4) ±0.252(1)4.50\!±\!0.252^(1) ±0.299(1)4.89\!±\!0.299^(1) HyperIMTS 5.94¯±0.068(2) 5.94\!±\!0.068^(2) ±0.068(1)5.92\!±\!0.068^(1) ±0.437(3)3.72\!±\!0.437^(3) ±0.444(3)3.97\!±\!0.444^(3) ±0.067(1)6.54\!±\!0.067^(1) ±0.062(1)6.56\!±\!0.062^(1) ±1.837(8)6.67\!±\!1.837^(8) ±2.073(8)7.30\!±\!2.073^(8) KAFNet ±0.021(3)6.41\!±\!0.021^(3) ±0.022(3)6.47\!±\!0.022^(3) ±3.034(6)5.72\!±\!3.034^(6) ±3.061(6)5.98\!±\!3.061^(6) ±0.023(4)6.66\!±\!0.023^(4) ±0.024(5)6.72\!±\!0.024^(5) ±0.418(9)7.14\!±\!0.418^(9) ±0.406(9)7.82\!±\!0.406^(9) NeuralFlow ±0.042(8)7.05\!±\!0.042^(8) ±0.044(9)7.19\!±\!0.044^(9) ±0.283(10)10.13\!±\!0.283^(10) ±0.320(10)10.44\!±\!0.320^(10) ±0.218(10)7.83\!±\!0.218^(10) ±0.212(10)7.86\!±\!0.212^(10) ±0.720(5)6.07\!±\!0.720^(5) ±0.764(6)6.74\!±\!0.764^(6) SeFT ±0.060(9)7.10\!±\!0.060^(9) ±0.060(8)7.15\!±\!0.060^(8) ±0.454(8)9.13\!±\!0.454^(8) ±0.463(8)9.44\!±\!0.463^(8) ±1.706(11)8.94\!±\!1.706^(11) ±1.700(11)8.95\!±\!1.700^(11) ±0.137(10)7.64\!±\!0.137^(10) ±0.170(10)8.41\!±\!0.170^(10) Warpformer ±0.041(4)6.42\!±\!0.041^(4) ±0.045(5)6.50\!±\!0.045^(5) ±0.555(5)5.04\!±\!0.555^(5) ±0.598(5)5.33\!±\!0.598^(5) ±0.065(7)6.78\!±\!0.065^(7) ±0.066(7)6.80\!±\!0.066^(7) ±0.168(6)6.13\!±\!0.168^(6) ±0.195(4)6.67\!±\!0.195^(4) mTAN ±0.034(10)7.47\!±\!0.034^(10) ±0.032(10)7.54\!±\!0.032^(10) ±0.438(9)10.09\!±\!0.438^(9) ±0.443(9)10.39\!±\!0.443^(9) ±0.069(9)7.32\!±\!0.069^(9) ±0.067(9)7.34\!±\!0.067^(9) 4.84¯±0.174(2) 4.84\!±\!0.174^(2) 5.23¯±0.232(2) 5.23\!±\!0.232^(2) tPatchGNN ±0.029(5)6.44\!±\!0.029^(5) ±0.021(4)6.48\!±\!0.021^(4) ±0.322(4)4.55\!±\!0.322^(4) ±0.335(4)4.82\!±\!0.335^(4) ±0.059(5)6.68\!±\!0.059^(5) ±0.060(6)6.73\!±\!0.060^(6) ±1.306(7)6.36\!±\!1.306^(7) ±1.384(7)7.05\!±\!1.384^(7) Appendix D Benchmark Details D.1 Synthetic and Semi-Synthetic Data Common setup. The controlled benchmark contains two fully synthetic datasets, Synthetic-Regime and Synthetic-Multiscale, and two semi-synthetic datasets constructed from ETTm1 and Weather. Each dataset contains 512 training, 128 validation, and 128 test trajectories. For each trajectory, 16 historical observations are provided as input. Training and validation use 30 future target points, whereas each test realization contains 128 future query timestamps. For fully synthetic data, the historical and future intervals are [0,0.5][0,0.5] and (0.5,1](0.5,1], respectively. Gaussian noise with standard deviation 0.0020.002 is added only to the historical observations. Synthetic-Regime. Each trajectory is controlled by z1,z2∼i.i.d.(0,1)z_1,z_2 i.i.d. N(0,1). We first define r(t) r(t) =clip(t−0.580.68−0.58,0,1), =clip ( t-0.580.68-0.58,0,1 ), (42) w(t) w(t) =r(t)23−2r(t). =r(t)^2\3-2r(t)\. The smooth and periodic components are b1(t) b_1(t) =0.10z1+0.015sin(2πt), =0.10z_1+0.015 (2π t), (43) b2(t) b_2(t) =0.10z2+0.012cos(2πt), =0.10z_2+0.012 (2π t), b3(t) b_3(t) =0.07(z1−z2)+0.010sin(4πt), =0.07(z_1-z_2)+0.010 (4π t), p1(t) p_1(t) =(0.48+0.04z2)sin(12π(t−0.64)), =(0.48+0.04z_2) \! (12π(t-0.64) ), p2(t) p_2(t) =(0.44+0.04z1)cos(10π(t−0.64)), =(0.44+0.04z_1) \! (10π(t-0.64) ), p3(t) p_3(t) =0.40+0.03(z1+z2)sin(8π(t−0.64)). =\0.40+0.03(z_1+z_2)\ \! (8π(t-0.64) ). The final three-dimensional trajectory is yj(t)=bj(t)+w(t)pj(t),j∈1,2,3.y_j(t)=b_j(t)+w(t)p_j(t), j∈\1,2,3\. (44) The transition function w(t)w(t) continuously changes the process from a low-amplitude smooth regime to a high-frequency periodic regime, producing different levels of forecasting difficulty across the future interval. Synthetic-Multiscale. For z1,z2,z3∼i.i.d.(0,1)z_1,z_2,z_3 i.i.d. N(0,1), we define g0(t) g_0(t) =0.25z1+0.18t+0.32sin(2π⋅1.25t+0.20z2), =0.25z_1+0.18t+0.32 (2π\!·\!1.25t+0.20z_2), (45) gh(t) g_h(t) =0.14(1+0.08z1)sin(2π⋅8t+0.25z3), =0.14(1+0.08z_1) (2π\!·\!8t+0.25z_3), μ μ =0.72+0.025tanh(z3), =0.72+0.025 (z_3), gb(t) g_b(t) =0.30(1+0.06z2)exp−(t−μ)22(0.055)2. =0.30(1+0.06z_2) \! \- (t-μ)^22(0.055)^2 \. The output channels are y1(t) y_1(t) =g0(t)+gh(t)+gb(t), =g_0(t)+g_h(t)+g_b(t), (46) y2(t) y_2(t) =0.75g0(t)−0.80gh(t)+0.65gb(t) =0.75g_0(t)-0.80g_h(t)+0.65g_b(t) +0.08cos(5πt), +0.08 (5π t), y3(t) y_3(t) =−0.45g0(t)+0.55gh(t)−0.70gb(t) =-0.45g_0(t)+0.55g_h(t)-0.70g_b(t) +0.12t+0.10sin(1.5πt+0.30z1). +0.12t+0.10 (1.5π t+0.30z_1). This construction combines trend, low- and high-frequency periodicity, and a localized transient component. Semi-synthetic datasets. For ETTm1 and Weather, the original regularly sampled series are chronologically divided into training, validation, and test segments with a ratio of 60%/20%/20%60\%/20\%/20\%. Standardization statistics are fitted using the training segment. Each sample contains 96 historical grid points and 512 future grid points. We retain 16 historical points as input and sample 30 future points for training and validation. At test time, 128 future points are selected from the complete 512-point future grid, which is retained to compute the reference continuous-time risk. For all four datasets, test timestamps are sampled from pα(u)=αBeta(u,10,60)+(1−α)Uniform(0,1),p_α(u)=α\,Beta(u;10,60)+(1-α)\,Uniform(0,1), (47) where α∈0,0.3,0.5,0.9,0.99α∈\0,0.3,0.5,0.9,0.99\ controls the degree of temporal non-uniformity. For each α, the trained model, historical observations, and underlying future trajectory are fixed, while only the test query timestamps are resampled. Each setting is evaluated using 20 timestamp resamples shared across all models. Table 5: Complete MAE and CAE results on CESNET, FNSPID, GDELT, and HumanActivity. Values are mean ± standard deviation over three random seeds; standard deviations are rounded to three decimal places. Model CESNET ×10−1× 10^-1 FNSPID ×10−1× 10^-1 GDELT ×10−1× 10^-1 HumanActivity ×10−1× 10^-1 MAE CAE MAE CAE MAE CAE MAE CAE APN ±0.003(5)6.69\!±\!0.003^(5) ±0.003(5)6.62\!±\!0.003^(5) 1.86¯±0.173(2) 1.86\!±\!0.173^(2) 1.85¯±0.173(2) 1.85\!±\!0.173^(2) ±0.066(1)6.59\!±\!0.066^(1) ±0.067(1)6.50\!±\!0.067^(1) ±0.031(6)1.43\!±\!0.031^(6) ±0.032(6)1.42\!±\!0.032^(6) ASTGI ±0.064(1)6.55\!±\!0.064^(1) ±0.067(1)6.48\!±\!0.067^(1) ±0.092(6)2.19\!±\!0.092^(6) ±0.093(6)2.18\!±\!0.093^(6) ±0.088(4)6.69\!±\!0.088^(4) ±0.088(4)6.60\!±\!0.088^(4) ±0.024(3)1.38\!±\!0.024^(3) ±0.024(3)1.38\!±\!0.024^(3) GRU-D ±0.249(10)7.45\!±\!0.249^(10) ±0.262(10)7.37\!±\!0.262^(10) ±0.286(8)2.76\!±\!0.286^(8) ±0.284(8)2.75\!±\!0.284^(8) ±0.033(10)6.81\!±\!0.033^(10) ±0.032(10)6.72\!±\!0.032^(10) ±0.251(9)2.06\!±\!0.251^(9) ±0.253(9)2.06\!±\!0.253^(9) GraFITi 6.56¯±0.033(2) 6.56\!±\!0.033^(2) 6.49¯±0.033(2) 6.49\!±\!0.033^(2) ±0.097(4)1.96\!±\!0.097^(4) ±0.096(4)1.95\!±\!0.096^(4) ±0.146(3)6.66\!±\!0.146^(3) ±0.146(3)6.57\!±\!0.146^(3) 1.37¯±0.012(2) 1.37\!±\!0.012^(2) 1.37¯±0.013(2) 1.37\!±\!0.013^(2) HyperIMTS ±0.056(7)6.79\!±\!0.056^(7) ±0.063(7)6.75\!±\!0.063^(7) ±0.391(3)1.91\!±\!0.391^(3) ±0.390(3)1.91\!±\!0.390^(3) ±0.072(8)6.74\!±\!0.072^(8) ±0.074(8)6.65\!±\!0.074^(8) ±0.013(1)1.36\!±\!0.013^(1) ±0.013(1)1.36\!±\!0.013^(1) KAFNet ±0.027(4)6.68\!±\!0.027^(4) ±0.031(4)6.61\!±\!0.031^(4) ±0.055(1)1.73\!±\!0.055^(1) ±0.056(1)1.73\!±\!0.056^(1) 6.60¯±0.092(2) 6.60\!±\!0.092^(2) 6.51¯±0.093(2) 6.51\!±\!0.093^(2) ±0.010(4)1.39\!±\!0.010^(4) ±0.010(4)1.39\!±\!0.010^(4) NeuralFlow ±0.122(8)7.41\!±\!0.122^(8) ±0.127(8)7.35\!±\!0.127^(8) ±0.234(10)3.46\!±\!0.234^(10) ±0.233(10)3.45\!±\!0.233^(10) ±0.064(9)6.77\!±\!0.064^(9) ±0.064(9)6.68\!±\!0.064^(9) ±0.305(11)3.11\!±\!0.305^(11) ±0.304(11)3.11\!±\!0.304^(11) SeFT ±0.118(9)7.44\!±\!0.118^(9) ±0.115(9)7.36\!±\!0.115^(9) ±0.828(11)4.33\!±\!0.828^(11) ±0.827(11)4.32\!±\!0.827^(11) ±0.062(6)6.72\!±\!0.062^(6) ±0.059(5)6.61\!±\!0.059^(5) ±0.041(10)2.12\!±\!0.041^(10) ±0.041(10)2.12\!±\!0.041^(10) Warpformer ±0.054(3)6.63\!±\!0.054^(3) ±0.059(3)6.56\!±\!0.059^(3) ±0.053(7)2.20\!±\!0.053^(7) ±0.054(7)2.19\!±\!0.054^(7) ±0.110(11)7.48\!±\!0.110^(11) ±0.109(11)7.35\!±\!0.109^(11) ±0.003(5)1.41\!±\!0.003^(5) ±0.003(5)1.41\!±\!0.003^(5) mTAN ±0.010(11)7.61\!±\!0.010^(11) ±0.016(11)7.53\!±\!0.016^(11) ±0.314(9)2.96\!±\!0.314^(9) ±0.314(9)2.94\!±\!0.314^(9) ±0.062(5)6.70\!±\!0.062^(5) ±0.063(6)6.63\!±\!0.063^(6) ±0.018(8)1.85\!±\!0.018^(8) ±0.018(8)1.85\!±\!0.018^(8) tPatchGNN ±0.070(6)6.75\!±\!0.070^(6) ±0.082(6)6.70\!±\!0.082^(6) ±0.247(5)2.16\!±\!0.247^(5) ±0.246(5)2.16\!±\!0.246^(5) ±0.083(7)6.73\!±\!0.083^(7) ±0.083(7)6.64\!±\!0.083^(7) ±0.035(7)1.45\!±\!0.035^(7) ±0.036(7)1.44\!±\!0.036^(7) For fully synthetic data, the reference RctR_ct is computed by trapezoidal integration over a dense common time grid. For semi-synthetic data, it is computed as the equally weighted squared error over all 512 future grid points. D.2 Real-World Datasets We evaluate models on eight real-world datasets spanning healthcare, climate science, human sensing, finance, software engineering, international events, and network systems. These datasets exhibit diverse irregularity mechanisms, including event-triggered logging, activity-dependent collection, operational constraints, human scheduling, missing observations, and system jitter (4). Only the numerical time-series modalities are used. The processed data are represented either as wide tables with one row per record and timestamp or as sparse triplets of timestamp, variable index, and value. Sparse observations are converted into multivariate tensors with binary masks. Filled values at unobserved positions are used only for tensor storage and do not contribute to training or evaluation. Record identifiers are deterministically divided into training, validation, and test partitions, and forecasting windows are constructed without crossing record boundaries. CESNET. CESNET contains network traffic measurements collected from 11 devices, with 10 traffic variables. Its irregular timestamps mainly arise from scheduling delay and logging jitter. Raw timestamps are aggregated to minute resolution, and duplicate observations within the same minute are averaged. We use 1,440 minutes of history to predict the following 1,440 minutes, producing 252, 84, and 126 training, validation, and test windows, respectively. FNSPID. FNSPID contains six financial variables for 10 entities, including open, high, low, close, adjusted close, and trading volume. Observations are restricted by market operating hours, creating structured gaps between trading periods. Timestamps are converted into relative calendar days without filling non-trading days. Both the historical and forecasting windows contain 30 days, resulting in 496/171/176 train/validation/test samples. GDELT. GDELT records international events and contains five numerical variables for eight event streams. Since observations are generated when geopolitical events occur, the sampling process is naturally event-driven. Individual events and fractional-day timestamps are retained without daily aggregation. We use 14-day historical and forecasting windows, yielding 208/52/156 train/validation/test samples. HumanActivity. HumanActivity contains 12 location variables obtained from sensors placed on the body during different human activities. The dataset contains 25 records, and the sensors are observed at asynchronous and slightly different timestamps. Relative timestamps are quantized, while duplicate measurements from the same sensor and timestamp are averaged. The history and forecasting lengths are 3,000 and 1,000 time units, producing 735/264/361 train/validation/test samples. MIMIC-I. MIMIC-I is a large critical-care database containing de-identified clinical records from ICU patients (9). We combine laboratory, input, output, and prescription events from the first 48 hours of each admission and aggregate them into half-hour bins. The processed data contain 22,406 admissions and 96 variables. The first 72 bins are used as history and the following 24 bins as the forecasting horizon, resulting in 11,156/3,699/3,703 train/validation/test samples. Its irregularity mainly reflects clinician-driven measurements and heterogeneous clinical recording frequencies. Table 6: Complete MAE and CAE results on MIMIC-I, RepoHealth, StudentLife, and USHCN. Values are mean ± standard deviation over three random seeds; standard deviations are rounded to three decimal places. Model MIMIC-I ×10−1× 10^-1 RepoHealth ×10−1× 10^-1 StudentLife ×10−1× 10^-1 USHCN ×10−1× 10^-1 MAE CAE MAE CAE MAE CAE MAE CAE APN ±0.034(6)4.78\!±\!0.034^(6) ±0.034(6)4.79\!±\!0.034^(6) 3.47¯±0.049(2) 3.47\!±\!0.049^(2) 3.55¯±0.051(2) 3.55\!±\!0.051^(2) 5.75¯±0.017(2) 5.75\!±\!0.017^(2) ±0.017(3)5.78\!±\!0.017^(3) ±0.235(4)3.56\!±\!0.235^(4) ±0.291(5)3.70\!±\!0.291^(5) ASTGI ±0.001(11)6.74\!±\!0.001^(11) ±0.001(11)6.76\!±\!0.001^(11) ±1.554(7)5.27\!±\!1.554^(7) ±1.584(7)5.36\!±\!1.584^(7) ±0.075(3)5.77\!±\!0.075^(3) 5.77¯±0.075(2) 5.77\!±\!0.075^(2) ±0.080(3)3.11\!±\!0.080^(3) ±0.085(3)3.21\!±\!0.085^(3) GRU-D ±0.018(7)4.98\!±\!0.018^(7) ±0.020(7)5.00\!±\!0.020^(7) ±0.614(11)7.56\!±\!0.614^(11) ±0.600(11)7.65\!±\!0.600^(11) ±0.023(8)6.15\!±\!0.023^(8) ±0.024(8)6.16\!±\!0.024^(8) ±0.408(11)4.53\!±\!0.408^(11) ±0.408(11)4.69\!±\!0.408^(11) GraFITi ±0.072(1)4.40\!±\!0.072^(1) 4.47¯±0.073(2) 4.47\!±\!0.073^(2) ±0.588(1)3.21\!±\!0.588^(1) ±0.587(1)3.26\!±\!0.587^(1) ±0.035(6)5.83\!±\!0.035^(6) ±0.034(4)5.79\!±\!0.034^(4) ±0.177(1)2.97\!±\!0.177^(1) ±0.191(1)3.08\!±\!0.191^(1) HyperIMTS 4.46¯±0.035(2) 4.46\!±\!0.035^(2) ±0.035(1)4.41\!±\!0.035^(1) ±0.370(3)3.89\!±\!0.370^(3) ±0.374(3)3.97\!±\!0.374^(3) ±0.024(1)5.74\!±\!0.024^(1) ±0.024(1)5.76\!±\!0.024^(1) ±0.967(8)4.15\!±\!0.967^(8) ±0.996(8)4.29\!±\!0.996^(8) KAFNet ±0.035(3)4.60\!±\!0.035^(3) ±0.033(3)4.61\!±\!0.033^(3) ±1.978(6)5.17\!±\!1.978^(6) ±1.997(6)5.25\!±\!1.997^(6) ±0.023(4)5.78\!±\!0.023^(4) ±0.023(5)5.83\!±\!0.023^(5) ±0.291(9)4.31\!±\!0.291^(9) ±0.288(9)4.47\!±\!0.288^(9) NeuralFlow ±0.015(8)5.11\!±\!0.015^(8) ±0.015(8)5.12\!±\!0.015^(8) ±0.510(10)7.36\!±\!0.510^(10) ±0.531(10)7.46\!±\!0.531^(10) ±0.075(10)6.49\!±\!0.075^(10) ±0.073(10)6.51\!±\!0.073^(10) ±0.615(5)3.58\!±\!0.615^(5) ±0.638(6)3.77\!±\!0.638^(6) SeFT ±0.048(9)5.15\!±\!0.048^(9) ±0.046(9)5.16\!±\!0.046^(9) ±0.152(8)6.63\!±\!0.152^(8) ±0.152(8)6.72\!±\!0.152^(8) ±0.851(11)6.93\!±\!0.851^(11) ±0.849(11)6.94\!±\!0.849^(11) ±0.008(10)4.43\!±\!0.008^(10) ±0.004(10)4.59\!±\!0.004^(10) Warpformer ±0.025(4)4.63\!±\!0.025^(4) ±0.023(4)4.64\!±\!0.023^(4) ±0.257(5)4.86\!±\!0.257^(5) ±0.260(5)4.98\!±\!0.260^(5) ±0.052(7)5.85\!±\!0.052^(7) ±0.052(7)5.86\!±\!0.052^(7) ±0.026(6)3.66\!±\!0.026^(6) ±0.022(4)3.67\!±\!0.022^(4) mTAN ±0.043(10)5.38\!±\!0.043^(10) ±0.045(10)5.39\!±\!0.045^(10) ±0.242(9)6.96\!±\!0.242^(9) ±0.229(9)7.05\!±\!0.229^(9) ±0.127(9)6.19\!±\!0.127^(9) ±0.125(9)6.20\!±\!0.125^(9) 3.01¯±0.106(2) 3.01\!±\!0.106^(2) 3.13¯±0.107(2) 3.13\!±\!0.107^(2) tPatchGNN ±0.023(5)4.65\!±\!0.023^(5) ±0.021(5)4.65\!±\!0.021^(5) ±0.322(4)4.36\!±\!0.322^(4) ±0.324(4)4.43\!±\!0.324^(4) ±0.081(5)5.82\!±\!0.081^(5) ±0.081(6)5.84\!±\!0.081^(6) ±0.725(7)3.82\!±\!0.725^(7) ±0.738(7)3.94\!±\!0.738^(7) RepoHealth. RepoHealth describes issue, pull-request, and development activities from four software repositories using 10 numerical variables. Observations become denser during active development periods and sparser during inactive periods. Timestamps are converted into relative days, and duplicated daily observations are averaged. We use 30-day history and forecasting windows, with 97/59/50 train/validation/test samples. StudentLife. StudentLife contains daily mobile-sensing measurements from 20 students, including activity, location, sleep, and phone-use variables. Its nine variables are irregularly observed because data availability depends on individual routines and device usage. Both the historical and forecasting windows contain 30 days, producing 271/89/90 train/validation/test samples. USHCN. USHCN contains long-term climate observations from 1,114 weather stations across the United States (16). Five meteorological variables are used, and missing daily observations make the resulting trajectories sporadic. Following the commonly used irregular forecasting protocol, the first 150 time units are used as history and the following 50 units as the forecasting horizon. The resulting split contains 666/220/224 train/validation/test samples. D.3 Baselines and Hyperparameters We compare eleven representative models covering recurrent, continuous-time, set-based, attention-based, graph-based, and pre-alignment forecasting paradigms. To ensure consistent comparisons, each model uses the same configuration across all real-world datasets. APN. APN aggregates irregular observations into learnable temporal patches through soft boundaries and predicts values at arbitrary query timestamps (15). We use hidden dimension d=56d=56, two layers, eight attention heads, time-embedding dimension 8, two patches, dropout 0.1, and batch size 32. ASTGI. ASTGI directly represents each discrete observation as a point in a learnable spatio-temporal space and adaptively constructs causal nearest-neighbor graphs for information propagation (14). We use hidden dimension 128, two propagation layers, time-embedding dimension 128, channel-embedding dimension 64, 96 candidate neighbors, MLP ratio 4, dropout 0.1, and batch size 8. GRU-D. GRU-D extends recurrent models with trainable decay mechanisms that explicitly incorporate missingness masks and elapsed time between observations (5). We use one recurrent layer with hidden dimension 128 and batch size 32. GraFITi. GraFITi represents observed and queried values as nodes in a bipartite graph and performs graph-based information propagation for irregular forecasting (23). We use hidden dimension 128, four graph layers, one attention head, zero dropout, and batch size 32. HyperIMTS. HyperIMTS constructs hypergraph interactions to model higher-order dependencies among asynchronously observed variables (13). We use hidden dimension 128, three layers, one attention head, dropout 0.1, and batch size 32. KAFNet. KAFNet employs kernel aggregation to pre-align irregular observations before temporal modeling (27). We use hidden dimension 16, pre-convolution dimension 32, one layer, one attention head, four Gaussian kernels, time-embedding dimension 10, and batch size 32. NeuralFlow. NeuralFlow parameterizes continuous latent dynamics through an invertible flow, providing an efficient alternative to numerical ODE solvers (2). We use hidden dimension 100, two coupling-flow layers, latent dimension 20, time-encoding hidden dimension 8, three hidden layers, reconstruction dimension 30, and batch size 32. SeFT. SeFT treats all irregular observations as an unordered set and uses learned set functions and attention pooling to obtain a fixed-dimensional representation (8). We use hidden dimension 128, three ϕ/ρφ/ρ layers, two ψ layers, four attention heads, positional dimension 4, dropout 0.1, and batch size 32. (a) Synthetic-Regime (b) Synthetic-Multiscale (c) ETTm1 (d) Weather Figure 5: Relative estimation errors of MSE and CSE with respect to RctR_ct under increasing sampling non-uniformity on all synthetic and semi-synthetic datasets. Warpformer. Warpformer applies multiscale temporal warping to transform irregular observations into representations suitable for Transformer-based modeling (24). We use hidden dimension 128, three layers, inner dimension 64, one attention head, key/value dimensions of 8, and batch size 32. mTAN. mTAN employs multi-time attention to map irregular observations onto continuous reference timestamps (20). We use hidden dimension 128, latent dimension 32, one attention head, 128 reference points, observation standard deviation 0.01, and batch size 32. tPatchGNN. tPatchGNN organizes irregular observations into transformable temporal patches and models variable interactions using graph neural networks (25). We use hidden dimension 32, one graph layer, one attention head, one Transformer layer, one-hop propagation, node-embedding dimension 10, time-embedding dimension 10, a linear output head, and batch size 32. Controlled-experiment baselines. The controlled experiments use SeFT, CRU, GraFITi, tPatchGNN, and KAFNet. SeFT uses hidden dimension 64, two ϕ/ρφ/ρ layers, two ψ layers, two attention heads, and dropout 0.02. CRU uses hidden dimension 20 and 20 basis functions (19). GraFITi uses hidden dimension 28 and two graph layers. tPatchGNN uses hidden dimension 32, patch length 4, one graph layer, and one Transformer layer. KAFNet uses hidden dimension 64, pre-convolution dimension 32, two layers, eight Gaussian kernels, and time-embedding dimension 16. Table 7: Complete model-evaluation results on the two synthetic datasets under α=0.9α=0.9. Values are mean ± standard deviation over 20 timestamp resamples; RctR_ct is computed from the complete trajectory. Superscripts denote ranks. Model Synthetic-Regime ×10−2× 10^-2 Synthetic-Multiscale ×10−2× 10^-2 MSE CSE RctR_ct MSE CSE RctR_ct SeFT ±0.0228(3)1.0884\!±\!0.0228^(3) ±0.3760(1)8.6994\!±\!0.3760^(1) 6.2552(1)6.2552^(1) ±0.0058(3)1.2655\!±\!0.0058^(3) ±0.0632(3)1.1169\!±\!0.0632^(3) 1.1659(3)1.1659^(3) CRU ±0.0232(1)0.6527\!±\!0.0232^(1) 8.8105¯±0.3814(2) 8.8105\!±\!0.3814^(2) 6.4831¯(2) 6.4831^(2) ±0.0055(4)1.3118\!±\!0.0055^(4) ±0.0380(1)0.8085\!±\!0.0380^(1) 0.9575(1)0.9575^(1) GraFITi ±0.0233(5)2.5741\!±\!0.0233^(5) ±0.4644(5)9.8776\!±\!0.4644^(5) 7.6588(5)7.6588^(5) ±0.0053(5)1.6822\!±\!0.0053^(5) ±0.0477(5)1.4708\!±\!0.0477^(5) 1.5059(5)1.5059^(5) tPatchGNN 0.9132¯±0.0241(2) 0.9132\!±\!0.0241^(2) ±0.3603(4)9.3798\!±\!0.3603^(4) 6.9119(4)6.9119^(4) 0.8939¯±0.0052(2) 0.8939\!±\!0.0052^(2) ±0.0687(4)1.1272\!±\!0.0687^(4) 1.1789(4)1.1789^(4) KAFNet ±0.0232(4)1.5435\!±\!0.0232^(4) ±0.3602(3)9.0437\!±\!0.3602^(3) 6.7512(3)6.7512^(3) ±0.0057(1)0.8088\!±\!0.0057^(1) 1.0751¯±0.0511(2) 1.0751\!±\!0.0511^(2) 1.1567¯(2) 1.1567^(2) D.4 Training and Evaluation Details Training protocol. All models are trained using the Adam optimizer and masked observation-point MSE. For prediction mask mbidm_bid, the training objective is ℒtrain=∑b,i,dmbid(y^bid−ybid)2∑b,i,dmbid.L_train= _b,i,dm_bid ( y_bid-y_bid )^2 _b,i,dm_bid. (48) All baselines are trained for 100 epochs, and the checkpoint with the lowest validation MSE is retained for testing. The learning rate is 10−310^-3 for all models. NeuralFlow and mTAN use weight decay 10−410^-4, while the remaining real-world baselines use no weight decay. Models that originally employ probabilistic objectives, including mTAN and NeuralFlow, are also optimized using the common MSE objective to ensure a unified training protocol. The controlled experiments use learning rate 10−310^-3, weight decay 10−510^-5, batch size 64, and global gradient-norm clipping at 1.0. Training and validation target timestamps are uniformly sampled from the future interval. All model–dataset experiments are independently repeated with three random seeds, 2024, 2025, and 2026. We report the mean and standard deviation across the three runs. In the controlled sampling-shift experiment, each trained checkpoint is further evaluated using 20 independent timestamp resamples for each sampling strength α. Temporal-density estimation. Future timestamps are mapped to the common forecasting interval [0,1][0,1] before density estimation. We use a trajectory-specific leave-one-out Gaussian kernel density estimator. The bandwidth follows a robust Silverman rule: hs=max0.9min(σs,IQRs1.34)Ls−1/5,1max(20,2Ls),10−6.h_s= \0.9 ( _s, IQR_s1.34 )L_s^-1/5, 1 (20,2L_s),10^-6 \. (49) Non-positive scale estimates are excluded when computing the minimum. Estimated densities are lower-bounded by 10−810^-8, and no importance-weight clipping or boundary reflection is applied. Conditional and marginal timestamp densities are estimated only from test timestamps and are shared by all models evaluated on the same dataset. Metric computation. For each model–dataset–seed combination, the selected checkpoint is loaded once, and MSE, GSE, and CSE are computed from the same test predictions. Let msidm_sid indicate whether variable d is observed at future timestamp ts,it_s,i and let ws,iw_s,i denote the corresponding temporal weight. The trajectory-level weighted error is R^s(w)=∑i,dmsidws,i(y^sid−ysid)2∑i,dmsidws,i. R_s(w)= _i,dm_sidw_s,i ( y_sid-y_sid )^2 _i,dm_sidw_s,i. (50) The observation-distribution, global-time, and continuous-time weights are ws,iobs w^obs_s,i =1, =1, (51) ws,iglobal w^global_s,i =p^T(ts,i)p^(ts,i∣Xs), = p_T(t_s,i) p(t_s,i X_s), ws,ict w^ct_s,i =1p^(ts,i∣Xs). = 1 p(t_s,i X_s). Each trajectory is normalized independently, after which the trajectory-level errors are averaged with equal weight. Padding positions and timestamps without any observed target variables are excluded from both density estimation and metric computation. Implementation environment. All experiments are implemented in Python 3.12 and PyTorch 2.10.0 with CUDA 12.8. Training and evaluation are conducted on NVIDIA GeForce RTX 4090 GPUs. Appendix E Additional Experimental Results E.1 Complete Results on Synthetic and Semi-Synthetic Data Figure 5 reports the relative estimation errors of MSE and CSE with respect to RctR_ct on all four controlled datasets. The curves show the mean over 20 independent test-timestamp resamples, and the shaded regions denote one standard deviation. Because the trained models and underlying future trajectories remain fixed across different values of α, the observed changes are caused solely by the shift in the test-time sampling distribution. Table 8: Complete model-evaluation results on the two semi-synthetic datasets under α=0.9α=0.9. Values are mean ± standard deviation over 20 timestamp resamples; RctR_ct is computed on the complete future grid. Superscripts denote ranks. Model ETTm1 ×100× 10^0 Weather ×10−1× 10^-1 MSE CSE RctR_ct MSE CSE RctR_ct SeFT ±0.0059(4)1.1032\!±\!0.0059^(4) ±0.0153(4)1.1145\!±\!0.0153^(4) 1.1130(4)1.1130^(4) ±0.0146(5)4.4046\!±\!0.0146^(5) ±0.0629(5)4.9629\!±\!0.0629^(5) 4.8493(5)4.8493^(5) CRU ±0.0058(5)1.3394\!±\!0.0058^(5) ±0.0153(3)1.1001\!±\!0.0153^(3) 1.0994(3)1.0994^(3) ±0.0144(1)3.0839\!±\!0.0144^(1) ±0.0667(3)4.6087\!±\!0.0667^(3) 4.4631(3)4.4631^(3) GraFITi ±0.0056(3)1.0886\!±\!0.0056^(3) ±0.0163(5)1.4140\!±\!0.0163^(5) 1.4016(5)1.4016^(5) ±0.0147(4)3.8918\!±\!0.0147^(4) ±0.0665(4)4.9221\!±\!0.0665^(4) 4.7612(4)4.7612^(4) tPatchGNN 0.9251¯±0.0056(2) 0.9251\!±\!0.0056^(2) ±0.0151(1)0.9330\!±\!0.0151^(1) 0.9302(1)0.9302^(1) ±0.0181(3)3.8789\!±\!0.0181^(3) ±0.0645(1)3.6856\!±\!0.0645^(1) 3.5274(1)3.5274^(1) KAFNet ±0.0056(1)0.9093\!±\!0.0056^(1) 0.9640¯±0.0155(2) 0.9640\!±\!0.0155^(2) 0.9584¯(2) 0.9584^(2) 3.4839¯±0.0144(2) 3.4839\!±\!0.0144^(2) 3.7783¯±0.0607(2) 3.7783\!±\!0.0607^(2) 3.6496¯(2) 3.6496^(2) Under uniform sampling, the inverse-density weights are constant and CSE consequently reduces to MSE. As α increases, observations become increasingly concentrated in a limited portion of the future interval. MSE therefore places disproportionate emphasis on densely sampled regions, and its deviation from RctR_ct generally increases. In contrast, CSE compensates for this concentration through inverse-density weighting and remains substantially closer to RctR_ct under strong non-uniform sampling. Although density estimation and importance weighting may introduce additional finite-sample variability, the correction becomes increasingly beneficial as the mismatch between the observed and uniform temporal distributions grows. The consistent behavior on ETTm1 and Weather further shows that this advantage is not restricted to analytically generated trajectories. Tables 7 and 8 report the complete model-level results under α=0.9α=0.9. On all four datasets, CSE recovers exactly the same ranking of the five models as RctR_ct, whereas MSE produces at least one ranking inversion. This result indicates that the effect of temporal sampling is not limited to a uniform shift in metric values: different models are affected to different degrees because their prediction errors vary differently over time. Correcting the target-time distribution therefore improves both risk estimation and the recovery of model comparisons, supporting the main-text conclusion that CSE more accurately reflects continuous-time predictive performance under non-uniform sampling. E.2 Complete Results on Real-World Datasets Tables 3 and 4 report the complete MSE and CSE results on the eight real-world datasets, including standard deviations over three random seeds. Tables 5 and 6 provide the corresponding MAE and Continuous-Time Absolute Error (CAE) results. The best and second-best mean values are shown in bold and underlined, respectively, and superscripts denote ranks computed from the unrounded means. Consistent with the observations in the main text, CSE is neither a fixed-direction nor a fixed-ratio transformation of MSE. On CESNET, FNSPID, GDELT, and HumanActivity, CSE is generally lower than MSE, whereas the opposite trend is more common on MIMIC-I, RepoHealth, and USHCN. Moreover, even within the same dataset, the magnitude of the change varies across models. This confirms that the sampling distribution interacts with model-specific temporal error patterns rather than simply applying a dataset-level scaling factor. The reported standard deviations also distinguish this systematic metric effect from training variability: most model–dataset combinations remain relatively stable across seeds, although datasets such as RepoHealth and USHCN exhibit larger variation for several models. The MAE–CAE comparison exhibits a similar pattern. The direction of the correction is broadly consistent with the MSE–CSE comparison, while its magnitude remains dependent on both the dataset and the model. In particular, the best-performing model on MIMIC-I changes from GraFITi under MAE to HyperIMTS under CAE, accompanied by additional local rank changes on several datasets. Therefore, the discrepancies reported in the main text are not specific to squared error, but arise more generally from the temporal distribution under which pointwise prediction losses are aggregated. Since the complete continuous trajectories of real-world datasets are unavailable, these results do not directly identify which metric is exact; rather, the systematic differences suggest that relying solely on observation-point metrics may not fully characterize models’ continuous-time predictive performance. langley00