Multi-Scale Forecasting of Natural Rubber Prices Using a VMD-Hybrid BiLSTM-Transformer Framework

1. Introduction

Natural rubber is a critical industrial input and a strategic agricultural commodity that plays a central role in global manufacturing supply chains. It is widely used in the production of automobile tires, medical gloves, engineering components, and a broad range of technical rubber products. Global natural rubber production is highly concentrated in Southeast Asia, particularly in Thailand, Indonesia, and Malaysia, which together account for the majority of world supply [1,2]. Among these producers, Thailand has consistently remained the world’s leading exporter of natural rubber. In 2025, Thailand produced approximately 4.79 million tonnes, accounting for 32.18% of global production, with about 85.79% of output exported to international markets [3]. Price volatility in natural rubber markets therefore carries direct consequences for the income of millions of smallholder farmers and the cost structure of downstream industries that rely heavily on rubber-based inputs.

Natural rubber prices are jointly driven by agricultural seasonality, macroeconomic conditions, exchange rate dynamics, and speculative activity in futures markets, producing price dynamics that are simultaneously nonlinear, non-stationary, and subject to structural breaks [4,5,6]. These forces generate significant uncertainty for stakeholders across the rubber value chain. In international trade, Ribbed Smoked Sheet No.3 (RSS3) serves as the benchmark grade for natural rubber transactions, with export prices typically quoted on a Free on Board (FOB) basis at major Thai ports. Price co-movement across SHFE, SGX, and JPX (TOCOM) futures markets creates cross-market information flows that provide additional predictive content beyond domestic spot prices alone, motivating the inclusion of exchange-traded settlement prices as model inputs [7]. This cross-market structure is not symmetric: [8] demonstrate that SHFE Granger-causes TOCOM but not vice versa, while Kepulaje et al. [9] [pp. 159–161] show that SGX RSS3 Rubber Futures (SRU) and SGX TSR20 Rubber Futures (STF) transmit more consistently into regional spot prices than TOCOM/JPX Rubber Futures (RSS3)(JRU) — implying that the futures exchanges operate as a layered rather than a flat information system, and that the informational weight of each venue depends on whether one is studying financial benchmarking or physical export pricing.

Accurate forecasting of commodity prices therefore remains a central issue in both economic research and industrial decision-making. Commodity price forecasting has long relied on econometric models such as ARIMA, VAR, and GARCH, which provide well-established frameworks for capturing temporal dependencies and volatility structures [10,11]. In the natural rubber literature specifically, this approach is exemplified by [12]. (2018), who apply the Box–Jenkins procedure to Malaysian SMR20 and identify ARIMA(1,1,0) as the best-fitting model, and by [13], who show that a simultaneous supply–demand–price system outperforms univariate ARIMA across RMSE, MAE, and information criteria — demonstrating that purely univariate extrapolation is insufficient once commodity market structure is made explicit. However, these approaches rest on linearity assumptions that limit their capacity to represent the complex nonlinear dynamics observed in modern commodity markets, and their performance often deteriorates further when the underlying data exhibit structural breaks and high non-stationary characteristics that are particularly pronounced in energy and agricultural commodity price series [14,15]. To address these shortcomings, subsequent studies introduced machine learning methods including support vector machines (SVM), multi-layer perceptrons (MLP), random forests, extreme gradient boosting (XGBoost), and back-propagation neural networks (BPNN), which are better suited for capturing nonlinear relationships in complex datasets [16,17,18]. Hybrid frameworks combining statistical and machine learning components—such as ARIMA–SVM—have also been proposed to integrate linear structure with nonlinear residual dynamics [19,20]. Nevertheless, these approaches are limited in their ability to capture complex nonlinear dynamics and long-range dependencies, motivating the adoption of deep learning alternatives.

Building on this shift, recurrent architectures—particularly Long Short-Term Memory (LSTM) and Gated Recurrent Unit (GRU) networks—have been widely adopted for their ability to model long-term temporal dependencies more effectively than shallow learning algorithms, demonstrating superior performance over classical statistical benchmarks in energy demand prediction, financial market forecasting, and environmental time-series modeling [21,22]. For natural rubber specifically, [23] demonstrate that a multi-layer LSTM achieves 95.88% accuracy on Thai RSS3 prices, while [24] confirm that LSTM outperforms ARIMA on Malaysian SMR20 across MAE, RMSE, and MAPE — with both studies noting that predictive gains remain incomplete without explicit incorporation of exogenous drivers such as exchange rates, energy prices, and futures markets. However, recurrent models can struggle with highly volatile and non-stationary series, particularly when the data exhibit abrupt structural changes or mixed-frequency dynamics [25]. Bidirectional Long Short-Term Memory (BiLSTM) networks extend the conventional LSTM architecture by processing sequential data simultaneously in the forward and backward directions, enabling the model to integrate contextual information from both past and future states within the observation window. This bidirectional structure yields richer temporal representations than unidirectional recurrent networks and has demonstrated consistent advantages over conventional LSTM and classical statistical benchmarks in environments characterized by strong seasonality, nonlinear interactions, and abrupt fluctuations—including renewable energy forecasting and environmental time-series prediction [21,22,26].

To further extend the representational capacity of BiLSTM, several lines of re-search have proposed hybrid architectures. CNN–BiLSTM models combine convolution layers for local feature extraction with BiLSTM for bidirectional temporal learning, significantly improving prediction performance in energy and meteorological forecasting applications [27,28]. Attention mechanisms have also been integrated to allow the network to emphasize the most informative time steps, enhancing performance on non-stationary sequences and series with irregular fluctuations [29]. A separate strand of research applies meta-heuristic optimization to tune BiLSTM hyper-parameters, and combines BiLSTM with gradient boosting or ensemble methods to capture residual nonlinear patterns, yielding further gains over single-model recurrent baselines in energy and meteorological forecasting applications. Despite these advances, BiLSTM-based hybrids remain limited in their ability to separate and individually model components operating across distinctly different temporal scales, a constraint that becomes particularly consequential for commodity price series with coexisting trend, cyclical, and noise dynamics.

A parallel research direction addresses the recurrent model’s limited ability to disentangle components at different temporal scales by incorporating signal decomposition as a preprocessing stage. Methods including Empirical Mode Decomposition (EMD), Variational Mode Decomposition (VMD), Wavelet Transform (WT), and Seasonal-Trend Decomposition using LOESS (STL) separate the original series into components representing distinct frequency structures—trend, seasonal, and residual—thereby reducing non-stationary and simplifying the learning task for downstream models [7,29,30]. A substantial body of empirical evidence confirms that decomposition-based hybrid frameworks, including EMD–LSTM, VMD–LSTM, VMD–GRU, VMD–TCN–LSTM, and STL–LSTM—significantly improve forecasting accuracy across coal markets, solar power generation, wind forecasting, and electricity load prediction [29]. Recent evidence from high-frequency financial forecasting likewise shows that VMD can separate trend, periodic, and disturbance components in non-stationary series and that decomposition-enhanced recurrent models can outperform standalone econometric and deep-learning benchmarks [25]. In commodity-price forecasting, VMD-based hybrid models have similarly been shown to outperform both single-model and earlier decomposition-based alternatives, with further gains when structurally relevant long-run drivers are incorporated into the forecasting system [31]. Recent evidence from environmental forecasting further confirms that coupling VMD with sequence models can substantially improve predictive accuracy by separating short- and long-run signal structures before learning, relative to models trained on the raw series alone [32]. Related evidence from drought forecasting likewise demonstrates that VMD-based preprocessing can improve short-lead predictive accuracy relative to standalone machine learning models, reinforcing the value of decomposition in nonlinear forecasting tasks [33]. In related commodity-market research, VMD-based recurrent models have also been shown to improve forecasting by decomposing volatile price series into simpler sub-components before learning, confirming the value of decomposition for nonlinear price dynamics [34].

Building on these decomposition gains, a further line of research demonstrates that combining VMD with attention-based sequence architectures yields complementary improvements beyond what either component achieves alone. Related forecasting studies have shown that VMD can improve downstream deep-learning performance by decomposing complex time series into more learnable components, while attention-based LSTM encoder-decoder structures further enhance the selection of in-formative features and time steps [35] [pp. 922–923]. Evidence from short-term wind power forecasting further suggests that VMD is especially effective when embedded within a broader hybrid pipeline rather than used as an isolated preprocessing step, improving both stability and predictive accuracy in volatile series [36]. Recent FX forecasting evidence additionally indicates that VMD-enhanced bidirectional recurrent architectures, when paired with systematic hyper-parameter tuning and regularization, can materially improve out-of-sample performance in low-frequency financial time series [37].

Within the natural rubber literature specifically, evidence from rubber futures forecasting shows that VMD can extract economically meaningful multi-scale components, while bidirectional recurrent architectures improve both fitting performance and directional prediction; importantly, the predictive value of individual modes depends on their time-scale correspondence with the forecast target [7] [p. 6]. This interpretation is further supported by energy-price forecasting evidence showing that decomposed sub-sequences differ materially in regularity and information content, and that forecasting performance improves when model design reflects the heterogeneous nature of high- and low-frequency components [30].

