Paper deep dive
Discretizing Continuous Time Series for Imputation with Masked Diffusion Training
Dongbin Kim, Seungyun Lee, Geonwoo Shin, Jaewook Lee
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/22/2026, 1:38:53 AM
Summary
The paper introduces the Masked Diffusion Time-series Imputation Model (MDTIM), a framework that adapts masked diffusion models for continuous time series imputation. It addresses limitations in existing methods by using a structurally orthogonal [MASK] token and Stochastic Discretization to map continuous values to ordinal-aware tokens. The model employs Ordinal-Aware Soft Labeling and Spectral Consistency Regularization to preserve temporal dynamics and global coherence, outperforming state-of-the-art baselines in robustness and accuracy.
Entities (8)
Relation Signals (6)
MDTIM → isbasedon → Masked Diffusion Model
confidence 95% · we propose the Masked Diffusion Time-series Imputation Model (MDTIM), which leverages the training paradigm of masked diffusion model
MDTIM → uses → Stochastic Discretization
confidence 95% · we further introduce Stochastic Discretization, which maps continuous values to ordinal-aware tokens while preserving continuous dynamics.
MDTIM → uses → Ordinal-Aware Soft Labeling
confidence 95% · we employ Ordinal-Aware Soft Labeling with a truncated Gaussian kernel.
Stochastic Discretization → bridgesgapbetween → discrete masked diffusion
confidence 90% · To bridge the gap between discrete masked diffusion and the continuous, ordinal nature of time series, we further introduce Stochastic Discretization
MDTIM → outperforms → CSDI
confidence 90% · Compared to the diffusion-based CSDI [33], MDTIM consistently yields lower MAE across all datasets
MDTIM → outperforms → SAITS
confidence 85% · MDTIM achieves the lowest MAE on Energy, ETTh, and Weather... MDTIM is comparable to SAITS at L=48, the gap widens substant
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Time series imputation is a crucial area for reliable time series analysis, yet it remains challenging due to the complex temporal dynamics and noise of real-world data. Existing approaches, however, exhibit two limitations: missing and observed values are embedded within the same representation space without explicit structural separation, and continuous diffusion-based methods are trained to predict added noise rather than the original signal. To address these, we propose the Masked Diffusion Time-series Imputation Model (MDTIM), which leverages the training paradigm of masked diffusion model for imputation tasks. The MASK token is structurally orthogonal to valid observations, and the model directly predicts the original values, naturally aligning both the representation and the learning objective with the imputation task. To bridge the gap between discrete masked diffusion and the continuous, ordinal nature of time series, we further introduce Stochastic Discretization, which maps continuous values to ordinal-aware tokens while preserving continuous dynamics. Our experiments on diverse benchmarks confirm that MDTIM achieves superior robustness and scalability, consistently outperforming state-of-the-art deterministic and generative baselines across various missing scenarios.
Tags
Links
- Source: https://arxiv.org/abs/2608.19119v1
- Canonical: https://arxiv.org/abs/2608.19119v1
Trouble viewing inline? Open PDF directly →
Full Text
73,688 characters extracted from source content.
Expand or collapse full text
Discretizing Continuous Time Series for Imputation with Masked Diffusion Training Dongbin Kim Affiliation: [1ex] Seoul National University[1ex] dongbin413,rats96,shin0621,jaewook@snu.ac.kr Seungyun Lee Affiliation: [1ex] Seoul National University[1ex] dongbin413,rats96,shin0621,jaewook@snu.ac.kr Geonwoo Shin Affiliation: [1ex] Seoul National University[1ex] dongbin413,rats96,shin0621,jaewook@snu.ac.kr Jaewook Lee Abstract Time series imputation is a crucial area for reliable time series analysis, yet it remains challenging due to the complex temporal dynamics and noise of real-world data. Existing approaches, however, exhibit two limitations: missing and observed values are embedded within the same representation space without explicit structural separation, and continuous diffusion-based methods are trained to predict added noise rather than the original signal. To address these, we propose the Masked Diffusion Time-series Imputation Model (MDTIM), which leverages the training paradigm of masked diffusion model for imputation tasks. The [MASK] token is structurally orthogonal to valid observations, and the model directly predicts the original values, naturally aligning both the representation and the learning objective with the imputation task. To bridge the gap between discrete masked diffusion and the continuous, ordinal nature of time series, we further introduce Stochastic Discretization, which maps continuous values to ordinal-aware tokens while preserving continuous dynamics. Our experiments on diverse benchmarks confirm that MDTIM achieves superior robustness and scalability, consistently outperforming state-of-the-art deterministic and generative baselines across various missing scenarios. 11footnotetext: Corresponding author. 1 Introduction Time series imputation is a fundamental task in time series analysis, as real-world data are often partially observed due to sensor failures or transmission errors [16]. At its core, imputation aims to recover the underlying clean signal from partial observations, which requires the representation and learning objective of the model to be well-aligned with the structure of partial observability. From this perspective, we observe that existing approaches leave room for improvement along two axes. First, at the representation level, both discriminative Transformer-based methods [6] and continuous diffusion models [33] embed missing and observed values within the same data manifold. Although attention masks can prevent the model from attending to missing positions, placeholder values (e.g., zeros) still occupy the same embedding space as valid observations, potentially introducing spurious correlations in the input representation [8]. Second, at the objective level, continuous diffusion models such as CSDI [33] are typically trained to predict the Gaussian noise added to the data rather than the original values themselves, which is somewhat indirect for a task whose goal is to recover what was originally there. Masked Diffusion Models (MDM), recently developed in Natural Language Processing [28], offer a natural fit for both aspects. Unlike continuous diffusion that corrupts data with Gaussian noise, MDM replaces tokens with a special [MASK] symbol and learns to recover the original values. The [MASK] token is structurally orthogonal to any valid observation, which provides a clearer separation between missing and observed states in the representation. The model also directly predicts the original values rather than the noise, more closely matching the imputation goal. Despite this conceptual fit, directly applying MDM to time series imputation presents a fundamental challenge. Time series data are inherently continuous and ordinal, whereas language tokens are discrete and categorical. Naive quantization severs the ordinal relationships between adjacent values and loses information below the grid resolution, which is particularly problematic when fine-grained reconstruction is needed. This raises our central research question: How can we leverage the training paradigm of masked diffusion while preserving the continuous and ordinal nature of time series data? Figure 1: Imputation results of MDTIM on the ETTh dataset. Red regions indicate intervals where both channels are missing, while yellow regions represent single-channel missingness. The uncertainty (σ), derived from the output probability of MDTIM, is estimated to be higher in the red regions (information scarcity) compared to the yellow regions (richer information). To answer this, we propose the Masked Diffusion Time-series Imputation Model (MDTIM), a framework that bridges discrete masked diffusion and continuous time series modeling. We introduce Stochastic Discretization, which injects noise during tokenization to preserve information in expectation, and replace the standard one-hot objective with Ordinal-Aware Soft Labeling to capture the ordinal relationships between tokens. We further incorporate a spectral consistency objective and an expectation-based decoder for precise continuous reconstruction. Our main contributions are summarized as follows: • We propose a Masked Diffusion Framework that addresses the lack of structural separation between missing and observed representations in time series imputation. • We introduce a novel quantization pipeline by Stochastic Discretization and Ordinal-Aware Soft Labeling, which preserves ordinal relationships within a categorical token space without information loss. • We employ an Expectation-based Unmasking strategy and Spectral Consistency Regularization, allowing MDTIM to reconstruct continuous time series precisely. Our method consistently outperforms state-of-the-art deterministic and generative baselines across diverse benchmarks. 2 Related Work 2.1 Time Series Imputation Time series data are often partially observed due to sensor faults, irregular sampling, or data-collection constraints [35]. Naive deletion or mean/zero imputation can bias estimation and degrade downstream tasks [16], so time-series imputation aims to recover missing values by leveraging temporal dependencies, cross-variable correlations, and informative missingness patterns [4, 6]. Methodologically, deep learning approaches have progressed from RNN-based architectures with masking [4, 3] to Transformer-based models that leverage self-attention to capture long-range dependencies; SAITS [6], in particular, achieves strong performance via Diagonally-Masked Self-Attention. Beyond deterministic models, diffusion-based models such as CSDI [33] and SSSD [20] produce probabilistic imputations by conditioning on observed data and iterative denoising. SSSD further incorporates structured state space models to better capture long-range temporal dependencies. However, this conditional formulation is misaligned with the imputation task: the model learns “what noise was added” rather than “what was originally there”. Moreover, its training objective is decoupled from the masking ratio, applying uniform weight regardless of reconstruction difficulty. 2.2 Masked Diffusion Model Discrete diffusion models extend the diffusion paradigm to categorical state spaces [13]. Among them, D3PM [2] introduced an absorbing-state formulation, where tokens are progressively replaced by a special [MASK] symbol and recovered through iterative denoising. MDLM [28] adapted this framework to language modeling via a continuous-time weighted cross-entropy objective, demonstrating its effectiveness as a non-autoregressive alternative for conditional text generation. Building on MDLM, subsequent works have improved robustness through unmasking-order-aware training [15] and conditional fidelity by treating observed tokens as fixed anchors during sampling [17, 27]. Despite these advances, masked diffusion has been studied almost exclusively in discrete modalities such as NLP, and its extension to time series remains underexplored due to the fundamental gap between discrete tokens and continuous, ordinal temporal signals. Figure 2: Overview of the proposed Stochastic Discretization and Ordinal-Aware Soft Labeling framework. The process consists of three stages: (1) Masking Input Data, where missing values are identified; (2) Stochastic Discretization, which applies instance-adaptive normalization and injects stochastic noise ϵ∼(−0.5,0.5)ε (-0.5,0.5) to bridge the continuous-discrete gap; and (3) Ordinal-Aware Soft Labeling, which maps continuous values to discrete tokens while generating soft target distributions to preserve ordinal semantic relationships. 3 Discrete Representation for Time-Series Modeling Although the Masked Diffusion Model (MDM) offers a promising framework for imputation, directly applying it to time series presents a structural challenge: MDMs are inherently designed for discrete state spaces (e.g., text tokens), whereas time series data consist of continuous numerical values. To bridge this gap, we introduce a unified discretization framework illustrated in Figure 2, which combines Stochastic Discretization to mitigate quantization error and Ordinal-Aware Soft Labeling to preserve the temporal order of the original signals within the discrete vocabulary. 3.1 Stochastic Discretization Instance-Adaptive Normalization. To stabilize local statistics against distribution shifts, we normalize each input window x∈ℝT×Cx ^T× C into x~ x. Specifically, the normalization statistics are computed solely using the observed values (excluding masked positions), ensuring that there is no data-leakage in the training and imputation process. Stochastic Token Generation. We map the normalized values to a discrete vocabulary of size K. Specifically, we first project the continuous value x~t,c x_t,c onto the grid coordinates and inject uniform noise ϵ∼(−0.5,0.5)ε (-0.5,0.5) before rounding: ct,c=x~t,c−vminvmax−vmin⋅(K−1)+1,zt,c=⌊ct,c+ϵ⌉,c_t,c= x_t,c-v_ v_ -v_ ·(K-1)+1, z_t,c= c_t,c+ε , (1) where ⌊⋅⌉ · denotes the nearest integer function. This stochastic mechanism preserves information in expectation. For instance, a coordinate ct,c=3.6c_t,c=3.6 is assigned to token 44 with 60%60\% probability and token 33 with 40%40\%, ensuring [zt,c]≈ct,cE[z_t,c]≈ c_t,c. This allows the model to learn precise dynamics beyond the fixed grid resolution. Figure 3: Overview of the Masked Diffusion Time-series Imputation Model (MDTIM) framework. The architecture operates on the discretized input through three key phases: (a) Discretized Token Embedding: The discrete tokens generated by the stochastic tokenizer are mapped into dense vector representations. (b) Ordinal-Spectral Consistent Training: The backbone is trained with a dual objective combining Ordinal-Aware Soft Cross-Entropy and Spectral Consistency Loss. This ensures spectral consistency across global frequencies while maintaining local ordinal accuracy. (c) Continuous Reconstruction via Expectation: Finally, continuous values are recovered by computing the probability-weighted expectation of the predicted token distribution, enabling precise dense reconstruction from discrete outputs. 3.2 Ordinal-Aware Soft Labeling Treating discretized tokens as independent classes ignores their ordinal nature (e.g., token k is semantically closer to k+1k+1 than k+10k+10). To enforce continuity, we employ Ordinal-Aware Soft Labeling with a truncated Gaussian kernel. Given a ground-truth token index y∈1,…,Ky∈\1,…,K\, the target probability sis_i for the i-th class is defined as: si=1Zexp(−(i−y)2σ2)if |i−y|≤w,0otherwise,s_i= cases 1Z (- (i-y)^2σ^2 )&if |i-y|≤ w,\\ 0&otherwise, cases (2) where w is the truncation window size (e.g., w=2w=2), σ controls smoothness, and Z is the normalization constant. Crucially, we assign zero probability to the [MASK] token (s0=0s_0=0) and distant bins. This concentrates probability mass on valid ordinal neighbors while preserving the structural separation of the mask state. 4 Masked Diffusion Framework Building upon the discrete representation of time series data described in Section 3, we introduce the Masked Diffusion Time-series Imputation Model (MDTIM), a framework designed to model the joint distribution of multivariate time series and incorporate training techniques of masked diffusion models. Figure 3 illustrates the overall pipeline. 4.1 Factorized Temporal-Variate Backbone To effectively model the joint distribution of multivariate time series, we employ a Factorized Temporal-Variate Transformer based on the Diffusion Transformer (DiT) proposed by 25. This architecture processes the input through alternating factorization of the time (T) and feature (C) axes. Discretized Embedding. Given the masked input indices t∈0,…,KB×T×Cz_t∈\0,…,K\^B× T× C at diffusion step t, we utilize a distinct embedding matrix cE_c for each channel c to preserve semantic orthogonality. The input is projected to an initial hidden state 0∈ℝB×T×C×DH^0 ^B× T× C× D, where D is the hidden dimension. Interleaved DiT Blocks. The core network consists of L layers, each sequentially processing the temporal and feature axes via two DiT blocks conditioned on the diffusion timestep t, with Rotary Position Embeddings (RoPE) [32] encoding relative positions. The Temporal block reshapes (l−1)H^(l-1) to ℝ(B⋅C)×T×DR^(B· C)× T× D and captures time-dependencies along the temporal axis: ′(l)=TemporalEncoder((l−1),t).H (l)=TemporalEncoder(H^(l-1),t). (3) The Variate block then reshapes ′(l)H (l) to ℝ(B⋅T)×C×DR^(B· T)× C× D and symmetrically models cross-channel correlations: (l)=VariateEncoder(′(l),t).H^(l)=VariateEncoder(H (l),t). (4) The final output is projected to vocabulary size Kout=1.5K_out=1.5K, covering the extended range [−1.5,1.5][-1.5,1.5] to account for potential distribution shifts in unobserved regions. 4.2 Ordinal-Spectral Consistent Training Objective To simultaneously ensure local reconstruction fidelity and global temporal coherence, we formulate a dual-domain objective. This combines a discrete diffusion loss adapted for ordinal continuity with an auxiliary spectral regularization term. Ordinal-Aware Masked Diffusion. We adopt the continuous-time training paradigm of Masked Diffusion Language Models (MDLM) [2, 28], optimizing a weighted variational lower bound (NELBO). Standard MDLM employs a simple cross-entropy loss against one-hot targets, which treats all incorrect tokens equally, failing to capture the ordinal magnitude inherent in time-series data. To address this, we employ the Ordinal-Aware Soft Labels derived in Sec. 3.2 as the optimization targets. For a set of masked indices ℳtM_t at timestep t, the diffusion loss is formulated as: ℒdiff=t[w(t)∑l∈ℳt−⟨(l),logθ(t(l))⟩],L_diff=E_t [w(t) _l _t- ^(l), _θ(z_t^(l)) ], (5) where logθ(t(l))∈ℝK _θ(z_t^(l)) ^K is the predicted log-probability vector, and (l)∈ℝKs^(l) ^K denotes the soft target distribution. The detailed formulation of time-dependent weight w(t)w(t) and the corresponding noise schedule are detailed in Appendix C.2.2. Replacing the one-hot target ye_y with (l)s^(l) can be viewed as a smooth relaxation of the standard MDLM objective. Writing (l)=(1−ϵ)y+ϵ(l)s^(l)=(1-ε)e_y+ε\,u^(l) with ϵ=1−sy(l)ε=1-s_y^(l), the per-token loss decomposes as: −⟨(l),logθ⟩=(1−ϵ)ℒNELBO(l)+ϵDKL((l)∥θ)+const,- ^(l), _θ =(1-ε)L_NELBO^(l)+ε D_KL(u^(l)\|p_θ)+const, (6) so that the original NELBO is recovered when (l)s^(l) collapses to a one-hot (ϵ→0ε→ 0) and otherwise encourages the predicted distribution to align with the ordinal neighborhood of y. The detailed derivation is provided in Appendix B. Spectral Consistency Regularization. While the token-wise objective ensures local accuracy, it may neglect global temporal correlations and frequency structures. To strictly enforce global coherence, we incorporate a Spectral Consistency Loss, adapting the frequency-domain regularization proposed in DiffusionTS [41]. During training, we compute a continuous estimate of the full sequence, ^t x_t, by applying expectation-based decoding (Sec. 4.3). We then minimize the L1L_1 distance between the Fourier representations of the reconstructed and ground-truth signals: ℒFFT=t[‖ℱ(^t)−ℱ(0)‖1],L_FFT=E_t [\|F( x_t)-F(x_0)\|_1 ], (7) where ℱ(⋅)F(·) denotes the Fast Fourier Transform (FFT). Therefore, the total training objective is a weighted combination: ℒtotal=ℒdiff+λℒFFTL_total=L_diff+ _FFT. 4.3 Continuous Reconstruction via Expectation Since Stochastic Discretization injects random noise into the input, the model output depends on the specific noise realization. To obtain a robust estimate, we run inference M times (e.g., M=10M=10) with independent noise samples and average the predicted probability distributions, p¯t,c=1M∑m=1Mpt,c(m) p_t,c= 1M _m=1^Mp_t,c^(m). The continuous value is then recovered as the expectation over the averaged distribution: x^t,c=∑k=1Kp¯t,c(k)⋅vk, x_t,c= _k=1^K p_t,c^(k)· v_k, (8) where vkv_k is the center value of the k-th bin within [−1.5,1.5][-1.5,1.5]. The result is then denormalized using the instance-wise statistics to recover the original scale. Table 1: Quantitative comparison of multivariate time-series imputation performance (L=48L=48). We report the MAE of MDTIM and baselines averaged over 3 random seeds at missing ratios of 30% and 70% (50% in Appendix D.1). The best results are highlighted in bold. Dataset Energy ETTh Weather Sine Missing Type Uniform Geometric Uniform Geometric Uniform Geometric Uniform Geometric Model 30% 70% 30% 70% 30% 70% 30% 70% 30% 70% 30% 70% 30% 70% 30% 70% RNN BRITS 0.246 0.379 0.329 0.366 0.194 0.299 0.227 0.292 0.050 0.073 0.058 0.071 0.010 0.021 0.020 0.019 MRNN 1.086 1.143 1.087 1.143 0.743 0.782 0.751 0.781 0.651 0.661 0.654 0.659 0.169 0.170 0.169 0.170 GRUD 0.364 0.457 0.426 0.445 0.310 0.394 0.348 0.383 0.104 0.369 0.164 0.350 0.008 0.014 0.013 0.012 Transformer Transformer 0.323 0.401 0.354 0.393 0.168 0.252 0.183 0.247 0.076 0.070 0.083 0.068 0.039 0.046 0.042 0.045 Informer 0.344 0.379 0.371 0.373 0.207 0.302 0.226 0.296 0.047 0.059 0.052 0.058 0.027 0.039 0.036 0.037 PatchTST 0.586 0.529 0.528 0.517 0.202 0.272 0.231 0.262 0.074 0.082 0.081 0.077 0.015 0.017 0.019 0.015 SAITS 0.177 0.212 0.204 0.208 0.140 0.210 0.152 0.206 0.045 0.049 0.050 0.048 0.026 0.026 0.030 0.025 ImputeFormer 0.066 0.219 0.113 0.194 0.146 0.245 0.165 0.235 0.049 0.171 0.090 0.149 0.006 0.041 0.012 0.037 CNN TimesNet 0.617 0.818 0.689 0.808 0.593 0.719 0.628 0.717 0.211 0.351 0.225 0.348 0.167 0.220 0.189 0.217 SCINet 0.557 0.571 0.546 0.565 0.238 0.323 0.261 0.317 0.074 0.089 0.086 0.087 0.019 0.031 0.026 0.029 Linear DLinear 0.795 0.527 0.684 0.517 0.379 0.445 0.373 0.435 0.372 0.223 0.345 0.220 0.070 0.053 0.062 0.051 FiLM 0.877 0.520 0.791 0.512 0.696 0.627 0.707 0.623 0.380 0.209 0.360 0.206 0.127 0.109 0.132 0.108 FreTS 0.157 0.219 0.225 0.206 0.222 0.299 0.262 0.286 0.071 0.080 0.080 0.074 0.068 0.068 0.075 0.065 Generative GPVAE 0.476 0.730 0.501 0.727 0.333 0.449 0.369 0.440 0.158 0.255 0.172 0.252 0.152 0.162 0.153 0.162 USGAN 0.264 0.412 0.330 0.401 0.208 0.302 0.240 0.295 0.087 0.122 0.101 0.118 0.014 0.029 0.026 0.026 CSDI 0.094 0.139 0.110 0.136 0.160 0.244 0.177 0.240 0.039 0.049 0.044 0.048 0.003 0.004 0.004 0.004 FGTI 0.050 0.100 0.065 0.097 0.218 0.349 0.293 0.328 0.038 0.051 0.046 0.049 0.001 0.003 0.002 0.003 MDTIM (Ours) 0.044 0.085 0.053 0.082 0.127 0.211 0.146 0.205 0.032 0.044 0.036 0.043 0.002 0.004 0.003 0.003 5 Experiments We empirically present a comprehensive evaluation of the proposed Masked Diffusion Time-series Imputation Model (MDTIM), comparing it against state-of-the-art deterministic and generative baselines. Our experiments focus on imputation accuracy, robustness to complex missing patterns, and computational efficiency. 5.1 Experimental Setup We evaluate our method on widely used benchmarks: Energy, ETTh, Weather, and the synthetic Sine dataset. To rigorously test the models, we simulate data corruption using two distinct mechanisms: Uniform masking, where observations are dropped randomly to mimic independent failures, and Geometric masking, which removes consecutive time steps to simulate prolonged sensor malfunctions or transmission errors. For both scenarios, we assess performance across varying missing rates of 30%, 50%, and 70%. We compare MDTIM against a comprehensive suite of 19 baselines spanning various architectural paradigms. These include RNN-based methods such as BRITS[3], MRNN[40], and GRU-D[4]; and Transformer-based models including the canonical Transformer[34], Informer[43], PatchTST[24], SAITS[6], and ImputeFormer[23]. Also, we include recent high-performance convolutional and linear architectures: TimesNet[36], SCINet[19], DLinear[42], FILM[44], and FreTS[39]. Finally, we evaluate against probabilistic frameworks including GP-VAE[9], US-GAN[22], and the diffusion-based CSDI[33] and FGTI[38]. To account for the stochastic nature of the missing data patterns, we evaluate each trained model across three independent inference trials using distinct random seeds for mask generation, and report the averaged Mean Absolute Error (MAE). 5.2 Imputation Performance Main Results. As shown in Table 1, MDTIM achieves the lowest MAE on Energy, ETTh, and Weather across both Uniform and Geometric missing patterns. FGTI [38], a frequency-aware diffusion model, attains the best results on the synthetic Sine dataset whose signal is dominated by a few well-defined frequencies, while MDTIM trails by only a marginal gap (e.g., 0.002 vs. 0.001 at 30% Uniform). On the more complex real-world ETTh dataset, however, FGTI degrades substantially (0.218 vs. 0.127 at 30% Uniform), suggesting that frequency-domain priors generalize less reliably to irregular and non-stationary signals. Compared to the diffusion-based CSDI [33], MDTIM consistently yields lower MAE across all datasets, owing to our Expectation-based Reconstruction that recovers continuous values deterministically from the predicted token distribution and mitigates the sampling variance inherent in continuous diffusion. MDTIM also exhibits notable robustness under geometric masking, where most baselines suffer significant degradation; for instance, on Energy at 30%, MDTIM increases only from 0.044 to 0.053, indicating that the masked diffusion objective encourages inference of global temporal structure rather than reliance on local interpolation. Scalability to Long-Term Dependencies. To verify that MDTIM captures long-range temporal dependencies, we extend the evaluation to longer sequence lengths L∈96,192L∈\96,192\. As shown in Table 2, the performance advantage of MDTIM grows with the sequence length: while MDTIM is comparable to SAITS at L=48L=48, the gap widens substantially at longer horizons. At L=192L=192 with 30% Uniform missing, MDTIM achieves an MAE of 0.123, considerably lower than BRITS (0.192) and SAITS (0.145), indicating that MDTIM captures global context without suffering from the attention dilution often observed in baselines. Table 2: Robustness analysis on extended sequence lengths (L∈96,192L∈\96,192\). We report the MAE averaged over 3 random seeds. Length 96 192 Missing Type Uniform Geometric Uniform Geometric Model 30% 50% 70% 30% 50% 70% 30% 50% 70% 30% 50% 70% BRITS 0.191 0.229 0.287 0.220 0.250 0.280 0.192 0.230 0.287 0.222 0.250 0.280 SAITS 0.139 0.162 0.204 0.150 0.172 0.200 0.145 0.162 0.198 0.156 0.171 0.195 CSDI 0.149 0.177 0.222 0.162 0.187 0.218 0.170 0.199 0.245 0.183 0.209 0.241 MDTIM 0.124 0.152 0.199 0.138 0.163 0.194 0.123 0.148 0.190 0.135 0.157 0.186 Probabilistic Imputation Performance Beyond point estimation accuracy, we evaluate the quality of the predictive distributions generated by the models using the Continuous Ranked Probability Score (CRPS). CRPS estimates the calibration (reliability) and sharpness (precision) of the probabilistic inference. Table 3: Quantitative comparison of probabilistic imputation performance (CRPS, L=48L=48). We compare MDTIM with the diffusion-based baseline CSDI. Lower is better. Energy ETTh Sine Missing Type Rate CSDI MDTIM CSDI MDTIM CSDI MDTIM Uniform 30% 0.0629 0.0429 0.1120 0.0863 0.0023 0.0020 50% 0.0747 0.0541 0.1356 0.1078 0.0025 0.0020 70% 0.0927 0.0771 0.1748 0.1520 0.0030 0.0020 Geometric 30% 0.0732 0.0547 0.1262 0.1017 0.0027 0.0020 50% 0.0818 0.0626 0.1460 0.1197 0.0028 0.0020 70% 0.0905 0.0733 0.1708 0.1462 0.0029 0.0020 As shown in Table 3, MDTIM consistently outperforms CSDI, the diffusion-based baseline, across all datasets and missing scenarios. On the highly periodic Sine dataset, MDTIM achieves a near-constant CRPS regardless of the missing rate, indicating that it captures deterministic periodicity with high confidence. While continuous diffusion models often exhibit excessive variance due to sampling noise, MDTIM mitigates such ambiguity through its expectation-based decoding over discrete distributions. Furthermore, on the complex real-world ETTh dataset, MDTIM achieves significantly lower CRPS (e.g., 0.0863 vs. 0.1120 at 30% uniform missing), demonstrating that our Ordinal-Aware Soft Labeling guides the model to produce well-calibrated uncertainty estimates. Imputation under Naturally Missing Real-World Data. To assess MDTIM under realistic missingness patterns, we further evaluate on PhysioNet 2012 [29], a multivariate clinical time series benchmark with ∼ 80% natural missingness due to irregular ICU sampling, and apply additional uniform and geometric masking on top of the existing missingness. Since the highly sparse and irregular nature of PhysioNet yields no dominant spectral structure, we disable the spectral consistency loss for this dataset; all other settings remain identical. Table 4: Imputation performance on PhysioNet 2012. PhysioNet2012 Uniform Geometric Model 30% 50% 70% 30% 50% 70% BRITS 0.333 0.364 0.408 0.356 0.379 0.405 SAITS 0.283 0.314 0.360 0.298 0.324 0.358 CSDI 0.303 0.332 0.376 0.320 0.342 0.374 ImputeFormer 0.289 0.329 0.389 0.311 0.343 0.388 FGTI 0.338 0.376 0.425 0.370 0.397 0.423 MDTIM (Ours) 0.236 0.279 0.339 0.257 0.291 0.337 As shown in Table 4, MDTIM achieves the lowest MAE across all settings, outperforming SAITS by a clear margin. Notably, FGTI exhibits the largest degradation here, even underperforming BRITS. This indicates that frequency-domain priors become unreliable for highly irregular signals where dominant spectral components are weak or unstable, whereas our ordinal-aware discrete formulation generalizes effectively to real-world sparse and irregular observations. 5.3 Computational Efficiency and Model Scalability Diffusion-based methods often trade efficiency for accuracy. To examine this trade-off, we compare MDTIM against representative baselines at two model scales on the Energy dataset with uniform missing (Table 5). Table 5: Imputation performance and efficiency on Energy (uniform missing). Inference time is measured on the full test set (MDTIM: M=10M=10 noise samples; CSDI: T=50T=50 diffusion steps). Scale Small Large Model 30% 50% 70% Params (M) Time (s) 30% 50% 70% Params (M) Time (s) BRITS 0.246 0.296 0.379 2.18 4.88 0.245 0.293 0.371 8.55 4.96 SAITS 0.212 0.216 0.240 25.25 0.41 0.177 0.185 0.212 88.24 0.79 CSDI 0.094 0.112 0.139 1.19 401.69 0.074 0.093 0.128 4.49 894.76 MDTIM (Ours) 0.046 0.065 0.092 0.65 4.05 0.044 0.061 0.085 8.33 12.63 MDTIM addresses the latency bottleneck of continuous diffusion: while CSDI requires up to 894s at the Large scale, the Small MDTIM completes inference in 4.05s while halving the MAE (0.074→ 0.046). It is also notably parameter-efficient—the 0.65M Small variant surpasses the 88.24M Large SAITS, suggesting that the masked diffusion objective captures temporal dynamics more densely than scaling deterministic Transformers. 5.4 Downstream Forecasting Task Table 6: Forecasting MAE on Energy (L=48L=48) using a fixed PatchTST predictor with imputed inputs from each model. Energy Forecast Horizon Imputer H=4H=4 H=6H=6 H=8H=8 H=12H=12 SAITS 0.187 0.207 0.224 0.254 CSDI 0.171 0.190 0.208 0.237 MDTIM (Ours) 0.138 0.153 0.170 0.200 To assess whether the imputation accuracy of MDTIM translates to downstream tasks, we evaluate forecasting performance when the input window contains missing values. Each imputation model reconstructs the partially observed input, after which a fixed PatchTST [24] forecaster, pre-trained on clean sequences, predicts the subsequent H steps. Since the forecaster is identical across all settings, differences in forecasting MAE directly reflect the quality of imputation. As shown in Table 6, MDTIM yields substantially lower forecasting error than both baselines across all horizons on Energy, reducing MAE by 19–26% over SAITS and 16–19% over CSDI. These results indicate that the imputation accuracy of MDTIM propagates to downstream forecasting, and that high-fidelity imputation is essential for tasks that rely on partially observed inputs. 5.5 Ablation Studies To rigorously evaluate the proposed framework, we conduct ablation studies on the contribution of our expectation-based unmasking and soft labeling strategies, and the sensitivity to the vocabulary size. Impact of Unmasking and Labeling Strategies. Table 7: Ablation study on unmasking and labeling strategies. All results are reported as MAE averaged over 3 random seeds under the uniform missing scenario. Components ETTh Energy Sine Exp. Unmasking Soft-Label 30% 50% 70% 30% 50% 70% 30% 50% 70% ✗ ✗ 0.134 0.166 0.229 0.053 0.073 0.101 0.0071 0.0072 0.0075 ✗ ✓ 0.131 0.162 0.222 0.056 0.079 0.115 0.0069 0.0070 0.0074 ✓ ✗ 0.129 0.159 0.212 0.044 0.060 0.085 0.0026 0.0030 0.0036 ✓ ✓ 0.127 0.157 0.211 0.044 0.061 0.085 0.0022 0.0026 0.0035 Table 7 summarizes the effect of our two proposed components: Expectation-based Unmasking and Ordinal-Aware Soft Labeling. Expectation-based Unmasking consistently reduces error across all datasets, confirming that computing the expected value over discrete bins mitigates quantization error inherent in Argmax selection. The contribution of Soft Labeling, in contrast, scales with the structural regularity of the data. On the highly periodic Sine dataset, Soft Labeling yields a substantial additional gain on top of Expectation-based Unmasking, as ordinal-aware guidance aligns naturally with the smooth, predictable transitions of periodic signals. A similar but milder gain is observed on ETTh, whose complex temporal dynamics still benefit from ordinal regularization under higher uncertainty. On Energy, where the dynamics are less structured, Soft Labeling offers no further improvement beyond Expectation-based Unmasking. These results indicate that Soft Labeling is most effective when the underlying signal exhibits clear ordinal or periodic structure that ordinal-aware guidance can exploit. Sensitivity to Vocabulary Size (K). Table 8: MAE on ETTh with varying vocabulary size (K). ETTh Uniform Geometric Bins 30% 50% 70% 30% 50% 70% 20 0.138 0.166 0.217 0.157 0.181 0.211 40 0.127 0.157 0.211 0.146 0.173 0.205 60 0.129 0.160 0.214 0.148 0.175 0.208 We further investigate the impact of the vocabulary size K on imputation accuracy (Table 8), which exhibits a trade-off between discretization resolution and classification complexity. A small vocabulary (K=20K=20) suffers from high quantization error, failing to capture fine-grained fluctuations. Conversely, K=60K=60 slightly degrades performance, as a larger vocabulary increases the difficulty of the discrete classification task. K=40K=40 achieves the best balance, providing sufficient resolution to approximate the continuous signal while remaining tractable to optimize. 6 Conclusion and Limitations We proposed MDTIM, a framework that adapts the masked diffusion training paradigm to continuous time series imputation. By introducing Stochastic Discretization and Ordinal-Aware Soft Labeling, MDTIM bridges the gap between continuous dynamics and categorical tokenization, while the Expectation-based Unmasking strategy enables precise continuous reconstruction. Extensive experiments demonstrate that MDTIM consistently outperforms both deterministic and generative baselines in reconstruction accuracy and robustness across diverse missing scenarios, while requiring substantially less inference time than continuous diffusion baselines. Limitations and Future Work. While MDTIM demonstrates strong performance across diverse missing scenarios, several aspects warrant further investigation. First, the optimal vocabulary size K may vary across datasets depending on the range and granularity of the underlying signal; a data-adaptive discretization scheme that automatically selects the resolution would further improve generalization. Second, our spectral consistency regularization assumes the presence of stable frequency components and is therefore disabled on benchmarks dominated by irregular sampling such as PhysioNet. Designing an adaptive frequency objective that modulates its influence based on local signal regularity is a promising direction. Finally, our evaluation focuses on standard imputation benchmarks of moderate length; extending MDTIM to extremely long sequences and to multi-scale or non-stationary domains such as financial tick data remains for future work. References [1] J. M. L. Alcaraz and N. Strodthoff (2022) Diffusion-based time series imputation and forecasting with structured state space models. arXiv preprint arXiv:2208.09399. Cited by: §A.1.2, §A.1.3, §A.1. [2] J. Austin, D. D. Johnson, J. Ho, D. Tarlow, and R. Van Den Berg (2021) Structured denoising diffusion models in discrete state-spaces. Advances in neural information processing systems 34, p. 17981–17993. Cited by: §A.2, §2.2, §4.2. [3] W. Cao, D. Wang, J. Li, H. Zhou, L. Li, and Y. Li (2018) Brits: bidirectional recurrent imputation for time series. Advances in neural information processing systems 31. Cited by: §2.1, §5.1. [4] Z. Che, S. Purushotham, K. Cho, D. Sontag, and Y. Liu (2018) Recurrent neural networks for multivariate time series with missing values. Scientific reports 8 (1), p. 6085. Cited by: §2.1, §2.1, §5.1. [5] P. Dhariwal and A. Nichol (2021) Diffusion models beat gans on image synthesis. Advances in neural information processing systems 34, p. 8780–8794. Cited by: §A.1.2, §A.1. [6] W. Du, D. Côté, and Y. Liu (2023) Saits: self-attention-based imputation for time series. Expert Systems with Applications 219, p. 119619. Cited by: §1, §2.1, §2.1, §5.1. [7] W. Du, J. Wang, L. Qian, Y. Yang, Z. Ibrahim, F. Liu, Z. Wang, H. Liu, Z. Zhao, Y. Zhou, et al. (2024) Tsi-bench: benchmarking time series imputation. arXiv preprint arXiv:2406.12747. Cited by: §C.1. [8] G. Fadlon, I. Arbiv, N. Berman, and O. Azencot (2025) A diffusion model for regular time series generation from irregular data with completion and masking. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: §1. [9] V. Fortuin, D. Baranchuk, G. Rätsch, and S. Mandt (2020) Gp-vae: deep probabilistic time series imputation. In International conference on artificial intelligence and statistics, p. 1651–1661. Cited by: §5.1. [10] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. Advances in neural information processing systems 33, p. 6840–6851. Cited by: §A.1.1, §A.1.3, §A.1. [11] J. Ho, T. Salimans, A. Gritsenko, W. Chan, M. Norouzi, and D. J. Fleet (2022) Video diffusion models. Advances in neural information processing systems 35, p. 8633–8646. Cited by: §A.1. [12] J. Ho and T. Salimans (2022) Classifier-free diffusion guidance. arXiv preprint arXiv:2207.12598. Cited by: §A.1.2. [13] E. Hoogeboom, D. Nielsen, P. Jaini, P. Forr, and M. Welling (2021) Argmax flows and multinomial diffusion: learning categorical distributions. Advances in neural information processing systems 34, p. 12454–12465. Cited by: §2.2. [14] T. Karras, M. Aittala, T. Aila, and S. Laine (2022) Elucidating the design space of diffusion-based generative models. Advances in neural information processing systems 35, p. 26565–26577. Cited by: §A.1. [15] J. Kim, K. Shah, V. Kontonis, S. M. Kakade, and S. Chen (2025) Train for the worst, plan for the best: understanding token ordering in masked diffusions. In Forty-second International Conference on Machine Learning, Cited by: §2.2. [16] S. Kim, H. Kim, E. Yun, H. Lee, J. Lee, and J. Lee (2023) Probabilistic imputation for time-series classification with missing data. In International Conference on Machine Learning, p. 16654–16667. Cited by: §1, §2.1. [17] H. Koh, M. Jhang, D. Kim, S. Lee, and K. Jung (2025) Conditional [mask] discrete diffusion language model. In Proceedings of the 2025 Conference on Empirical Methods in Natural Language Processing, p. 8910–8934. Cited by: §2.2. [18] Z. Kong, W. Ping, J. Huang, K. Zhao, and B. Catanzaro (2020) Diffwave: a versatile diffusion model for audio synthesis. arXiv preprint arXiv:2009.09761. Cited by: §A.1. [19] M. Liu, A. Zeng, M. Chen, Z. Xu, Q. Lai, L. Ma, and Q. Xu (2022) Scinet: time series modeling and forecasting with sample convolution and interaction. Advances in Neural Information Processing Systems 35, p. 5816–5828. Cited by: §5.1. [20] J. M. Lopez Alcaraz and N. Strodthoff (2023) Diffusion-based time series imputation and forecasting with structured atate apace models. Transactions on machine learning research, p. 1–36. Cited by: §2.1. [21] J. E. Matheson and R. L. Winkler (1976) Scoring rules for continuous probability distributions. Management science 22 (10), p. 1087–1096. Cited by: §C.3. [22] X. Miao, Y. Wu, J. Wang, Y. Gao, X. Mao, and J. Yin (2021) Generative semi-supervised learning for multivariate time series imputation. In Proceedings of the AAAI conference on artificial intelligence, Vol. 35, p. 8983–8991. Cited by: §5.1. [23] T. Nie, G. Qin, W. Ma, Y. Mei, and J. Sun (2024) ImputeFormer: low rankness-induced transformers for generalizable spatiotemporal imputation. In Proceedings of the 30th ACM SIGKDD conference on knowledge discovery and data mining, p. 2260–2271. Cited by: §5.1. [24] Y. Nie, N. H. Nguyen, P. Sinthong, and J. Kalagnanam (2023) A time series is worth 64 words: long-term forecasting with transformers. In The Eleventh International Conference on Learning Representations, External Links: Link Cited by: §5.1, §5.4. [25] W. Peebles and S. Xie (2023) Scalable diffusion models with transformers. In Proceedings of the IEEE/CVF international conference on computer vision, p. 4195–4205. Cited by: §4.1. [26] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer (2022) High-resolution image synthesis with latent diffusion models. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, p. 10684–10695. Cited by: §A.1. [27] L. Rout, C. Caramanis, and S. Shakkottai (2025) Anchored diffusion language model. Advances in Neural Information Processing Systems. Cited by: §2.2. [28] S. Sahoo, M. Arriola, Y. Schiff, A. Gokaslan, E. Marroquin, J. Chiu, A. Rush, and V. Kuleshov (2024) Simple and effective masked diffusion language models. Advances in Neural Information Processing Systems 37, p. 130136–130184. Cited by: §A.2, Appendix B, §1, §2.2, §4.2. [29] I. Silva, G. Moody, D. J. Scott, L. A. Celi, and R. G. Mark (2012) Predicting in-hospital mortality of icu patients: the physionet/computing in cardiology challenge 2012. In 2012 computing in cardiology, p. 245–248. Cited by: §5.2. [30] J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli (2015) Deep unsupervised learning using nonequilibrium thermodynamics. In International conference on machine learning, p. 2256–2265. Cited by: §A.1. [31] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole (2020) Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456. Cited by: §A.1. [32] J. Su, M. Ahmed, Y. Lu, S. Pan, W. Bo, and Y. Liu (2024) Roformer: enhanced transformer with rotary position embedding. Neurocomputing 568, p. 127063. Cited by: §4.1. [33] Y. Tashiro, J. Song, Y. Song, and S. Ermon (2021) Csdi: conditional score-based diffusion models for probabilistic time series imputation. Advances in neural information processing systems 34, p. 24804–24816. Cited by: §A.1.2, §A.1.3, §A.1, §1, §2.1, §5.1, §5.2. [34] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin (2017) Attention is all you need. Advances in neural information processing systems 30. Cited by: §5.1. [35] J. Wang, W. Du, Y. Yang, L. Qian, W. Cao, K. Zhang, W. Wang, Y. Liang, and Q. Wen (2025) Deep learning for multivariate time series imputation: a survey. In Proceedings of the Thirty-Fourth International Joint Conference on Artificial Intelligence, IJCAI ’25. External Links: ISBN 978-1-956792-06-5, Link, Document Cited by: §2.1. [36] H. Wu, T. Hu, Y. Liu, H. Zhou, J. Wang, and M. Long (2023) TimesNet: temporal 2d-variation modeling for general time series analysis. In The Eleventh International Conference on Learning Representations, External Links: Link Cited by: §5.1. [37] H. Wu, J. Xu, J. Wang, and M. Long (2021) Autoformer: decomposition transformers with auto-correlation for long-term series forecasting. Advances in neural information processing systems 34, p. 22419–22430. Cited by: §C.1. [38] X. Yang, Y. Sun, X. Yuan, and X. Chen (2024) Frequency-aware generative models for multivariate time series imputation. Advances in Neural Information Processing Systems 37, p. 52595–52623. Cited by: §5.1, §5.2. [39] K. Yi, Q. Zhang, W. Fan, S. Wang, P. Wang, H. He, N. An, D. Lian, L. Cao, and Z. Niu (2023) Frequency-domain mlps are more effective learners in time series forecasting. Advances in Neural Information Processing Systems 36, p. 76656–76679. Cited by: §5.1. [40] J. Yoon, W. R. Zame, and M. Van Der Schaar (2018) Estimating missing data in temporal data streams using multi-directional recurrent neural networks. IEEE Transactions on Biomedical Engineering 66 (5), p. 1477–1490. Cited by: §5.1. [41] X. Yuan and Y. Qiao (2024) Diffusion-TS: interpretable diffusion for general time series generation. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: §4.2. [42] A. Zeng, M. Chen, L. Zhang, and Q. Xu (2023) Are transformers effective for time series forecasting?. In Proceedings of the AAAI conference on artificial intelligence, Vol. 37, p. 11121–11128. Cited by: §5.1. [43] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang (2021) Informer: beyond efficient transformer for long sequence time-series forecasting. In Proceedings of the AAAI conference on artificial intelligence, Vol. 35, p. 11106–11115. Cited by: §C.1, §5.1. [44] T. Zhou, Z. Ma, Q. Wen, L. Sun, T. Yao, W. Yin, R. Jin, et al. (2022) Film: frequency improved legendre memory model for long-term time series forecasting. Advances in neural information processing systems 35, p. 12677–12690. Cited by: §5.1. Appendix A Preliminaries on Diffusion Models A.1 Diffusion Models Diffusion models [30, 10, 31] are a class of latent-variable generative models that learn to invert a fixed noise-injection process. Since their introduction, they have become the dominant paradigm in many continuous-data generative tasks, achieving state-of-the-art results in image [5, 26, 14], audio [18], and video [11] synthesis, and have also been adopted as strong probabilistic models for time series forecasting and imputation [33, 1]. Notation. Let x0∈ℝdx_0 ^d denote a clean data sample drawn from the data distribution q(x0)q(x_0), and let xt∈ℝdx_t ^d for t∈1,…,Tt∈\1,…,T\ denote its noisy latents along the forward chain. A noise schedule βtt=1T⊂(0,1)\ _t\_t=1^T⊂(0,1) controls the amount of corruption at each step, with the per-step retention factor αt=1−βt _t=1- _t and the cumulative retention factor α¯t=∏s≤tαs α_t= _s≤ t _s; by construction α¯0=1 α_0=1 and α¯T≈0 α_T≈ 0. A.1.1 Diffusion Framework Forward (Gaussian corruption). The forward process gradually corrupts a data sample into pure Gaussian noise over T steps: q(xt∣xt−1)=(xt,1−βtxt−1,βtI).q(x_t x_t-1)=N(x_t; 1- _tx_t-1, _tI). (9) A useful property of this Markov chain is that its t-step marginal admits a closed form: q(xt|x0)=(xt,α¯tx0,(1−α¯t)I),xt=α¯tx0+1−α¯tϵ,ϵ∼(0,I).q(x_t|x_0)=N(x_t; α_tx_0,(1- α_t)I),\ x_t= α_tx_0+ 1- α_tε,\ ε (0,I). (10) That is, xtx_t at any timestep t can be sampled directly from x0x_0 without simulating the full chain, which makes per-step training tractable. As t→Tt→ T, the schedule is chosen so that α¯T≈0 α_T≈ 0 and xTx_T becomes approximately isotropic Gaussian, independent of x0x_0. Reverse (learned denoising). The reverse process incrementally denoises noise back into the data manifold. While the marginal q(xt−1|xt)q(x_t-1|x_t) is intractable, the conditional posterior q(xt−1|xt,x0)q(x_t-1|x_t,x_0) is Gaussian with closed-form mean and variance depending only on βt _t and α¯t α_t [10]. A neural network parameterizes the reverse kernel pθ(xt−1|xt)=(xt−1,μθ(xt,t),Σθ(xt,t)),p_θ(x_t-1|x_t)=N(x_t-1; _θ(x_t,t), _θ(x_t,t)), (11) where the mean is typically expressed through a noise predictor ϵθ(xt,t) _θ(x_t,t): μθ(xt,t)=1αt(xt−βt1−α¯tϵθ(xt,t)). _θ(x_t,t)= 1 _t (x_t- _t 1- α_t _θ(x_t,t) ). (12) A.1.2 Conditional Diffusion The diffusion framework introduced in Appendix A.1 defines an unconditional generative model. Many practical applications, including the time series imputation setting of this paper, instead require sampling from a conditional distribution p(x0∣y)p(x_0 y). Let ⊆1,…,dO \1,…,d\ index observed coordinates and ℳ=¯M= O the missing ones. The goal is to sample x0ℳ∼p(x0ℳ∣x0).x_0^M p(x_0^M x_0^O). (13) While external-conditioning paradigms such as classifier guidance [5] and classifier-free guidance [12] are effective when y is an external attribute such as a class label or text prompt, they do not apply directly to time series imputation, where the conditioning signal is a subset of the data itself. Diffusion-based imputation methods therefore commonly adopt the following setting [33, 1]. The forward process applies Gaussian corruption only to missing coordinates while preserving observed ones throughout all steps, q(xtℳ|x0ℳ)=(xtℳ,α¯tx0ℳ,(1−α¯t)I),xt=x0∀t.q(x_t^M|x_0^M)=N(x_t^M; α_t\,x_0^M,(1- α_t)I), x_t^O=x_0^O\ ∀ t. (14) The noise predictor receives the partially-noisy tensor x~t x_t (noisy on ℳM, clean on O) together with the missingness mask ℳM and the timestep t, ϵθ(x~t,ℳ,t) _θ( x_t,M,t), giving the network direct access to which coordinates are observed and what their values are. A.1.3 Training objective Training maximizes a variational lower bound (ELBO) on logpθ(x0) p_θ(x_0). The bound decomposes into a sum of KL divergences between the tractable forward posterior q(xt−1∣xt,x0)q(x_t-1 x_t,x_0) and the learned reverse kernel pθ(xt−1∣xt)p_θ(x_t-1 x_t) at each step [10]. Reparameterizing the network as a noise predictor ϵθ _θ and substituting the closed-form expression for xtx_t in Equation (10), the per-step KL reduces (up to a t-dependent reweighting) to the simple denoising loss of 10: ℒsimple(θ)=t∼1,…,T,x0,ϵ∼(0,I)[‖ϵ−ϵθ(α¯tx0+1−α¯tϵ,t)‖22].L_simple(θ)=E_t \1,…,T\,x_0,ε (0,I)[\|ε- _θ( α_tx_0+ 1- α_tε,t)\|_2^2]. (15) For the conditional setting in Appendix A.1.2, the noise predictor ϵθ(x~t,ℳ,t) _θ( x_t,M,t) takes the partially noisy tensor and the missingness mask as input, and the loss is restricted to the missing coordinates [33, 1]: ℒimp(θ)=t∼1,…,T,x0,ϵ∼(0,I),ℳ[‖ϵℳ−ϵθℳ(x~t,ℳ,t)‖22],L_imp(θ)=E_t \1,…,T\,x_0,ε (0,I),M [\,\|ε^M- _θ^M( x_t,M,t)\|_2^2\, ], (16) where the missingness mask ℳM is sampled from a distribution of patterns during training, encouraging the model to generalize across missingness configurations. A.2 Masked Diffusion Language Models Masked Diffusion Language Models (MDLM) [28] adapt the diffusion framework of Appendix A.1 to discrete sequences by replacing Gaussian corruption with an absorbing-state Markov chain [2]. The role of pure noise (0,I)N(0,I) in continuous diffusion is played here by a single distinguished symbol [MASK]: the forward process gradually replaces tokens with [MASK], the prior at t=1t=1 is the all-[MASK] sequence, and the reverse process iteratively unmasks positions back to data tokens. Notation. Let t(i)=iTt(i)= iT and s(i)=i−1Ts(i)= i-1T for i∈1,…,Ti∈\1,…,T\. Tokens take values in 0,1,…,K\0,1,…,K\ where 00 denotes [MASK]. We represent tokens as one-hot vectors in ≔e0,…,eK⊂0,1K+1V \e_0,…,e_K\⊂\0,1\^K+1 and denote the [MASK] one-hot by m≔e0m e_0. For a length-L sequence x=(xl)l=1Lx=(x^l)_l=1^L, let zt=(ztl)l=1Lz_t=(z_t^l)_l=1^L be the latent at time t. We write Cat(⋅,π)Cat(·;π) for a categorical distribution and ⟨⋅,⋅⟩ ·,· for the dot product. A.2.1 MDLM Framework Forward (absorbing mask corruption). The forward process factorizes across token positions and replaces each token by [MASK] with probability 1−αt1- _t: q(zt∣x)=∏l=1LCat(ztl,αtxl+(1−αt)m),q(z_t x)= _l=1^LCat(z_t^l; _tx^l+(1- _t)m), (17) with αt→1 _t→ 1 as t→0t→ 0 and αt→0 _t→ 0 as t→1t→ 1. We parameterize the noise schedule via αt=e−σ(t) _t=e^-σ(t). Our experiments use the log-linear schedule defined by σ(t)=−log(1−t)σ(t)=- (1-t), which yields a linear decay αt=1−t _t=1-t. The defining property of this process is that [MASK] is an absorbing state: once ztl=mz_t^l=m, it remains masked for all later times. Consequently, q(zt|x)q(z_t|x) depends only on which positions have been masked and the reverse process is naturally structured around iterative unmasking. Reverse (learned denoising kernel). MDLM parameterizes the reverse kernel pθ(zs∣zt)=∏l=1Lpθ(zsl∣zt)p_θ(z_s z_t)= _l=1^Lp_θ(z_s^l z_t) via a predicted mixture token xθl(zt)x_θ^l(z_t) over the K non-mask vocabulary entries. Two structural constraints mirror the absorbing forward process: (i) zero-masking, ⟨xθl(zt),m⟩=0 x_θ^l(z_t),m =0, so the predictor never re-introduces [MASK] and (i) carry-over, ztl≠m⇒zsl=ztlz_t^l≠ m z_s^l=z_t^l, so already-unmasked tokens are preserved. As a result, the reverse transition is pθ(zsl∣zt)=Cat(zsl,ztl),ztl≠m,Cat(zsl,αs−αt1−αtxθl(zt)+1−αs1−αtm),ztl=m.p_θ(z_s^l z_t)= casesCat(z_s^l;z_t^l),&z_t^l≠ m,\\ Cat(z_s^l; _s- _t1- _tx_θ^l(z_t)+ 1- _s1- _tm),&z_t^l=m. cases (18) For masked positions, the update is a mixture: with probability αs−αt1−αt _s- _t1- _t the token is unmasked according to xθlx_θ^l, and otherwise it remains [MASK]. A.2.2 Conditional Masking For text infilling, observed tokens can be treated as anchors and excluded from corruption. Let ⊆1,…,LA \1,…,L\ be the anchored index set and ¯ A its complement. Anchored positions are preserved exactly throughout the forward process while non-anchored positions undergo standard absorbing mask corruption: q(zt∣x)=∏l∈[ztl=xl]∏l∈¯Cat(ztl;αtxl+(1−αt)m),q(z_t x)= _l I[z_t^l=x^l] _l∈ ACat\! (z_t^l;\; _tx^l+(1- _t)m ), (19) so that ztl=xlz_t^l=x^l for all t and all l∈l (hard preservation). The reverse kernel mirrors this structure. Anchored positions are copied from ztz_t unchanged, and only non-anchored masked positions are denoised via the standard reverse kernel in Equation (18). This guarantees that observed tokens remain unchanged throughout generation while the model fills in the unobserved ones. This anchored formulation directly matches the imputation setting in Appendix A.1.2. The anchored set A corresponds to the observed coordinates O, and the complement ¯ A corresponds to the missing coordinates ℳM. A.2.3 Training Objective MDLM trains the denoising kernel by minimizing a negative ELBO (NELBO). Owing to the absorbing structure of the forward process, each per-step KL between the forward posterior and the learned reverse kernel reduces in closed form to a weighted cross-entropy over currently masked positions, wt∑l∈ℳtCE(xθl(zt),xl)w_t _l _tCE(x_θ^l(z_t),x^l) where wt=αt−αs1−αtw_t= _t- _s1- _t and ℳt=l∣ztl=mM_t=\l z_t^l=m\ is the set of currently masked positions. The denoiser xθl(zt)x_θ^l(z_t) predicts the clean token at position l directly, paralleling the prediction of clean signal x0x_0 in continuous diffusion. Taking T→∞T→∞, the discrete-time sum converges to a continuous-time objective that samples t∼(0,1)t (0,1) and uses the weight w(t)=αt′1−αtw(t)= _t 1- _t: ℒCT(x,θ)=t[αt′1−αt∑l∈ℳtCE(xθl(zt),xl)],L_CT(x;θ)=E_t [ _t 1- _t _l _tCE\! (x_θ^l(z_t),x^l ) ], (20) which integrates the per-step loss over a continuous time grid and is the form actually used in our experiments. With anchors (Section A.2.2), observed positions are excluded from the forward process and therefore never appear masked at any time t. The training loss is consequently computed only on non-anchor masked positions ℳtanch=l∈¯∣ztl=mM_t^anch=\l∈ A z_t^l=m\, i.e., ℳtM_t is replaced by ℳtanchM_t^anch in Equation (20). The denoiser is trained to reconstruct only the missing coordinates conditional on the observed ones, which is exactly the imputation objective pθ(x0ℳ∣x0)p_θ(x_0^M x_0^O) targeted in Appendix A.1.2. Appendix B Derivation of the Soft-Label Loss Decomposition We provide the derivation of the decomposition stated in Eq. 6. For brevity, we suppress the position index l and the timestep weight w(t)w(t), and consider a single per-token loss ℒ=−⟨,logθ⟩=−∑i=1Ksilogpθ(i),L=- , _θ =- _i=1^Ks_i p_θ(i), (21) where s is the ordinal-aware soft label centered at the ground-truth index y (Eq. 2) and θ∈ΔK−1p_θ∈ ^K-1 is the predicted distribution. Step 1: Decomposition of the soft label. Let ϵ:=1−sy∈[0,1)ε:=1-s_y∈[0,1) denote the total mass that s assigns outside the ground-truth index y. We define :=−syy1−sy,u:= s-s_ye_y1-s_y, (22) where ye_y denotes the one-hot vector at position y. By construction, uy=0u_y=0 and u is a valid probability distribution supported on the ordinal neighborhood i:0<|i−y|≤w\i:0<|i-y|≤ w\. The soft label can then be written as a convex combination =(1−ϵ)y+ϵ.s=(1-ε)\,e_y+ε\,u. (23) Step 2: Splitting the cross-entropy. Substituting Eq. 23 into ℒL and applying the linearity of the inner product gives ℒ =−⟨(1−ϵ)y+ϵ,logθ⟩ =- (1-ε)\,e_y+ε\,u,\ _θ (24) =−(1−ϵ)⟨y,logθ⟩−ϵ⟨,logθ⟩. =-(1-ε)\, _y, _θ -ε\, , _θ . (25) The first term is the standard MDLM cross-entropy against the one-hot target, −⟨y,logθ⟩=−logpθ(y)=ℒNELBO,- _y, _θ =- p_θ(y)=L_NELBO, (26) recovering the per-token NELBO contribution of vanilla MDLM [28]. Step 3: KL form of the second term. The second term is the cross-entropy H(,θ)H(u,p_θ), which by the standard identity decomposes as H(,θ)=−⟨,logθ⟩=DKL(∥θ)+H(),H(u,p_θ)=- , _θ =D_KL(u\,\|\,p_θ)+H(u), (27) where H()=−∑iuiloguiH(u)=- _iu_i u_i is the entropy of u and depends only on s (i.e., independent of θp_θ). Step 4: Combining. Substituting back into Eq. 25 yields ℒ=(1−ϵ)ℒNELBO+ϵDKL(∥θ)+ϵH()⏟const,L=(1-ε)\,L_NELBO+ε\,D_KL(u\,\|\,p_θ)+ ε\,H(u)_const, (28) which is exactly the decomposition in Eq. 6. Two limiting cases are worth noting: • As the soft label sharpens (→ys _y, equivalently ϵ→0ε→ 0), the KL term vanishes and ℒ→ℒNELBOL _NELBO, recovering the vanilla MDLM objective. • For ϵ>0ε>0, the additional non-negative term ϵDKL(∥θ)ε\,D_KL(u\,\|\,p_θ) penalizes probability mass placed outside the ordinal neighborhood of y, since u is supported only on i:0<|i−y|≤w\i:0<|i-y|≤ w\. The same decomposition holds inside the expectation t[w(t)⋅]E_t[w(t)·] in Eq. 5 by linearity, giving the full-objective form ℒdiff=t[w(t)∑l∈ℳt((1−ϵ(l))ℒNELBO(l)+ϵ(l)DKL((l)∥θ(t(l))))]+const.L_diff=E_t\! [w(t) _l _t ((1-ε^(l))\,L_NELBO^(l)+ε^(l)\,D_KL\! (u^(l)\,\|\,p_θ(z_t^(l)) ) ) ]+const. (29) Appendix C Experimental Settings C.1 Datasets Table 9: Statistics of benchmark datasets. Datasets Features Frequency Samples Domain Energy 28 10 min. 19,735 Weather ETTh 7 60 min. 17,420 Temperature Weather 21 10 min. 52,696 Weather Sine 5 - 10000 Simulation We evaluate the proposed MDTIM on four time-series datasets: three real-world benchmarks (Energy, ETTh, Weather) to assess practical validity, and one synthetic Sine dataset for controlled evaluation of periodic dynamics. These datasets cover a broad range of dynamics observed in long-horizon multivariate time-series. The Energy11 1 https://archive.ics.uci.edu/dataset/374/appliances+energy+prediction dataset contains electricity-related measurements with strong seasonality and occasional irregular variations. The ETTh22 2 https://github.com/zhouhaoyi/ETDataset dataset is a standard benchmark derived from Electricity Transformer Dataset [43] and is commonly used to evaluate long-range dependency modeling. The Weather33 3 https://w.bgc-jena.mpg.de/wetter/ [37] dataset includes multivariate meteorological observations characterized by nonlinear interactions, seasonal trends, and stochastic fluctuations. In addition, the synthetic Sine dataset offers a controlled environment with known periodic structure, which helps validate fundamental forecasting behavior under minimal uncontrolled factors. Overall, these datasets were selected to ensure domain diversity and to evaluate performance under both structured periodic signals and noisy, non-stationary real-world patterns. Dataset specifications are summarized in Table 9. For the implementation and training of the baseline models, we adopted the codebase curated by 744 4 https://github.com/WenjieDu/Awesome_Imputation. C.2 Training Configurations We provide the detailed hyperparameter settings used to train MDTIM. To ensure reproducibility, we list the common architectural parameters and training schemes applied across all datasets. All experiments are conducted with dual Intel Xeon Gold 6444Y CPUs and a single NVIDIA H100 PCIe GPU (80GB). C.2.1 Hyperparameters The architecture of MDTIM is based on the Factorized Temporal-Variate Backbone described in Section 4. Table 10 summarizes detailed hyperparameters. We utilized consistent hyperparameter settings across all datasets to ensure reproducibility. Table 10: Hyperparameters for MDTIM architecture and training. Category Parameter Value Description Model Hidden Size (D) 256 Dimension of hidden states Attention Heads 16 Number of heads in MSA DiT Blocks (L) 5 Number of factorized layers Dropout 0.2 Dropout probability Bins (K) 40 Vocabulary size for discretization Output Range 1.5 Value range [−1.5,1.5][-1.5,1.5] Cond. Dim 16 Time-step embedding dimension Training Batch Size 256 - Learning Rate 3×10−43× 10^-4 - Scheduler Const. w/ Warmup - Warmup Steps 2,500 - Max Steps 10,000 - Gradient Clip 1.0 Norm clipping value EMA Decay 0.995 Exponential moving average FFT Weight (λ) 1.0 Spectral loss weight C.2.2 Noise Schedule and Time-Dependent Weighting We employ a Log-Linear Noise Schedule to compute the continuous-time diffusion process. For a time step t∈[0,1]t∈[0,1] and a small constant ϵ=10−3ε=10^-3, the total noise σ(t)σ(t) and its rate of change are defined as: σ(t)=−log(1−(1−ϵ)t),dσdt=1−ϵ1−(1−ϵ)tσ(t)=- (1-(1-ε)t ), dσdt= 1-ε1-(1-ε)t (30) To ensure balanced training across varying noise levels, we apply importance sampling via a time-dependent loss weight w(t)w(t). This weight effectively normalizes the contribution of each time step to the objective function: w(t)=dσ(t)/dteσ(t)−1w(t)= dσ(t)/dte^σ(t)-1 (31) This weighting scheme assigns higher importance to low-noise regions (t≈0t≈ 0), prioritizing the learning of fine-grained details from clean data, while down-weighting highly corrupted states (t≈1t≈ 1) where reconstruction is ambiguous. C.3 Evaluation Metrics We evaluate both point imputation accuracy and the quality of predictive distributions. Let yi,t,cy_i,t,c and y^i,t,c y_i,t,c denote the ground truth and the imputed value for sample i, time step t, and channel c. All metrics are computed exclusively on the masked positions, denoted by the set of indices ℳM. MAE and MSE. We report Mean Absolute Error (MAE) and Mean Squared Error (MSE), averaged over the masked entries ℳM: MAE=1|ℳ|∑(i,t,c)∈ℳ|y^i,t,c−yi,t,c|,MSE=1|ℳ|∑(i,t,c)∈ℳ(y^i,t,c−yi,t,c)2.MAE= 1|M| _(i,t,c) | y_i,t,c-y_i,t,c |, = 1|M| _(i,t,c) ( y_i,t,c-y_i,t,c )^2. (32) In probabilistic settings, we set y^i,t,c y_i,t,c to the sample mean of M generated trajectories, y^i,t,c=1M∑m=1My~i,t,c(m) y_i,t,c= 1M _m=1^M y_i,t,c^(m). CRPS. To assess the probabilistic imputation performance, we use the Continuous Ranked Probability Score (CRPS) [21]. For a predictive CDF Fi,t,cF_i,t,c at a masked position, CRPS is defined as CRPS(Fi,t,c,yi,t,c)=∫−∞(Fi,t,c(z)−[z≥yi,t,c])2dz.CRPS\! (F_i,t,c,y_i,t,c )= _-∞^∞ (F_i,t,c(z)-I[z≥ y_i,t,c] )^2\,dz. (33) Lower CRPS indicates better probabilistic imputation, rewarding both calibration and sharpness. Sample-based CRPS Estimation. Since the model produces M stochastic samples y~i,t,c(m)m=1M∼Fi,t,c\ y_i,t,c^(m)\_m=1^M F_i,t,c, we approximate CRPS by CRPS¯i,t,c=1M∑m=1M|y~i,t,c(m)−yi,t,c|−12M2∑m=1M∑n=1M|y~i,t,c(m)−y~i,t,c(n)|. CRPS_i,t,c= 1M _m=1^M | y_i,t,c^(m)-y_i,t,c |- 12M^2 _m=1^M _n=1^M | y_i,t,c^(m)- y_i,t,c^(n) |. (34) We report the average of CRPS¯i,t,c CRPS_i,t,c over all masked entries in ℳM. Appendix D Additional Experimental Results In this section, we provide supplementary experimental results that complement the main paper. Section D.1 reports the full quantitative comparison including the 50% missing ratio omitted from Table 1 for space. Section D.2 presents an ablation comparing our discrete formulation against a continuous-modeling counterpart, isolating the contribution of Stochastic Discretization. D.1 Full Imputation Performance Across All Missing Ratios Table 11 reports the complete results of our main imputation experiment, including the 30%, 50%, and 70% missing ratios under both Uniform and Geometric scenarios. The 50% column extends the trend reported in the main paper: MDTIM consistently maintains the lowest MAE across the majority of settings, and the relative improvement over baselines remains stable across missing ratios. This confirms that the performance advantage of MDTIM is not specific to any particular corruption level, but holds robustly throughout the regime of partial observation. Table 11: Quantitative comparison of multivariate time-series imputation performance (L=48L=48). We report the MAE of MDTIM and baselines averaged over 3 random seeds. The best results are highlighted in bold. Dataset Energy ETTh Weather Sine Missing Type Uniform Geometric Uniform Geometric Uniform Geometric Uniform Geometric Model 30% 50% 70% 30% 50% 70% 30% 50% 70% 30% 50% 70% 30% 50% 70% 30% 50% 70% 30% 50% 70% 30% 50% 70% RNN BRITS 0.246 0.296 0.379 0.329 0.344 0.366 0.194 0.235 0.299 0.227 0.258 0.292 0.050 0.057 0.073 0.058 0.063 0.071 0.010 0.013 0.021 0.020 0.019 0.019 MRNN 1.086 1.107 1.143 1.087 1.108 1.143 0.743 0.762 0.782 0.751 0.765 0.781 0.651 0.653 0.661 0.654 0.656 0.659 0.169 0.169 0.170 0.169 0.169 0.170 GRUD 0.364 0.387 0.457 0.426 0.425 0.445 0.310 0.336 0.394 0.348 0.362 0.383 0.104 0.176 0.369 0.164 0.231 0.350 0.008 0.008 0.014 0.013 0.012 0.012 Transformer Transformer 0.323 0.335 0.401 0.354 0.356 0.393 0.168 0.201 0.252 0.183 0.212 0.247 0.076 0.067 0.070 0.083 0.072 0.068 0.039 0.041 0.046 0.042 0.044 0.045 Informer 0.344 0.338 0.379 0.371 0.357 0.373 0.207 0.240 0.302 0.226 0.256 0.296 0.047 0.046 0.059 0.052 0.050 0.058 0.027 0.025 0.039 0.036 0.032 0.037 PatchTST 0.586 0.336 0.529 0.528 0.378 0.517 0.202 0.218 0.272 0.231 0.243 0.262 0.074 0.061 0.082 0.081 0.069 0.077 0.015 0.011 0.017 0.019 0.014 0.015 SAITS 0.177 0.185 0.212 0.204 0.200 0.208 0.140 0.166 0.210 0.152 0.176 0.206 0.045 0.044 0.049 0.050 0.048 0.048 0.026 0.023 0.026 0.030 0.026 0.025 Imputeformer 0.066 0.100 0.219 0.113 0.143 0.194 0.146 0.178 0.245 0.165 0.196 0.235 0.049 0.068 0.171 0.090 0.111 0.149 0.006 0.012 0.041 0.012 0.019 0.037 CNN TimesNet 0.617 0.683 0.818 0.689 0.726 0.808 0.593 0.653 0.719 0.628 0.666 0.717 0.211 0.202 0.351 0.225 0.223 0.348 0.167 0.188 0.220 0.189 0.200 0.217 SCINet 0.557 0.418 0.571 0.546 0.442 0.565 0.238 0.252 0.323 0.261 0.273 0.317 0.074 0.068 0.089 0.086 0.076 0.087 0.019 0.018 0.031 0.026 0.024 0.029 Linear DLinear 0.795 0.351 0.527 0.684 0.388 0.517 0.379 0.319 0.445 0.373 0.357 0.435 0.372 0.149 0.223 0.345 0.161 0.220 0.070 0.032 0.053 0.062 0.041 0.051 FiLM 0.877 0.380 0.520 0.791 0.413 0.512 0.696 0.608 0.627 0.707 0.627 0.623 0.380 0.156 0.209 0.360 0.165 0.206 0.127 0.102 0.109 0.132 0.110 0.108 FreTS 0.157 0.147 0.219 0.225 0.184 0.206 0.222 0.238 0.299 0.262 0.266 0.286 0.071 0.059 0.080 0.080 0.067 0.074 0.068 0.037 0.068 0.075 0.049 0.065 Generative GPVAE 0.476 0.509 0.730 0.501 0.532 0.727 0.333 0.349 0.449 0.369 0.385 0.440 0.158 0.153 0.255 0.172 0.170 0.252 0.152 0.150 0.162 0.153 0.151 0.162 USGAN 0.264 0.323 0.412 0.330 0.362 0.401 0.208 0.243 0.302 0.240 0.266 0.295 0.087 0.094 0.122 0.101 0.104 0.118 0.014 0.018 0.029 0.026 0.026 0.026 CSDI 0.094 0.112 0.139 0.110 0.122 0.136 0.160 0.192 0.244 0.177 0.206 0.240 0.039 0.043 0.049 0.044 0.046 0.048 0.003 0.004 0.004 0.004 0.004 0.004 FGTI 0.050 0.071 0.100 0.065 0.080 0.097 0.218 0.255 0.349 0.293 0.311 0.328 0.038 0.043 0.051 0.046 0.047 0.049 0.001 0.002 0.003 0.002 0.003 0.003 MDTIM (Ours) 0.044 0.061 0.085 0.053 0.067 0.082 0.127 0.157 0.211 0.146 0.173 0.205 0.032 0.036 0.044 0.036 0.039 0.043 0.002 0.003 0.004 0.003 0.003 0.003 D.2 Comparison with Continuous Modeling To isolate the contribution of our discrete formulation, we compare MDTIM (Disc), the proposed model, against a continuous-modeling variant MDTIM (Cont). MDTIM (Cont) shares the same Factorized Temporal-Variate backbone and masked diffusion training paradigm, but operates directly on continuous input values and is optimized with a Mean Squared Error (MSE) loss. In contrast, MDTIM (Disc) employs our Stochastic Discretization to map continuous signals into discrete tokens and optimizes the ordinal-aware discrete diffusion objective described in Section 4. (a) Uniform Missing (b) Geometric Missing Figure 4: Comparison between MDTIM (Disc) and MDTIM (Cont) on the Energy dataset, alongside representative baselines (SAITS, CSDI). Results are reported under both uniform and geometric missing scenarios across varying missing rates. As shown in Figure 4, MDTIM (Cont) already performs competitively with strong baselines, indicating that the masked diffusion paradigm of learning to reconstruct from varying corruption levels with timestep conditioning is itself effective even in continuous spaces. Nevertheless, MDTIM (Disc) consistently achieves lower MAE across all missing rates and missing types. This performance gap confirms that the explicit structural separation between valid observations and missing placeholders, achieved through our orthogonal tokenization, provides a benefit beyond the masked diffusion training scheme alone, and is crucial for time series imputation. Appendix E Visualization of Imputation Results We provide comprehensive visualizations of the imputation results across all benchmark datasets: ETTh, Energy, Weather, and Sine with MDTIM. Figures 5, 6, 7, and 8 illustrate the reconstructed time series for all channels (ETTh, Sine) or subset of channels (Energy, Weather). In each figure, the left column displays the results under 50% uniform masking, while the right column shows under 50% geometric masking. The blue lines represent the ground truth values, while the green lines denote the imputed values (x x) reconstructed by MDTIM. The green shaded areas indicate the estimated uncertainty intervals. As observed, MDTIM effectively captures the complex temporal dynamics and periodicity of the multivariate time series. (a) Uniform Missing (b) Geometric Missing Figure 5: Visualization of imputation results on the ETTh dataset. (a) Uniform Missing (b) Geometric Missing Figure 6: Visualization of imputation results on the Energy dataset. (a) Uniform Missing (b) Geometric Missing Figure 7: Visualization of imputation results on the Weather dataset. (a) Uniform Missing (b) Geometric Missing Figure 8: Visualization of imputation results on the Sine dataset.