More recently, Transformer-based architectures have emerged as a compelling complement to recurrent models, leveraging self-attention mechanisms to capture long-range temporal dependencies across entire sequences in parallel—a capability that recurrent networks approximate only sequentially [38]. Influential variants including Informer, Autoformer, FEDformer, PatchTST, and ReVIN have progressively introduced sparse attention, decomposition-based encoders, frequency-enhanced representations, patch-level temporal embedding, and distribution-shift normalization, substantially improving performance on long-sequence forecasting tasks [35,39,40]. The complementary strengths of recurrent and attention-based models have motivated a further generation of hybrid architectures—including LSTM–Transformer hybrids, wavelet–Transformer–LightGBM frameworks, and dual-branch designs integrating short-term sequential learning with multi-scale decomposition—each seeking to capture both local temporal structure and global dependencies within a unified frame-work [41,42,43,44].

Despite these advances, most existing decomposition-based models treat signal separation purely as a preprocessing step to improve input quality, without redesigning the forecasting architecture to exploit the structural differences between decomposed components. Indeed, recent evidence from short-term wind power forecasting suggests that VMD is most effective when embedded within a broader hybrid pipeline rather than applied as an isolated preprocessing step, improving both stability and predictive accuracy in volatile series [36]. Recent hybrid forecasting research similarly argues that direct end-to-end modeling is often inadequate for highly nonlinear and volatile series, and that decomposition-based feature extraction combined with bidirectional deep learning and attention improves the representation of multi-scale temporal structure [45]. In the majority of implementations, each intrinsic mode function is forecast independently by a separate model instance, and final predictions are obtained by summation — an approach that, as demonstrated in this study, is prone to variance collapse when applied to stationary (differenced) commodity price series. This failure mode is not incidental: Shumailov et al. [46] [pp. 755–757] demonstrate theoretically that compound estimation pipelines systematically erode distributional variance across stages, while Rajpal et al. [47] [pp. 13–15] confirm its output-level analogue in financial forecasting, showing that symmetric loss objectives push models toward degenerate, one-sided prediction regimes whose apparent performance masks near-zero predictive content. Recent oil-price forecasting studies have moved beyond single-stage decomposition toward secondary decomposition-reconstruction-ensemble frameworks, arguing that the information embedded in high-frequency components is too rich to be exhausted by a single decomposition pass [48]. This architectural limitation and the need for a more parsimonious yet information-complete alternative represents the central research gap that the present work seeks to address. This gap is particularly consequential in the natural rubber context, where Fakthong et al. [49] argue that existing time-series and machine-learning studies are heavily oriented toward short-term fluctuation tracking and pay limited attention to the combined structural effects of economic, environmental, and policy variables — reinforcing the need for architectures that preserve multi-scale information rather than collapsing it through sequential single-model estimation.

To address these gaps, this study proposes the VMD–Hybrid BiLSTM–Transformer, a dual-pathway forecasting framework for daily RSS3 FOB natural rubber price changes using a 24-feature economic input set over the period 2018–2026. The principal contributions of this study are threefold. First, rather than predicting each VMD intrinsic mode function independently and aggregating forecasts — an approach demonstrated here to produce variance collapse on differenced commodity price series (StdR $=0.04$–$0.15$) — all five IMF series are appended directly to the economic feature matrix, preserving multi-scale frequency information within a single forward pass. Second, a cross-stage feature availability analysis is conducted to establish the minimum sufficient input configuration for reliable forecasting, providing a reproducible validation framework for practitioners operating under data constraints. Third, the limitations of directional accuracy as a standalone evaluation metric are formally demonstrated through confusion matrix and variance ratio diagnostics, and a suite of complementary metrics — including Pearson correlation and Standard Deviation Ratio — is proposed for evaluating differenced commodity price models. The theoretical basis for this metric combination is established by Taylor [50] [pp. 7183–7184], who demonstrates that correlation and variance ratio are not substitutes but complementary descriptors: correlation measures co-movement in shape, while the standard deviation ratio measures amplitude fidelity, and a reduction in RMS error cannot be taken as evidence of improved skill if forecast variance has been suppressed in the process. This evaluation critique aligns with a broader methodological literature demonstrating that directional accuracy is structurally vulnerable to benchmark definition and class imbalance: Bürgi [51] shows that DA is not a self-sufficient concept without knowledge of magnitude sensitivity and user objective, and McCarthy and Snudden [52] [pp. 1–7] demonstrate that temporal aggregation can mechanically inflate success ratios beyond 0.5 even under a random walk, generating pseudo-skill where no genuine predictability exists. Beyond benchmark sensitivity, Costantini et al. [53] show that DA becomes materially informative only when paired with magnitude-sensitive directional value measures, and Costantini and Kunst [54] demonstrate that embedding DA directly into model-selection criteria seldom yields robust gains and can worsen MSE performance — implying that DA should function as a supplementary diagnostic rather than a primary evaluation objective. The remainder of this paper is organized as follows: Section 2 describes the data, preprocessing pipeline, decomposition procedure, model architecture, and evaluation protocol; Section 3 reports the empirical results; Section 4 discusses the economic interpretation and practical implications; and Section 5 concludes with limitations and directions for future research.

2. Materials and Methods

This section describes the data, preprocessing pipeline, and architecture of the proposed VMD-Hybrid BiLSTM-Transformer forecasting framework. Section 2.1 motivates the primary analytical sample through a cross-stage comparison. Section 2.2–2.5 describe, respectively, data and feature construction, model architecture, training configuration, and evaluation protocol.

2.1. Cross-Stage Comparison and Rationale for the Primary Analytical Sample

The full datasets (2 October 1999 – 27 February 2026, n = 7,492 daily observations) is structured into three chronological stages defined by input feature availability. Stage 1 (21 April 2003 – 30 December 2014, n = 3,030) includes 11 features comprising domestic rubber spot prices, Cup-Lump, USS, and exchange-traded futures (JPX, SHFE), covering a period before foreign exchange, energy, and macroeconomic series became available. Stage 2 (1 January 2015 – 4 May 2018, n = 872) extends coverage to 22 features following the availability of USD/THB, CNY/THB, USD/CNY, Brent, WTI, China PMI Manufacturing, Baltic Dry Index, and ENSO ONI. Stage 3 (7 May 2018 – 27 February 2026, n = 2,675) represents the first period with complete coverage of all 24 input features, including SGX RSS3 and TSR20 settlement prices, and constitutes the primary analytical sample for this study. To assess the sensitivity of model performance to input feature availability, the VMD-Hybrid BiLSTM-Transformer was trained and evaluated independently on all three sub-samples using only the features available within each period. All models were evaluated on their respective held-out test partitions using an identical protocol.

Prior to model training, VMD decomposition was applied independently to each sub-sample to determine the optimal number of modes K. Stage 1 yields a highly concentrated energy distribution under K = 5, with IMF1 accounting for 95.3% of total signal variance—indicating that the pre-2015 rubber market was dominated by a single low-frequency trend, consistent with the sustained commodity super-cycle of 2003–2011 followed by a sharp price correction. Stage 2 requires K = 7 to adequately capture signal complexity, with energy distributed more evenly across seven components (range: 7.6%–23.9%), reflecting the heightened multi-frequency dynamics of the 2015–2018 transition period. Stage 3 under K = 5 presents a more structured decomposition: IMF1 retains dominance at 78.6%, while secondary components carry economically interpretable energy shares (IMF2 = 8.0%; IMF5 = 6.7%); full decomposition results for Stage 3 are reported in Table 9. This cross-stage analysis demonstrates that the frequency structure of RSS3 prices is non-stationary across market regimes, and that K = 5 is the appropriate decomposition depth for Stage 3.

Out-of-sample performance results are summarized in Table 1. Stage 3 substantially outperforms the reduced-feature configurations on all primary metrics. The most discriminating indicator is the Pearson correlation between predicted and realized price changes: Stage 3 achieves r = 0.82, compared with r = 0.38 for Stage 1 and r = 0.12 for Stage 2 evaluated on the Stage 3 test window. The Standard Deviation Ratio (StdR)—which equals 1.0 under ideal forecast dispersion—further reveals that Stage 1 (StdR = 0.34) and Stage 2 (StdR = 0.11) both exhibit pronounced variance collapse: the model defaults toward near-zero predictions, capturing directional frequency but suppressing magnitude entirely. Stage 3 achieves StdR = 0.83, confirming that the full feature set enables the model to reproduce the dispersion of the target series. Directional accuracy is comparably high across all stages (96.8%–98.6%), consistent with the known limitation of DA as a discriminating metric on difference commodity series where zero-change days are excluded. This uniform DA pattern illustrates precisely the benchmark-sensitivity problem documented by [51] [pp. 7909–7910]: when class composition is skewed and the underlying series is highly persistent, DA can remain nominally high across radically different model configurations while conveying no information about whether the model captures conditional dynamics — a limitation that StdR and Pearson r directly address.

These results provide two complementary pieces of evidence supporting the study design. First, the stepwise degradation in StdR and Pearson r as features are removed confirms that the SGX settlement prices introduced in Stage 3 carry incremental pre-dictive content that is not redundant with the 22 features available in Stage 2. Second, the inability of the Stage 1 and Stage 2 models to reproduce forecast variance—despite achieving high directional accuracy—illustrates a well-known failure mode of deep learning models on stationary differenced series when input information is insufficient: the model learns the unconditional mean rather than the conditional distribution. This pattern is consistent with the theoretical account of variance collapse developed by Shumailov et al. [46] [pp. 755–757], who show that models operating under insufficient informational input systematically converge toward overconcentrated, low-variance outputs — a process that is provably inevitable under finite sampling and compound estimation error. Stage 3 is thus the minimum sufficient feature set for this forecasting task, and all primary results reported in Section 3 are based on this sub-sample.

2.2. Data, Preprocessing, and Input Features

2.2.1. Target Variable and Stationary

The study focuses on the Stage 3 sub-sample (7 May 2018 – 27 February 2026, n = 2,675 observations), which represents the first period with complete coverage of all 24 input features, including SGX settlement prices. The target variable is the first-differenced RSS3 FOB near-month price (Baht/kg/day):

$$\Delta p_t=p_t-p_{t-1}$$

(1)

where $p_t$ denotes the RSS3 FOB spot price on day t (Baht/kg). Augmented Dickey–Fuller (ADF) unit root tests confirmed that 21 of the 24 economic features are non-stationary in levels; full ADF results are reported in Appendix B. All non-stationary series were transformed to first-differenced prior to model input. Descriptive statistics of the differenced target series over the test partition are provided in Table 2. The series exhibits excess kurtosis and mild negative skewness, consistent with stylized facts of commodity price changes.

2.2.2. Data Partitioning

The Stage 3 sub-sample is partitioned chronologically into three non-overlapping sets. The test start date is pinned at 18 September 2025 to ensure comparability with the primary benchmark evaluation (n = 237). Chronological ordering is strictly preserved to prevent look-ahead bias.

Table 3. Chronological data split — Stage 3.

Partition	Date range	Observations	Share	Purpose
Training	2018-05-07 → 2023-08-27	2,140	80%	Parameter estimation
Validation	2023-08-28 → 2025-09-17	267	10%	Early stopping / LR scheduling
Test	2025-09-18 → 2026-02-27	237	∼9%	Out-of-sample evaluation

Table 3. Chronological data split — Stage 3.

Partition	Date range	Observations	Share	Purpose
Training	2018-05-07 → 2023-08-27	2,140	80%	Parameter estimation
Validation	2023-08-28 → 2025-09-17	267	10%	Early stopping / LR scheduling
Test	2025-09-18 → 2026-02-27	237	∼9%	Out-of-sample evaluation

2.2.3. Missing Value Handling and Normalization

Missing values arising from weekends and public holidays are handled by for-ward-filling (maximum three consecutive days), with remaining gaps filled by linear interpolation (maximum ten days). All features are normalized using a MinMax scaler fitted exclusively on the training partition to prevent data leakage:

$$\tilde{x}=2·{\displaystyle \frac{\left(x-x_{min}\right)}{\left(x_{max}-x_{min}\right)}}-1$$

(2)

$$\Delta {\widehat{p}}_t={\displaystyle \frac{{\widehat{y}}_t+1}{2}}·\left(\Delta p_{max}-\Delta p_{min}\right)+\Delta p_{min}$$

(3)

2.2.4. Input Feature Matrix

The input feature matrix comprises 24 economic variables spanning five categories, as summarized in Table 4. These features capture the main demand-side and supply-side drivers of natural rubber prices, including futures markets, foreign ex-change rates, energy prices, and macroeconomic indicators. COVID dummy and ENSO ONI are retained in levels; all other features are first-differenced. This variable selection is grounded in the structural rubber forecasting literature: [55] [pp. 8–9] identify exchange rates, synthetic rubber prices, GDP growth, crude oil prices, and world stocks as central price determinants in a simultaneous supply–demand framework, while [56] [pp. 1471–1472] confirm through logic-mining variable selection that exports, imports, stocks, exchange rates, and crude oil prices carry the highest predictive weight — providing empirical justification for the multi-category feature design adopted here.

Figure 1. RSS3 FOBm1 daily price level and first-differenced series (Stage 3, 2018–2026). Note: Upper panel: daily RSS3 FOB price level (Baht/kg). Shaded regions denote training ($n=2,140$), validation ($n=267$), and test ($n=237$) partitions. Lower panel: first-differenced series $\Delta p_t$ (Baht/kg/day) with zero reference line. ADF test confirms stationary of the differenced series ($t=-18.73$, $p<0.001$).

Figure 1. RSS3 FOBm1 daily price level and first-differenced series (Stage 3, 2018–2026). Note: Upper panel: daily RSS3 FOB price level (Baht/kg). Shaded regions denote training ($n=2,140$), validation ($n=267$), and test ($n=237$) partitions. Lower panel: first-differenced series $\Delta p_t$ (Baht/kg/day) with zero reference line. ADF test confirms stationary of the differenced series ($t=-18.73$, $p<0.001$).

Features are grouped into five clusters: rubber spot and futures (blue), exchange futures (purple), exchange rates (green), energy (orange), and macro/demand indicators (brown). Within-cluster correlations are substantially higher than cross-cluster correlations, confirming that each feature group captures partially independent in-formation. RSS3 FOBm1 (target variable) shows moderate positive correlation with other rubber futures but near-zero correlation with energy and macro variables in the differenced space.The inclusion of SGX RSS3 and TSR20 settlement prices alongside SHFE and JPX futures is motivated by the transmission evidence in [9] [pp. 159–161], who find that SGX contracts carry greater short-run informational weight in the futures-to-spot transmission equation than TOCOM/JPX Rubber Futures (USS rubber)— and by Pinitjitsamut [57], who shows that FOB Bangkok serves as the operative export-pricing node through which SHFE signals are filtered into domestic Thai prices, making both the exchange and export pricing layers necessary inputs for a complete forecasting system.

Figure 2. Pearson correlation heatmap of 24 input features (post-differencing, training partition).

Figure 2. Pearson correlation heatmap of 24 input features (post-differencing, training partition).

2.2.5. VMD Feature Augmentation

Variational Mode Decomposition (VMD), proposed by [58], decomposes a time series into a finite set of band-limited intrinsic mode functions (IMFs) by solving a constrained variational optimization problem. Unlike Empirical Mode Decomposition (EMD), VMD simultaneously estimates all modes and their center frequencies by minimizing total bandwidth while ensuring their sum reconstructs the original signal. The optimization is solved using the Alternating Direction Method of Multipliers (ADMM) [59] [pp. 13–20]. Empirical evidence across financial and commodity markets confirms that VMD can separate trend, periodic, and disturbance components in non-stationary series, enabling downstream recurrent models to outperform standalone econometric and deep-learning benchmarks [25]. In commodity-price forecasting specifically, VMD-based hybrid models have been shown to outperform both single-model and earlier decomposition-based alternatives, with further gains when structurally relevant long-run drivers are incorporated into the forecasting system [31].

Let $\Delta p_t$ denote the differenced target series. The VMD objective is to find K modes ${\left\{u_k\right\}}_{k=1}^K$ and their center frequencies ${\left\{{\omega }_k\right\}}_{k=1}^K$ by solving:

$$\begin{array}{cc}\hfill \underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{min}& \sum _{k=1}^K{∥{\partial }_t\left[\left(\delta \left(t\right)+{\displaystyle \frac{j}{\pi t}}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}∥}_2^2\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \sum _{k=1}^K{u}_k\left(t\right)=f\left(t\right)\hfill \end{array}$$

(4)

Recent evidence from environmental forecasting further confirms that coupling VMD with sequence models can substantially improve predictive accuracy by separating short- and long-run signal structures before learning, relative to models trained on the raw series alone [32]. Related work on wind-power forecasting similarly demonstrates that VMD can improve downstream deep-learning performance by decomposing complex time series into more learnable components, while attention-based encoder-decoder structures further enhance the selection of informative features and time steps [35] [pp. 922–923]. (Full decomposition results are presented in Section 3.2.)

VMD-as-features design. Unlike the conventional decomposition-and-ensemble pipeline — which trains a separate model per IMF and aggregates predictions, an approach prone to cumulative estimation error and variance collapse on differenced series [7] [pp. 6] Shumailov et al. [46] [pp. 755–757] provide theoretical grounding for this failure mode, demonstrating that compound pipelines in which successive model stages each introduce finite-sample approximation error will systematically erode distributional variance — a result formalized in their Gaussian collapse theorem showing that output variance converges to zero almost surely as the number of estimation stages increases. Therefore, this study appends the five IMF series directly to the input feature matrix. The augmented input at time step t is:

$$x_t=\left[e_t\parallel u_{1,t}\parallel u_{2,t}\parallel u_{3,t}\parallel u_{4,t}\parallel u_{5,t}\right]\in {\mathbb{R}}^{29}$$

(5)

where $e_t\in {\mathbb{R}}^{24}$ is the vector of 24 economic features and $u_{k,t}$ is the value of the k-th IMF at time t. This design is consistent with recent forecasting evidence showing that VMD can enhance downstream sequence models by isolating heterogeneous temporal components and reducing the burden placed on a single learner to model mixed-frequency dynamics [32]. This architectural choice follows a broader pattern in hybrid forecasting research, where decomposition, feature reorganization, and downstream nonlinear learners are jointly designed to improve performance on complex non-stationary signals [36,45].

This design preserves multi-scale frequency information within a single forward pass; ablation experiments reported in Section 3.3 demonstrate that the conventional per-IMF approach yields substantially lower StdR on the differenced series. Prior VMD-based forecasting research also suggests that different frequency components are not equally suited to a single predictive mechanism, motivating architectures that can adapt to heterogeneous temporal structures [30]. Recent FX forecasting evidence further indicates that VMD-enhanced bidirectional recurrent architectures, when paired with systematic hyper-parameter tuning and regularization, can materially improve out-of-sample performance in low-frequency financial time series [37]. Evidence from natural rubber futures forecasting specifically shows that VMD can extract economically meaningful multi-scale components, while bidirectional recurrent architectures improve both fitting performance and directional prediction; importantly, the predictive value of individual modes depends on their time-scale correspondence with the forecast target [7]. The use of VMD before sequence learning is also consistent with prior work showing that decomposition can enrich the input representation for attention-based recurrent architectures in non-stationary forecasting tasks [35]. Recent oil-price forecasting studies have furthermore moved beyond single-stage decomposition toward secondary decomposition-reconstruction-ensemble frameworks, arguing that the information embedded in high-frequency components is too rich to be exhausted by a single decomposition pass [48]; the present VMD-as-features design offers a more parsimonious alternative while still exploiting multi-scale structure within a single forward pass.

2.3. Dual-Pathway Model Architecture

Given the augmented input matrix $X=\left[x_{t-L+1},\dots ,x_t\right]\in {\mathbb{R}}^{L\times 29}$ with look-back window $L=30$ trading days, two independent encoding pathways operate in parallel. Their outputs are fused by a learned Fusion Head to generate the scalar price-change forecast. The complete architecture is summarized in Table 5 and illustrated in Figure 3.

2.3.1. Path 1: BiLSTM with Temporal Attention

A two-layer Bidirectional LSTM [60] with hidden size $H=128$ per direction processes the input sequence simultaneously in the forward and backward temporal directions. The concatenated hidden state at each step $\tau$ is:

$$h_{\tau }=\left[h_{\tau }^{\to };h_{\tau }^{\leftarrow }\right]\in {\mathbb{R}}^{2H}$$

(6)

where $h_{\tau }^{\to }$ and $h_{\tau }^{\leftarrow }$ are the forward and backward hidden states respectively, and $H=128$, yielding a 256-dimensional representation. A learned temporal attention mechanism aggregates hidden states into a context vector:

$${\alpha }_{\tau }=\mathrm{softmax}\left(tanh\left(W_a·h_{\tau }+b_a\right)\right)\in {\mathbb{R}}^L$$

(7)

$$c=\sum _{\tau =1}^L{\alpha }_{\tau }·h_{\tau }\in {\mathbb{R}}^{2H}$$

(8)

where $W_a\in {\mathbb{R}}^{2H\times 1}$ and $b_a\in \mathbb{R}$ are learnable parameters. The context vector $c\in {\mathbb{R}}^{256}$ summarizes the entire sequence with emphasis on the most informative time steps.

2.3.2. Path 2: Transformer Encoder

A two-layer Transformer encoder [38] operates on a linearly projected version of the input sequence. Each input token is first projected to a $d_{\mathrm{model}}=128$ dimensional space and augmented with a learnable positional embedding:

$$z_{\tau }=W_{\mathrm{proj}}·x_{\tau }+p_{\tau }\in {\mathbb{R}}^{d_{\mathrm{model}}}$$

(9)

where $x_{\tau }$ is the input feature vector, $W_{\mathrm{proj}}\in {\mathbb{R}}^{d_{\mathrm{model}}\times 29}$ is the projection matrix, $p_{\tau }\in {\mathbb{R}}^{d_{\mathrm{model}}}$ is the learnable positional encoding for position $\tau$, and $d_{\mathrm{model}}=128$. Learnable positional encodings are preferred over fixed sinusoidal encodings because financial time series exhibit irregular temporal patterns that may not conform to fixed periodic representations.

Each Transformer encoder layer applies multi-head self-attention followed by a position-wise feed-forward network with residual connections and layer normalization:

$$\mathrm{Attention}(Q,K,V)=\mathrm{softmax}\left({\displaystyle \frac{Q{K}^{\top }}{\sqrt{d_k}}}\right)V$$

(10)

where $Q\in {\mathbb{R}}^{L\times d_k}$, $K\in {\mathbb{R}}^{L\times d_k}$, and $V\in {\mathbb{R}}^{L\times d_v}$ are the query, key, and value matrices derived from linear projections of the token embeddings, and $d_k$ is the dimension of the key vectors. The model uses $n_{\mathrm{heads}}=4$ attention heads, feed-forward dimension 512, and dropout rate 0.1. The final-step representation $z_L\in {\mathbb{R}}^{128}$ (last token of the encoder output) is taken as the Transformer pathway output, encoding long-range temporal dependencies across the entire look-back window. This final-position representation encodes a summary of the full look-back window, playing a role similar to the [CLS] classification token, which serves as a sequence-level summary vector in BERT-style encoders.

Figure 3. VMD-Hybrid BiLSTM-Transformer Architecture Diagram

Figure 3. VMD-Hybrid BiLSTM-Transformer Architecture Diagram

2.3.3. Fusion Head and Output Layer

The outputs of the two encoding pathways are concatenated to form a joint representation:

$$v=\left[c;z_L\right]\in {\mathbb{R}}^{384}$$

(11)

The fused vector is passed through a multi-layer projection network to generate the scalar price-change forecast: LayerNorm(384) → Linear($384\to 128$) → GELU → Dropout($p=0.2$) → Linear($128\to 64$) → GELU → Linear($64\to 1$). GELU activation is adopted because it provides smoother gradient flow than ReLU for regression tasks on continuous financial targets, consistent with its use in recent Transformer-based forecasting architectures [41,42,43,44].

The scalar output ${\widehat{y}}_t$ represents the predicted normalized price change $\Delta {\tilde{p}}_t$. The actual price change is recovered via inverse MinMax transformation (Equation 3), and the corresponding price level is reconstructed using a rolling one-step scheme: ${\widehat{p}}_t=p_{t-1}+\Delta {\widehat{p}}_t$.

2.4. Model Training

The model is trained end-to-end on the Stage 3 training partition (n = 2,140 observations, 2018–2023) using the following configuration:

Table 6. Training hyperparameters.

Hyperparameter	Value	Rationale
Loss function	Huber Loss ($\delta =0.5$)	Robust to outlier spikes
Optimiser	AdamW (lr $=5\times {10}^{-4}$, wd $={10}^{-4}$)	Adaptive lr + weight decay
LR scheduler	ReduceLROnPlateau ($\times 0.5$, patience 10)	Reduce on validation plateau
Early stopping	Patience $=30$ epochs	Prevent overfitting
Max epochs	300	Upper bound
Batch size	32	Mini-batch SGD
Lookback L	30 trading days	$\sim 6$ trading weeks of context
Total parameters	1,053,698

Table 6. Training hyperparameters.

Hyperparameter	Value	Rationale
Loss function	Huber Loss ($\delta =0.5$)	Robust to outlier spikes
Optimiser	AdamW (lr $=5\times {10}^{-4}$, wd $={10}^{-4}$)	Adaptive lr + weight decay
LR scheduler	ReduceLROnPlateau ($\times 0.5$, patience 10)	Reduce on validation plateau
Early stopping	Patience $=30$ epochs	Prevent overfitting
Max epochs	300	Upper bound
Batch size	32	Mini-batch SGD
Lookback L	30 trading days	$\sim 6$ trading weeks of context
Total parameters	1,053,698

The Huber loss with threshold $\delta$ is defined as:

$${\mathcal{L}}_{\delta }(y,\widehat{y})=\left\{\begin{array}{cc}{\displaystyle \frac{1}{2}}{(y-\widehat{y})}^2\hfill & \mathrm{if}|y-\widehat{y}|\le \delta \hfill \\ \delta \left(|y-\widehat{y}|-{\displaystyle \frac{1}{2}}\delta \right)\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(12)

The choice of $\delta =0.5$ was determined by grid search over $\delta \in \{0.1,0.3,0.5,1.0\}$ on the validation set; $\delta =0.5$ provides quadratic sensitivity for small prediction errors while bounding the influence of large price shocks on gradient updates. The lookback window $L=30$ was selected to cover approximately six trading weeks, capturing both intra-month cyclical patterns and cross-market lag effects without introducing excessive computational cost; sensitivity to L is reported in Appendix F.

2.5. Evaluation Protocol, Baselines, and Multi-step Extension

Model performance is assessed across two complementary evaluation spaces. Primary metrics are computed on the predicted and actual first-differenced series ($\Delta {\widehat{p}}_t$ vs. $\Delta p_t$) in real Baht/kg/day units after inverse transformation (Equation 3). These directly measure the model’s ability to predict daily price movements. Supplementary price-level metrics (MAE, RMSE, MAPE, $R^2$) are computed on the reconstructed series ${\widehat{p}}_t=p_{t-1}+\Delta {\widehat{p}}_t$ and are reported for completeness only; their values are inflated by an autoregressive anchor effect and do not reflect additional predictive skill in the differenced space.

In Table 7, Pearson r and StdR are reported as co-primary metrics following the evaluation framework of Taylor [50] [pp. 7183–7184], who shows that neither metric is sufficient alone — high correlation does not guarantee correct amplitude, and low RMS error does not imply that forecast variance faithfully reproduces the target series — and consistent with Mehdiyev et al. [61] [pp. 264–265] and Kyriakidis et al. [62] [pp. 104–105], who demonstrate that robust forecast evaluation requires concurrent consideration of multiple criteria because model rankings vary substantially across metrics.

DA is computed directly on the differenced series — $\mathrm{sign}(\Delta {\widehat{p}}_t)$ vs. $\mathrm{sign}(\Delta p_t)$ — not on the reconstructed price level. Days on which $\Delta p_t=0$ are excluded from the denominator ($n^{*}$). Following Zhu et al. [7] [p. 6], DA $>0.6$ is regarded as the threshold for practically meaningful directional prediction. However, this threshold should be interpreted alongside confusion-matrix and variance-ratio diagnostics: Costantini et al. [53] show that raw DA becomes materially informative only when paired with magnitude-sensitive directional value measures, while Rajpal et al. [47] [pp. 13–15] demonstrate that symmetric loss objectives can produce apparently acceptable DA through one-sided prediction that is fully revealed only by specificity and StdR diagnostics.

The use of a ratio-based dispersion metric as a collapse diagnostic is consistent with the Variability Collapse Index proposed by Xu and Liu [63], who derive a formally analogous quantity from the minimum MSE linear probing loss and demonstrate that ratio-based metrics are invariant to invertible linear transformations of the feature space — a property that point-estimate collapse measures such as within-class distance do not possess. Four baseline models are evaluated on the identical test partition ($n=237$ observations, 18 September 2025 – 27 February 2026) using the same protocol.

Table 8. Baseline model specifications.

Baseline	Specification	Purpose
Naive No-Change	$\Delta {\widehat{p}}_t=0\forall t$	Lower bound; tests if any model beats zero-change forecast
Naive Random Walk	$\Delta {\widehat{p}}_t=\Delta p_{t-1}$	Standard persistence baseline for differenced series
ARIMA(2,0,2)	AIC-selected via grid search on training set; rolling 1-step refitting on test set	Linear time-series benchmark*
Vanilla LSTM	Unidirectional 2-layer LSTM (hidden $=128$); same 29-dim input; no VMD decomposition stage	Ablates the VMD-as-features contribution

Notes: * ARIMA serves as the natural linear benchmark because it dominates the early rubber forecasting literature — [12] identify ARIMA(1,1,0) as the best-fitting model for Malaysian SMR20, and [13] use ARIMA as the univariate baseline against which structural supply–demand models are evaluated — making it the standard against which forecasting advances in this commodity are measured.

Table 8. Baseline model specifications.

Baseline	Specification	Purpose
Naive No-Change	$\Delta {\widehat{p}}_t=0\forall t$	Lower bound; tests if any model beats zero-change forecast
Naive Random Walk	$\Delta {\widehat{p}}_t=\Delta p_{t-1}$	Standard persistence baseline for differenced series
ARIMA(2,0,2)	AIC-selected via grid search on training set; rolling 1-step refitting on test set	Linear time-series benchmark*
Vanilla LSTM	Unidirectional 2-layer LSTM (hidden $=128$); same 29-dim input; no VMD decomposition stage	Ablates the VMD-as-features contribution

Notes: * ARIMA serves as the natural linear benchmark because it dominates the early rubber forecasting literature — [12] identify ARIMA(1,1,0) as the best-fitting model for Malaysian SMR20, and [13] use ARIMA as the univariate baseline against which structural supply–demand models are evaluated — making it the standard against which forecasting advances in this commodity are measured.

Multi-Step Forecasting Extension. 

In addition to the one-step-ahead horizon ($h=1$), this study evaluates the proposed model under a direct multi-step forecasting strategy for horizons $h\in \{1,2,3,5,10,20,30\}$ trading days. This study adopts the direct forecasting strategy, training a separate model for each horizon h to avoid error accumulation. While this requires H independent models, it produces unbiased h-step-ahead forecasts without recursive error propagation [64]. Under the direct strategy, the training target is shifted accordingly:

$$y_i^{\left(h\right)}=\Delta p_{i+L+h-1},h=1,2,\dots ,H$$

(13)

where $y_i^{\left(h\right)}$ denotes the target value for the h-step-ahead forecast associated with sample i, defined as the realized price change $\Delta p_{i+L+h-1}$ at forecast horizon h.

The emphasis on short lead times is consistent with prior VMD-based forecasting studies, where predictive gains from decomposition are strongest at near-term horizons and tend to weaken as the forecast horizon extends [33]. For $h>1$, future exogenous features are unavailable at prediction time; the model receives only the values observed within the look-back window and makes no use of forward-looking information. This is equivalent to assuming that all exogenous conditions remain constant beyond the last observation, a limitation discussed in Section 4.

During the preparation of this section, the authors used Claude (Anthropic) to assist with manuscript drafting and editing, data preprocessing, and visualization. The authors have reviewed and edited all AI-generated content and take full responsibility for the content of this publication.

3. Results

This section reports the empirical results of the proposed VMD–Hybrid BiLSTM–Transformer framework. All primary evaluations are conducted on the Stage 3 out-of-sample test set ($n=237$ daily observations, 18 September 2025 – 27 February 2026). Section 3.2 presents the VMD signal decomposition; Section 3.3 reports the ablation analysis confirming the contribution of each architectural component; Section 3.4 reports out-of-sample forecast performance against benchmark models; Section 3.5 presents directional bias diagnostics; and Section 3.6 reports multi-step forecast skill across horizons $h\in \{1,2,3,5,10,20,30\}$.

3.1. Evaluation Framework

The primary evaluation is conducted in differenced space using five metrics. Directional Accuracy (DA) measures the proportion of days on which the predicted sign of the price change matches the realized sign; days on which $\Delta p_t=0$ are excluded from the denominator. Pearson correlation (r) quantifies the linear association between predicted and actual price changes. MAE and RMSE measure average and root-mean-squared prediction error in Baht/kg/day. Standard Deviation Ratio (StdR) is defined as:

$$\mathrm{StdR}={\displaystyle \frac{\mathrm{std}(\Delta \widehat{p})}{\mathrm{std}(\Delta p)}}$$

(14)

StdR equals 1.0 under ideal forecast dispersion and approaches 0 when a model produces near-constant predictions regardless of actual variation. A value below 0.20 is taken as evidence of variance collapse — a failure mode in which a model learns the unconditional mean rather than the conditional distribution of price changes.This convergence toward the unconditional mean is the output-level analogue of the distributional concentration documented by Shumailov et al. [46] [pp. 755–757], who prove that models operating under compound estimation error systematically suppress tail events and converge toward near-zero variance outputs — and of the one-sided prediction regimes identified by Rajpal et al. [47] [pp. 13–15], in which symmetric loss objectives produce degenerate directional outputs whose apparent performance masks near-zero predictive content. Rajpal et al. [47] further show that this degenerate regime produces severe imbalance between recall and specificity — a pathology that high DA alone cannot detect, reinforcing the need for confusion-matrix and variance-ratio diagnostics alongside directional hit-rate statistics. StdR is the most important diagnostic for differenced commodity price models because high DA can be achieved trivially by predicting the majority class, whereas reproducing the variance of the target series requires genuine predictive information.

Price-level metrics — specifically $R^2$ and MAPE — are included in Table 10 for completeness only. These statistics are inflated by the autoregressive anchor effect: because the reconstructed level ${\widehat{p}}_t=p_{t-1}+\Delta {\widehat{p}}_t$ inherits the previous day’s price regardless of forecast quality, $R^2$ and MAPE at the level will appear favourable even when the differenced-space forecast is uninformative. They therefore carry no additional diagnostic value beyond the five primary metrics reported in Table 7.

Figure 4. VMD decomposition of the differenced price series (K=5)

Figure 4. VMD decomposition of the differenced price series (K=5)

3.2. VMD Decomposition of the RSS3 Price Series

Variational Mode Decomposition was applied to the differenced RSS3 FOBm1 series ($\Delta p_t$) over the Stage 3 training and validation partitions ($n=2,407$ observations). The decomposition with $K=5$ modes (see Appendix A) achieved a reconstruction RMSE of 0.054 normalized units, representing 3.2% of the series standard deviation — confirming that the five components collectively reconstruct the original differenced series with minimal residual error ($\sum \mathrm{IMF}\approx \Delta p_t$).

Table 9. VMD decomposition results — Stage 3, $K=5$.

IMF	Label	Center freq.	Energy (%)	Economic interpretation
1	Trend	$f\approx 0.0004$	78.6%	Long-run structural price trend (macroeconomic cycle, supply shifts)
2	Low-freq.	$f\approx 0.069$	8.0%	Weekly–monthly market cycles (export dynamics, inventory)
3	Low-freq.	$f\approx 0.147$	3.9%	Bi-weekly fluctuations (tapping seasons, policy announcements)
4	High-freq.	$f\approx 0.328$	2.7%	Short-term trading noise (intraday, speculative positions)
5	High-freq.	$f\approx 0.471$	6.7%	High-frequency microstructure noise
Total		Reconstruction RMSE $=0.054$		Near-perfect reconstruction

Table 9. VMD decomposition results — Stage 3, $K=5$.

IMF	Label	Center freq.	Energy (%)	Economic interpretation
1	Trend	$f\approx 0.0004$	78.6%	Long-run structural price trend (macroeconomic cycle, supply shifts)
2	Low-freq.	$f\approx 0.069$	8.0%	Weekly–monthly market cycles (export dynamics, inventory)
3	Low-freq.	$f\approx 0.147$	3.9%	Bi-weekly fluctuations (tapping seasons, policy announcements)
4	High-freq.	$f\approx 0.328$	2.7%	Short-term trading noise (intraday, speculative positions)
5	High-freq.	$f\approx 0.471$	6.7%	High-frequency microstructure noise
Total		Reconstruction RMSE $=0.054$		Near-perfect reconstruction

Figure 5. Predicted vs Actual — Test Period (VMD-Hybrid BiLSTM-Transformer)

Figure 5. Predicted vs Actual — Test Period (VMD-Hybrid BiLSTM-Transformer)

Upper panel: daily first-differenced predictions vs. actual (Baht/kg/day). Lower panel: reconstructed price level via rolling one-step accumulation (${\widehat{p}}_t=p_{t-1}+\Delta {\widehat{p}}_t$).

IMF1 accounts for 78.6% of total signal energy and exhibits the lowest center frequency ($f\approx 0.0004$), consistent with the long-run structural trends documented in the natural rubber literature, including fluctuations in global automobile production and plantation supply cycles [7]. The two intermediate-frequency components (IMF2–IMF3) collectively capture 11.9% of energy and correspond to temporal scales of one to four weeks, consistent with the inventory adjustment cycles and seasonal tapping rhythms documented in rubber export data from Thailand, Indonesia, and Malaysia [7,30]. The structural dominance of IMF₁ is further consistent with the simultaneous supply–demand framework of Arunwarakorn et al. [55] [pp. 8–9], who show that equilibrium rubber prices are anchored by slow-moving fundamentals — plantation area, synthetic rubber prices, GDP growth, and crude oil — operating over multi-year cycles rather than at the weekly or daily frequencies captured by IMF₂–IMF₅. The high-frequency residuals (IMF4–IMF5) account for 9.4% of energy and reflect intraday speculative activity and short-term information shocks operating at sub-weekly scales.

The top panel shows the original $\Delta p_t$, while the remaining panels present IMF₁–IMF₅ with their centre frequencies and energy shares. IMF₁ dominates the signal (78.6% energy), capturing the low-frequency trend, whereas IMF₄–IMF₅ mainly represent high-frequency noise.

3.3. Ablation Study

To verify the contribution of each architectural component before reporting full model performance, an ablation analysis was conducted by training four variants of the proposed model on the Stage 3 training partition and evaluating each on the identical test set ($n=237$).

Table 10. Ablation study results — Test set ($n=237$).

Model variant	DA%	Corr	StdR	Key finding
Per-IMF BiLSTM+Transformer (conventional VMD pipeline)	$\sim 58\%$	$\sim 0.16$	$0.043$–$0.15$	Variance collapse on diff series
VMD-as-features + BiLSTM only (no Transformer)	$\sim 63\%$	$\sim 0.75$	$\sim 0.70$	Good but Transformer adds Corr
VMD-as-features + Transformer only (no BiLSTM)	$\sim 61\%$	$\sim 0.72$	$\sim 0.68$	Local patterns captured less well
VMD-as-features + BiLSTM + Transformer (full model) ★	$67.1\%★$	$0.812★$	$0.819★$	Best on all metrics

Table 10. Ablation study results — Test set ($n=237$).

Model variant	DA%	Corr	StdR	Key finding
Per-IMF BiLSTM+Transformer (conventional VMD pipeline)	$\sim 58\%$	$\sim 0.16$	$0.043$–$0.15$	Variance collapse on diff series
VMD-as-features + BiLSTM only (no Transformer)	$\sim 63\%$	$\sim 0.75$	$\sim 0.70$	Good but Transformer adds Corr
VMD-as-features + Transformer only (no BiLSTM)	$\sim 61\%$	$\sim 0.72$	$\sim 0.68$	Local patterns captured less well
VMD-as-features + BiLSTM + Transformer (full model) ★	$67.1\%★$	$0.812★$	$0.819★$	Best on all metrics

The ablation results confirm that all three components — the VMD-as-features design, the BiLSTM pathway, and the Transformer pathway — are necessary for optimal performance. The conventional per-IMF pipeline, which trains a separate model for each intrinsic mode function and sums the predictions, substantially mitigates the variance collapse problem relative to the full model: StdR drops from 0.819 to 0.087, confirming that the per-component summation approach fails to reproduce the dispersion of daily price changes on the differenced series.

Removing the Transformer pathway (BiLSTM-only variant) reduces Pearson r from 0.812 to 0.751, indicating that the global self-attention mechanism captures temporal dependencies that the recurrent pathway alone cannot represent. Removing the BiLSTM pathway (Transformer-only variant) reduces Pearson r to 0.718 and StdR to 0.682, confirming that local sequential modeling of adjacent time steps — a strength of recurrent architectures — provides complementary information to global attention. The full model’s dual-pathway design, processing the same 29-dimensional augmented input through parallel BiLSTM and Transformer encoders, yields the highest values on all three primary metrics, confirming that the two pathways are genuinely complementary rather than redundant.This result is consistent with the theoretical prediction of Shumailov et al. [46] [p. 757], whose Gaussian collapse theorem demonstrates that variance converges to zero almost surely when successive estimation stages each introduce finite-sample approximation error — a mechanism structurally equivalent to the cumulative error compounding that occurs when per-IMF forecasts are independently fitted and summed. This pattern is consistent with prior hybrid forecasting studies showing that performance gains arise from component complementary rather than from decomposition alone [36].

Figure 6. Ablation study — DA%, Pearson r, and StdR across model variants. The full VMD-Hybrid BiLSTM-Transformer (★) achieves the highest performance on all three metrics. The per-IMF conventional pipeline exhibits near-complete variance collapse (StdR $=0.043$–$0.15$).

Figure 6. Ablation study — DA%, Pearson r, and StdR across model variants. The full VMD-Hybrid BiLSTM-Transformer (★) achieves the highest performance on all three metrics. The per-IMF conventional pipeline exhibits near-complete variance collapse (StdR $=0.043$–$0.15$).

3.4. Out-of-Sample Forecast Performance

Having confirmed the contribution of each component through ablation, this section evaluates the full VMD-Hybrid model against four external benchmarks on the identical test partition ($n=237$).

The Naive No-Change model ($\Delta {\widehat{p}}_t=0$) serves as a lower bound; any useful forecasting model must exceed it on directional accuracy. The Naive Random Walk ($\Delta {\widehat{p}}_t=\Delta p_{t-1}$) tests whether recent momentum has predictive value. ARIMA(2,0,2), selected by AIC grid search, represents the linear time-series benchmark. Vanilla LSTM (unidirectional, 2-layer, same 29-dim input, no VMD) provides an additional ablation of the VMD-as-features contribution against a standard deep-learning baseline.

The proposed VMD-Hybrid model achieves the best performance across all primary metrics. Relative to ARIMA(2,0,2) — the strongest linear baseline — the hybrid model improves Pearson correlation by $+0.662$ (from 0.150 to 0.812), reduces MAE in the differenced space by 37.3% (from 0.249 to 0.156 Baht/kg/day), and reduces RMSE by 41.6% (from 0.443 to 0.258 Baht/kg/day). This margin is consistent with the broader rubber forecasting literature showing that purely univariate linear models are structurally insufficient once commodity-specific drivers are made explicit: Khin and Thambiah [13] demonstrate that a simultaneous supply–demand system outperforms ARIMA across all accuracy criteria, while Eng and Khalid [24] confirm that LSTM dominates ARIMA on Malaysian SMR20 — establishing that the gain over ARIMA reported here reflects a systematic pattern rather than a sample-specific outcome. Directional accuracy of 67.1% exceeds the random baseline (50%) by 17.1 percentage points, confirming genuine predictive skill beyond chance.

The model’s predictive skill is consistent with the layered cross-market structure documented by Kepulaje et al. [9] [pp. 159–161] and Pinitjitsamut [57]: because SHFE signals are incorporated into FOB Bangkok with a lag and pass-through remains incomplete and asymmetric in the short run, genuine one-day-ahead predictability is structurally limited — making a Pearson correlation of 0.812 on the differenced series a substantively strong result rather than an expected one.

Figure 7. Benchmark comparison of forecasting models. Panels report (a) directional accuracy (DA), (b) Pearson r, (c) ${\mathrm{MAE}}_{\mathrm{diff}}$ (Baht/kg/day), and (d) the variance ratio (StdR). Dashed lines denote the 50% random baseline (DA) and the ideal value StdR $=1$. The proposed VMD-Hybrid model consistently outperforms ARIMA, Vanilla LSTM, and naïve benchmarks. The StdR panel also reveals variance collapse in baseline models, clarifying why high DA in Vanilla LSTM can be misleading.

Figure 7. Benchmark comparison of forecasting models. Panels report (a) directional accuracy (DA), (b) Pearson r, (c) ${\mathrm{MAE}}_{\mathrm{diff}}$ (Baht/kg/day), and (d) the variance ratio (StdR). Dashed lines denote the 50% random baseline (DA) and the ideal value StdR $=1$. The proposed VMD-Hybrid model consistently outperforms ARIMA, Vanilla LSTM, and naïve benchmarks. The StdR panel also reveals variance collapse in baseline models, clarifying why high DA in Vanilla LSTM can be misleading.

The StdR of 0.819 indicates that the VMD-Hybrid model produces predictions whose variance is close to that of the actual series (StdR $=1.0$ ideal), demonstrating that the model captures the magnitude of daily price movements in addition to their direction. This contrasts sharply with the Naive No-Change (StdR $=0.000$) and Vanilla LSTM (StdR $=0.033$) baselines, which exhibit near-complete variance collapse — producing near-constant output regardless of actual market conditions.

The price-level metrics reported in Table 11 should be interpreted with caution. The Naive No-Change model, which predicts zero price change on every day, achieves $R^2=0.977$ and MAPE $=0.32\%$ — only marginally below the VMD-Hybrid ($R^2=0.991$, MAPE $=0.23\%$) — despite having zero predictive content (DA $=0\%$, Corr $=0.000$). This confirms that price-level metrics alone cannot distinguish genuine forecasting skill from a trivial baseline, as discussed in Section 3.1.

This is precisely the evaluation failure described by Taylor [50] [pp. 7183–7184]: a forecast can appear numerically accurate while being dynamically wrong if it suppresses variance rather than tracking the stochastic structure of the series. Ampountolas [65] [pp. 2, 14–18] illustrates the same gap in commodity forecasting practice — standard MAE, MSE, and RMSE comparisons identify the model with smallest average error but cannot reveal whether the winning model reproduces distributional scale, making StdR and Pearson r indispensable complements rather than optional additions.

3.5. Directional Bias Analysis

Having established that the VMD-Hybrid model outperforms all baselines on primary metrics, this section examines why Vanilla LSTM’s nominally higher directional accuracy of 78.6% does not reflect genuine forecasting skill.

Table 12. Directional confusion matrix — Test set ($n=237$).

Model	Pred ↑ (UP)	Pred ↓ (DOWN)	DA (%)	StdR
Naive No-Change	TP $=0$	TN $=46$	0.0%	0.000
ARIMA(2,0,2)	TP $=115$	TN $=26$	59.7%	0.339
Vanilla LSTM †	TP $=186$	TN $=1$	$78.6{\%}^{†}$	0.033
VMD-Hybrid ★	TP $=123$	TN $=36$	$67.1\%★$	$0.819★$

Note: Rows = actual direction; columns = predicted direction. DA computed on non-zero actual days only. TP = True Positive (Actual↑, Predicted↑). TN = True Negative (Actual↓, Predicted↓). Actual UP days $=191$ (80.6%); Actual DOWN days $=46$ (19.4%).

Table 12. Directional confusion matrix — Test set ($n=237$).

Model	Pred ↑ (UP)	Pred ↓ (DOWN)	DA (%)	StdR
Naive No-Change	TP $=0$	TN $=46$	0.0%	0.000
ARIMA(2,0,2)	TP $=115$	TN $=26$	59.7%	0.339
Vanilla LSTM †	TP $=186$	TN $=1$	$78.6{\%}^{†}$	0.033
VMD-Hybrid ★	TP $=123$	TN $=36$	$67.1\%★$	$0.819★$

Note: Rows = actual direction; columns = predicted direction. DA computed on non-zero actual days only. TP = True Positive (Actual↑, Predicted↑). TN = True Negative (Actual↓, Predicted↓). Actual UP days $=191$ (80.6%); Actual DOWN days $=46$ (19.4%).

During the test period, 80.6% of trading days saw a price increase (191 up-days vs. 46 down-days). Vanilla LSTM exploited this imbalance by predicting UP on 97.1% of all days (231 out of 237 non-zero days), achieving TN $=1$ — meaning the model correctly identified only one price decrease in the entire test set. With StdR $=0.033$, the Vanilla LSTM output is effectively a near-constant signal, consistent with a model that has learned the unconditional sample mean of the training set rather than dynamic conditional patterns.

Figure 8. Predicted versus actual price changes in the test set. Scatter plot of predicted daily price changes $\Delta {\widehat{p}}_t$ against observed changes $\Delta p_t$ for the Stage 3 test sample ($n=237$). The $45°$ reference line represents perfect prediction. Data points are colored by chronological quartile (Q1 = blue to Q4 = red). The strong linear association is reflected in $r=0.812$.

Figure 8. Predicted versus actual price changes in the test set. Scatter plot of predicted daily price changes $\Delta {\widehat{p}}_t$ against observed changes $\Delta p_t$ for the Stage 3 test sample ($n=237$). The $45°$ reference line represents perfect prediction. Data points are colored by chronological quartile (Q1 = blue to Q4 = red). The strong linear association is reflected in $r=0.812$.

By contrast, the VMD-Hybrid model predicted DOWN on 104 occasions (43.9% of days), achieving TN $=36$ out of 46 actual down-days — a recall for the down-direction of 78.3%. The balanced confusion matrix confirms that the model has learned genuine directional dynamics rather than a static positive bias. Collectively, the confusion matrix diagnostics and StdR comparison confirm that Vanilla LSTM’s superior DA reflects a degenerate positive-prediction strategy rather than genuine forecasting ability, and that Pearson $r=0.812$ is the appropriate primary performance indicator for this series. This finding corroborates Bürgi [51] [pp. 7909–7912], who demonstrates that DA is not a self-sufficient evaluation concept without knowledge of magnitude sensitivity and class composition, and McCarthy and Snudden [52] [pp. 1–7], who show that success ratios can exceed 0.5 mechanically under class imbalance, generating pseudo-skill where no genuine predictability exists — precisely the pathology illustrated here by Vanilla LSTM’s TN $=1$ result. Costantini et al. [53] sharpen this point further by showing that DA becomes materially informative only when paired with magnitude-sensitive directional value — a measure that weights correct sign calls by the size of the realized move rather than treating all correct predictions as equivalent. Under that criterion, Vanilla LSTM’s near-constant output would register near-zero directional value despite its nominally high hit rate, consistent with its StdR $=0.033$.

Figure 9. Directional bias analysis using confusion matrices. (a) Naive No-Change; (b) Naive Random Walk; (c) ARIMA(2,0,2); (d) Vanilla LSTM; (e) VMD-Hybrid

Figure 9. Directional bias analysis using confusion matrices. (a) Naive No-Change; (b) Naive Random Walk; (c) ARIMA(2,0,2); (d) Vanilla LSTM; (e) VMD-Hybrid

3.6. Multi-Step Forecast Skill Degradation

Having validated the model’s one-step-ahead performance and component contributions, this section assesses forecast skill across longer horizons to establish the practical range of reliable predictions.

This study adopts the direct forecasting strategy, training a separate model with identical architecture and hyper-parameters for each horizon h, to avoid the error accumulation associated with recursive multi-step forecasting [64]. While this requires H independent models, it produces unbiased h-step-ahead forecasts without recursive error propagation. A separate model was trained for each horizon $h\in \{1,2,3,5,10,20,30\}$.

Table 13. Multi-step direct forecast performance — VMD-Hybrid BiLSTM-Transformer.

h	DA%	Pearson r	MAE (Bt/kg/d)	RMSE	StdR
1 ★	67.1%	0.812	0.156	0.258	0.819
2	86.2%	0.928	0.622	—	0.836
3	82.1%	0.827	0.949	—	0.733
5	77.2%	0.590	1.257	—	0.581
10	81.7%	−0.007	1.136	—	0.006
20	81.7%	−0.057	1.151	—	0.011
30	81.7%	0.069	1.141	—	0.000

Table 13. Multi-step direct forecast performance — VMD-Hybrid BiLSTM-Transformer.

h	DA%	Pearson r	MAE (Bt/kg/d)	RMSE	StdR
1 ★	67.1%	0.812	0.156	0.258	0.819
2	86.2%	0.928	0.622	—	0.836
3	82.1%	0.827	0.949	—	0.733
5	77.2%	0.590	1.257	—	0.581
10	81.7%	−0.007	1.136	—	0.006
20	81.7%	−0.057	1.151	—	0.011
30	81.7%	0.069	1.141	—	0.000

Figure 10. Multi-step forecast skill degradation. Directional accuracy (DA, left axis) and Pearson r (right axis) across forecast horizons $h=1,2,3,5,10,20,30$. The dashed line indicates the random benchmark (DA $=50\%$). The logarithmic x-axis highlights the non-linear decline in predictive skill as the horizon increases. The $h=2$ anomaly reflects sustained directional trends in the test period rather than a systematic model property; see main text for discussion.

Figure 10. Multi-step forecast skill degradation. Directional accuracy (DA, left axis) and Pearson r (right axis) across forecast horizons $h=1,2,3,5,10,20,30$. The dashed line indicates the random benchmark (DA $=50\%$). The logarithmic x-axis highlights the non-linear decline in predictive skill as the horizon increases. The $h=2$ anomaly reflects sustained directional trends in the test period rather than a systematic model property; see main text for discussion.

The $h=2$ horizon produces higher DA (86.2%) and Pearson r (0.928) than $h=1$, which requires interpretation. This pattern most likely reflects the specific composition of the test period: during September 2025 – February 2026, the RSS3 price series exhibited sustained directional trends over multi-day windows, making two-day-ahead movements more predictable than single-day changes which capture more micro-structure noise. The $h=2$ model effectively filters intraday reversals that generate short-run unpredictability at $h=1$, resulting in smoother targets that are more learnable given the available features. This interpretation is consistent with the IMF energy structure reported in Section 3.2, where high-frequency components (IMF₄–IMF₅) that contribute noise at $h=1$ are less prominent in the two-day-ahead target series.

This pattern also illustrates the model-selection hazard documented by Costantini et al. [53] and Costantini and Kunst [54]: relying on DA as a selection criterion at longer horizons would favour the $h=10$–30 models — whose DA of 81.7% exceeds the $h=1$ model’s 67.1% — despite those models having zero predictive content by every variance-sensitive measure, a finding consistent with their evidence that DA-weighted selection rarely yields robust gains and can worsen performance in MSE terms. Beyond $h=2$, Pearson r declines monotonically: from 0.827 at $h=3$ to 0.590 at $h=5$, and to near-zero at $h=10$ ($r=-0.007$). StdR follows the same pattern, collapsing to approximately 0.006 at $h=10$, indicating that the model can no longer reproduce the variance of ten-day-ahead price changes. DA values at $h=10$–30 (approximately 81.7%) are misleading for the same reason identified in Section 3.5: they reflect a positive-class majority in the cumulative 10–30-day return distribution rather than genuine predictive content, as confirmed by StdR $\approx 0.000$–$0.011$. The emphasis on short lead times is consistent with prior VMD-based forecasting studies, where predictive gains from decomposition are strongest at near-term horizons and tend to weaken as the forecast horizon extends [33].

These results indicate that the VMD-Hybrid model provides reliable forecasting skill up to a five-day horizon ($h=5$, $r=0.590$, StdR $=0.581$), with meaningful but diminishing skill at $h=3$. Beyond five trading days, the model should not be used for directional price guidance.

4. Discussion

The VMD decomposition results provide structural insight into the multi-scale price dynamics of the natural rubber market, and their economic interpretation reinforces the theoretical rationale for the proposed architectural design. The dominant trend component (IMF₁, 78.6% of total signal energy) captures slow-moving, long-run structural forces — principally global tyre demand driven by automobile production, plantation area adjustment in Thailand, Indonesia, and Malaysia, and the substitution relationship between natural and synthetic rubber. The pronounced concentration of signal energy in this single low-frequency mode confirms that RSS3 price dynamics are fundamentally anchored by supply–demand fundamentals operating over multi-year cycles, consistent with the commodity economics literature and with prior VMD-based evidence from natural rubber futures markets showing that the dominant mode reflects broader structural market conditions and longer-horizon price behavior [7].

The intermediate-frequency components (IMF₂–IMF₃, 11.9% combined energy) are consistent with medium-term market cycles arising from seasonal tapping patterns, inventory accumulation and draw-down, and semi-annual export policy adjustments in major producing countries. These components exhibit a degree of regularity that makes them particularly amenable to modeling by the BiLSTM pathway, whose recurrent structure is well-suited to capturing sequentially persistent, periodic fluctuations. This interpretation aligns with financial and energy-price forecasting evidence demonstrating that decomposition-based models improve predictive accuracy precisely by separating lower-frequency structural movements — which possess greater regularity — from higher-frequency disturbance components before sequence learning is applied [25,30].

The high-frequency residuals (IMF₄–IMF₅, 9.4% combined energy) capture speculative trading noise and short-horizon information shocks, including geopolitical disruptions, weather-related supply shocks, and the COVID-19-induced market dislocations observed in 2020–2021. The complementary role of the Transformer encoder is most apparent here: its global self-attention mechanism operates across the entire lookback window simultaneously, enabling the model to identify non-contiguous historical observations that are most informative for predicting irregular,short-duration deviations — a capability that strictly sequential recurrent processing cannot replicate. This division of labor between pathways is further supported by prior rubber-futures evidence indicating that high-frequency IMFs carry the most information relevant to near-term directional prediction, whereas low-frequency modes dominate longer-horizon behavior [7], and by energy-price forecasting research establishing that model performance improves when architectural design is explicitly adapted to the heterogeneous information content of components across frequency bands [30].

These decomposition-level findings translate directly into aggregate forecasting outcomes. The VMD-Hybrid model achieves a Directional Accuracy of 67.1% and a Pearson correlation of 0.812 on the held-out test set, representing practically meaningful skill for a daily commodity price series in which the random benchmark is 50%. A market participant employing the model’s directional signal would correctly anticipate the direction of the next trading day’s price movement approximately two-thirds of the time — a margin sufficient to support statistically informed hedging decisions, inventory management, and procurement timing across the natural rubber supply chain. These results are consistent with the broader hybrid-forecasting literature demonstrating that VMD-enhanced architectures improve predictive performance by reducing the stationary burden on downstream sequence learners, with gains most pronounced when decomposition is integrated within a unified framework capable of exploiting both local sequential dependencies and global temporal structure [32,45]. In the commodity-price domain specifically, this pattern is corroborated by evidence that VMD improves model learnability through multi-frequency isolation, and that the further inclusion of structurally relevant long-run economic drivers yields additional gains in both level and directional accuracy [31].

5. Conclusions

This study proposed the VMD-Hybrid BiLSTM-Transformer, a unified forecasting framework for daily RSS3 FOB natural rubber price changes that integrates Variational Mode Decomposition with a dual-pathway deep learning architecture comprising a Bidirectional LSTM encoder with temporal attention and a Transformer encoder. The central architectural contribution is the treatment of VMD intrinsic mode functions as input features appended directly to the economic feature matrix — rather than as independent prediction targets in a per-component summation pipeline. This design decision proved consequential: the conventional per-IMF approach produced near-constant forecasts on the differenced series, with Standard Deviation Ratios of 0.04–0.15, whereas the proposed feature-augmentation design achieved a StdR of 0.819, confirming that multi-scale frequency information is effectively preserved within a single forward pass and that forecast variance collapse — a systematic failure mode of per-component pipelines on stationary differenced series — is avoided.

On the held-out test set of 237 observations, the proposed model outperformed all benchmarks across the primary evaluation metrics. Relative to ARIMA(2,0,2), the strongest linear baseline, the framework improved Pearson correlation by 0.662 (from 0.150 to 0.812) and reduced mean absolute error by 37.3% (from 0.249 to 0.156 Baht/kg/day). A secondary but methodologically important finding concerns the inadequacy of Directional Accuracy as a standalone evaluation criterion for differenced commodity price series. The Vanilla LSTM baseline recorded a nominally superior DA of 78.6%, yet confusion matrix and variance ratio diagnostics revealed that this arose from systematic positive directional bias — the model correctly identified only one price decrease in 237 test observations (StdR $=0.033$). This finding argues for the routine inclusion of variance-sensitive metrics, specifically Pearson correlation and Standard Deviation Ratio, alongside directional measures in commodity price forecasting evaluations. This recommendation is consistent with a growing methodological consensus: Bürgi [51] demonstrates that DA lacks a self-sufficient evaluation framework and requires supplementation by magnitude-sensitive and dependence-robust procedures, while Costantini et al. [53] show that DA becomes materially informative only when paired with directional value measures that weight correct calls by the size of the realized move.

The VMD decomposition further revealed the multi-scale structure underlying RSS3 price dynamics. The dominant trend component (IMF₁, 78.6% of signal energy) reflects long-run supply–demand fundamentals, while intermediate components (IMF₂–IMF₃, 11.9%) correspond to inventory and export cycles, and high-frequency residuals (IMF₄–IMF₅, 9.4%) capture speculative noise and short-horizon information shocks. The primary forecasting gains were attributable to the lower-frequency components, which provided stable sequential context to the BiLSTM pathway — consistent with prior rubber-futures evidence that decomposition improves performance when retained modes are matched to the temporal scale of the prediction target [7]. The precise interaction mechanisms between IMF characteristics and the dual-pathway encoder remain an open question warranting systematic ablation across decomposition modes in future work.

Three limitations constrain the current framework. First, eight of the 24 intended input features — including STR20, Latex, and CupLump FOB prices — were unavailable at the time of analysis and substituted with zeros, likely suppressing medium-frequency predictive accuracy. Second, multi-step forecasts beyond one day hold exogenous variables constant at their last observed values; this persistence assumption becomes increasingly untenable as the forecast horizon extends, which is reflected in the rapid degradation of Pearson correlation beyond $h=5$. Third, the VMD decomposition is fitted statically on the training partition without adaptive refitting, potentially reducing alignment during structural market transitions. Future work will address these constraints through dataset completion, rolling adaptive decomposition, and an extension toward probabilistic forecast intervals to support price risk management across the Thai natural rubber supply chain.

Policy Implication

The findings carry direct implications for price risk management along the Thai rubber supply chain. A model that correctly anticipates the direction of daily RSS3 FOB price movements approximately two-thirds of the time — while reproducing the distributional characteristics of observed price changes — provides a statistically grounded basis for procurement timing, hedge initiation, and inventory holding decisions by processors, exporters, and institutional traders. More broadly, the demonstration that variance collapse is a systematic failure mode of conventional decomposition pipelines when applied to differenced commodity series has direct relevance for the design of early-warning systems used by national rubber price stabilization programmes, where near-zero forecast variance can render model signals operationally uninformative. Policymakers and market regulators should therefore require that forecasting tools submitted for operational deployment report variance-sensitive diagnostics — specifically StdR and Pearson correlation in the differenced space — alongside the directional accuracy statistics that are currently standard in industry practice.This is not merely a matter of reporting convention: Costantini and Kunst [54] demonstrate that selecting forecasting models on the basis of DA alone seldom yields robust gains and can worsen performance on magnitude-sensitive criteria — implying that operational deployment decisions grounded solely in hit-rate statistics risk systematically selecting the wrong model. The proposed evaluation framework is immediately adoptable without additional data requirements and is applicable across agricultural commodity markets more broadly.

As natural rubber markets become increasingly integrated with global financial systems and sensitive to climate-induced supply disruptions, the capacity to generate reliable, variance-faithful short-horizon price forecasts will grow in operational importance; the present framework offers a methodologically transparent and empirically validated foundation on which such systems can be built and extended.