A Comparative Study of Univariate Models for Baltic Dry Index Forecasting

1. Introduction

The Baltic Dry Index (BDI) serves as a critical barometer for the global shipping industry, reflecting the costs of transporting dry bulk commodities such as coal, iron ore, and grain. As an economic indicator, the BDI provides valuable insights into the health of global trade, supply chain dynamics, and overall economic activity. Given its relevance, accurate forecasting of the BDI is of paramount importance for a diverse range of stakeholders, including shipping companies, commodity traders, investors, policymakers, and economists.

Historically, the BDI has exhibited significant volatility, influenced by a myriad of factors. These factors include global economic cycles, seasonal demand fluctuations, geopolitical events, and changes in shipping capacity. For instance, economic downturns often lead to reduced demand for raw materials, subsequently lowering shipping rates and the BDI. Conversely, periods of economic growth drive up demand for commodities, resulting in higher shipping rates. Additionally, events such as trade disputes, natural disasters, and changes in regulatory environments can cause sudden and unpredictable movements in the BDI. This volatility poses substantial challenges for forecasting but also presents opportunities for leveraging advanced analytical techniques to improve prediction accuracy.

Traditional forecasting methods, such as time series analysis and econometric models, have been employed with varying degrees of success in predicting BDI movements. These methods provide a foundation for understanding past trends and making short-term forecasts, they often struggle to capture the complex, non-linear relationships inherent in the factors influencing the BDI.

In recent years, advancements in machine learning and artificial intelligence have opened new avenues for improving forecast accuracy. Techniques such as artificial neural networks (ANNs), support vector machines (SVMs), and ensemble learning methods offer the potential to model complex patterns and interactions that traditional methods may overlook. These contemporary approaches can incorporate a wider range of variables and are capable of learning from large datasets, making them well-suited for forecasting highly volatile indices like the BDI.

Forecasting the Baltic Dry Index (BDI) is crucial for the shipping industry as it offers insights into global trade trends, aiding in fleet management, capacity planning, and route optimization. Accurate BDI forecasts help shipping companies reduce costs, improve efficiency, and make informed investment decisions. Additionally, understanding shipping demand trends aids in predicting commodity prices and managing risks, ultimately enhancing the competitiveness and operational efficiency of shipping companies. Thus, forecasting BDI has attracted considerable attention from both academic researchers and industry practitioners. The primary objective of this study is to compare various univariate forecasting techniques in order to develop a more accurate short-term prediction model for the BDI, offering practical insights for the shipping sector.

The structure of this paper is as follows: Section 2 reviews relevant literature, Section 3 details the research methodology, Section 4 analyzes and compares the forecasting results, and Section 5 concludes with key findings and implications.

2. Literature Review

BDI forecasting methods can be broadly categorized into three groups: conventional econometric models, nonlinear models and machine learning techniques, and ensemble machine learning methods. These approaches have become the most widely adopted methodologies for BDI analysis in academic research.

2.1. Econometric Model

Econometric models were among the first developed for BDI forecasting, and they have been widely employed in validation studies. Both univariate and multivariate econometric methodologies are commonly applied in this research. Previous studies frequently utilized models such as the vector error correction model (VECM), vector autoregressive (VAR), and autoregressive integrated moving average (ARIMA). For instance, Veenstra and Franses [1] developed a time series-based VAR model to forecast the BDI, identifying non-stationarity in the data, which raised concerns about the stability of various modeling techniques. Cullinane et al. [2] used the ARIMA model for modifications to the Baltic Freight Index (BFI), demonstrating its superior predictive power. Kavussanos and Alizadeh [3] examined the seasonality of dry bulk freight rates, employing seasonal ARIMA and VAR models to assess fluctuations across different vessel sizes. Tsioumas et al. [4] introduced a multivariate vector autoregressive model with exogenous variables (VARX), which outperformed ARIMA in improving BDI forecasting accuracy. Zhang et al. [5] compared ARIMA, ARIMA-GARCH models, and artificial neural networks (ANNs), noting that econometric models performed better in one-step-ahead forecasts, while ANN-based algorithms achieved a higher direction matching rate in weekly and monthly data. More recently, Katris et al. [6] proposed data-driven model selection using ARIMA, fractional ARIMA (FARIMA), and ARIMA/FARIMA models with GARCH and EGARCH errors, marking the first use of FARIMA and MARS models in BDI forecasting.

2.2. Nonlinear Models and Machine Learning

Traditional econometric models often struggle to effectively capture the nonstationary and nonlinear characteristics of BDI time series. As a result, recent years have seen the increasing use of nonlinear regression models and machine learning techniques in BDI prediction research. Common methods include artificial neural networks, support vector machines (SVM), and nonlinear regression. SVMs are praised for their ability to approximate nonlinear functions and generalize well. For instance, Yang et al. [7] showed that SVM models can successfully account for BDI’s nonlinear properties and developed a freight early warning system based on SVM to predict fluctuations in shipping market prices.

Similarly, Li and Parsons [8] employed neural networks to forecast short-and long-term monthly tanker freight rates. Their results demonstrated that neural networks outperformed traditional ARIMA models, particularly for long-term forecasts. Building on this, Sahin et al. [9] introduced three ANN techniques for predicting BDI, conducting detailed comparisons to identify the most effective model. In recent years, advanced models like convolutional neural networks (CNNs) and long short-term memory (LSTM) networks have gained recognition as top-performing predictive neural networks, though their use in BDI forecasting remains underexplored.

Additionally, both econometric models and AI methods have their individual limitations. Hybrid forecasting approaches have emerged as a way to combine these models, leveraging their strengths to offset individual weaknesses. Such integrated systems tend to offer higher accuracy. For example, Chou and Lin [10] developed a fuzzy neural network combined with technical indicators to predict the BDI, achieving better accuracy than either technique alone.

Most recently, Liu et al. [11] proposed a BiLSTM-based system for BDI forecasting, incorporating the gray relational degree analysis method to select seven highly correlated factors as inputs. This model outperformed common machine learning models like SVR and REG, as well as the standard LSTM neural network.

2.3. Ensemble Machine Learning Models

Ensemble learning has been extensively employed to enhance model performance, serving as an effective method for boosting the predictive power of individual models. Leonov and Nikolov [12] found that freight pricing in shipping markets is highly volatile. They analyzed fluctuations in the Baltic Panamax 2A and 3A freight rates using a novel approach in shipping economics: a hybrid wavelet-neural network model. In a study on Baltic Dry Index prediction, Zeng et al. [13] applied Empirical Mode Decomposition (EMD) to break down the original BDI series into several intrinsic mode functions (IMFs). Each IMF and combined component was modeled using an Artificial Neural Network, yielding a forecasting method based on EMD and ANN. Their results showed that this approach outperformed other techniques such as ANN and VAR.

Similarly, Kamal et al. [14] developed a deep learning ensemble model by integrating Recurrent Neural Networks (RNN), Long Short-Term Memory, and Gated Recurrent Units (GRU) to improve BDI forecasting. The ensemble method outperformed both a single econometric model and a single machine learning model, highlighting its superior predictive capabilities.

Most recently, Su et al. [15] predicted the Baltic Dry Index using a deep integrated model (CNN-BiLSTM-AM) comprising a convolutional neural network, bi-directional long short-term memory (BiLSTM), and the attention mechanism (AM). Their findings indicate that the integrated model CNN-BiLSTM-AM encompasses the nonlinear and nonstationary characteristics of the shipping industry, and it has a greater prediction accuracy than any single model.

Univariate forecasting is simple, efficient, and easy to interpret as it relies solely on past values of the variable being predicted, requiring less data and computational resources. It is robust against noise from external factors, and suitable for quick, short-term forecasts. Additionally, it serves as a strong baseline for comparing more complex models. Based on the above literature, this study proposes a comparison of univariate models, including GM(1,1), ARIMA, SVR, and EMD-SVR-GWO.

3. Methods

3.1. Grey Forecasting Model

Grey theory, initially proposed by Deng [16], is designed to analyze systems that have incomplete or limited information. A key advantage of grey models is their ability to generate forecasts using relatively small datasets.

The GM(1,1) model refers to a first-order, single-variable grey model. The procedure for constructing GM(1,1) is outlined as follows:

Let the original data sequence be

$$x^{(0)}=\left(x^{(0)}(1),x^{(0)}(2),x^{(0)}(3),\mathrm{...},x^{(0)}(n),\right)$$

(1)

where x⁽⁰⁾(k) stands for the original data sequence in period k. The following sequence $x^{(1)}$ is defined as

$$x^{(1)}=\left({\displaystyle \sum _{k=1}^1{x}^{(0)}}(k),{\displaystyle \sum _{k=1}^2{x}^{(0)}(k),\mathrm{...},{\displaystyle \sum _{k=1}^n{x}^{(0)}}(k)}\right)=\left(x^{(1)}(1),x^{(1)}(2),x^{(1)}(3),...,x^{(1)}(n).\right)$$

(2)

Equation (2) was generated based on the accumulated generating operation (AGO) of Equation (1).

Before building the model, a class ratio must be performed to ensure that the sequence is suitable for modeling. The class ratio$\sigma (k)$is defined as follows:

$${\sigma }^{(1)}(k)={\displaystyle \frac{x^{(1)}(k-1)}{x^{(1)}(k)}}, k\ge 2,$$

(3)

If ${\sigma }^{(1)}(k)$$\in$(0,1), then $x^{(1)}(k)$ is appropriate for modeling.

After finishing the class ratio test, the GM(1, 1) model is formulated as a first-order differential equation for $x^{(1)}(k)$ as:

$${\displaystyle \frac{d{x}^{(1)}}{dk}}+a{x}^{(1)}=b$$

(4)

where a and b are coefficients to be estimated.

Next, we applied the ordinary least squares method to Equation (4) to estimate the coefficients of a and b. Once $\widehat{a}$ and $\widehat{b}$ are estimated, we generate predictions by substituting $\widehat{a}$ and $\widehat{b}$ into the following equation:

$${\widehat{x}}^{(1)}(k+1)=\left(x^{(0)}(1)-{\displaystyle \frac{\widehat{b}}{\widehat{a}}}\right)e^{-\widehat{a}k}+{\displaystyle \frac{\widehat{b}}{\widehat{a}}}\mathrm{and}{\widehat{x}}^{(1)}(1)=x^{(0)}(1)$$

(5)

Finally, the forecasted values of the time series are obtained by applying the inverse accumulated generating operation (IAGO) to convert to$x^{(0)}(k)$as follows:

$${\widehat{x}}^{(0)}(k+1)={\widehat{x}}^{(1)}(k+1)-{\widehat{x}}^{(1)}(k),$$

(6)

3.2. SARIMA

The seasonal model employed in this study belongs to the general class of univariate models originally introduced by Box and Jenkins [17]. SARIMA has been widely applied across various disciplines and remains a cornerstone in time series forecasting. A fundamental requirement in constructing this model is that the time series must be stationary, meaning its probabilistic properties remain constant over time.

SARIMA extends the ARIMA framework by incorporating seasonal components, making it suitable for time series data exhibiting seasonality. The model is typically expressed as SARIMA(p, d, q)(P, D, Q)s, where

p represents the order of the non-seasonal autoregressive terms,

d denotes the degree of non-seasonal differencing,

q corresponds to the order of the non-seasonal moving average terms,

P, D, and Q represent the seasonal autoregressive, seasonal differencing, and seasonal moving average orders, respectively,

s indicates the length of the seasonal cycle.

The general representation of the SARIMA model can be written as:

$${\phi }_p(B){\Phi }_p(B^s){\nabla }^d{\nabla }_s^D{Z}_t={\theta }_q(B){\Theta }_Q(B^s){\epsilon }_t,$$

(7)

where

ψ_p(B) = (1－ψ₁B－ψ₂B²－…－ψ_p B^p),

Φ_p(B^s) is the seasonal operator of order P,

B is the backshift operator with B^m(Z_t)=Z_t-m,

s is the season length,

▽^d=(1－B) is the non-seasonal operator,

${\nabla }_s^D={(1-B^s)}^D$ is the seasonal differencing operator,

Z_t is the stationary data at time t,

θ_q(B) = (1－θ₁B－θ₂B²－…－θ_qB^q),

${\Theta }_Q(B^s)$ is the seasonal operators of order Q, and ε_t is the white noise with zero mean and variance.

The model-building process involves three main stages: identification, parameter estimation, and diagnostic checking. In the identification phase, the tentative structure of the AR and MA terms (p and q) is determined by analyzing the autocorrelation and partial autocorrelation patterns. Once an initial specification is chosen, parameter estimation is carried out, often using methods such as maximum likelihood or least squares. Finally, diagnostic checking assesses the adequacy of the fitted model by examining residuals to ensure they exhibit the properties of white noise. If the diagnostics confirm model suitability, the SARIMA model is then employed for forecasting.

3.3. The EMD-SVR-GWO Model

3.3.1. EMD Model

Empirical Mode Decomposition is a data-driven technique developed by Norden E. Huang [18] in the late 1990s for analyzing non-linear and non-stationary time series data. Unlike traditional signal processing methods that rely on predefined basis functions (e.g., Fourier or wavelet transforms), EMD adaptively decomposes a signal into a set of intrinsic mode functions based on the inherent characteristics of the data. This makes EMD particularly useful for a wide range of applications, including signal processing, economics, geophysics, and biomedical engineering.

The EMD process consists of the following steps:

(1) identify Extrema: Detect all local maxima and minima in the signal. These points indicate where the signal changes direction and are used to construct the upper and lower envelopes.

(2) Interpolate: Using the identified extrema, create smooth upper and lower envelopes by interpolating between the maxima and minima, respectively. Common interpolation methods include cubic splines or piecewise linear functions.

(3) Mean Calculation: Compute the mean of the upper and lower envelopes. This mean represents the local trend of the signal.

(4) Sifting Process: Subtract the mean from the original signal to obtain a candidate IMF. This component is called a “proto-IMF.” The sifting process involves repeating the steps of extrema identification, interpolation, and mean calculation on the proto-IMF until it meets the IMF criteria: (i) The number of zero crossings and extrema must either be equal or differ by at most one. (ii) The mean value of the envelopes should be zero at every point.

(5) Residual Calculation: Once a valid IMF is obtained, subtract it from the original signal to get a residual. This residual represents the remaining signal after extracting the first IMF.

(6) Repeat: Apply the EMD process recursively to the residual to extract further IMFs until the residual becomes a monotonic function or has no more than two extrema. The decomposition is complete when the residual is either a constant, a monotonic trend, or contains no further oscillatory modes.

EMD has been applied across various fields due to its ability to handle non-linear and non-stationary data. In signal processing, it facilitates denoising, feature extraction, and fault detection in mechanical and electrical systems. The technique supports the analysis of financial time series, market trends, and economic cycles in economics. Geophysics benefits from EMD through seismic signal analysis, ocean wave modeling, and climate data interpretation. In biomedical engineering, EMD is instrumental in processing physiological signals such as electroencephalograms, electrocardiograms, and speech signals.

The advantages of EMD include its data-driven nature, which eliminates the need for predefined basis functions and makes it versatile for various types of data. Additionally, its adaptability allows it to handle non-linear and non-stationary signals effectively.

In summary, Empirical Mode Decomposition is a powerful and adaptive method for analyzing complex signals. By decomposing signals into intrinsic mode functions, EMD provides a flexible approach for revealing underlying patterns and trends in non-linear and non-stationary data. Its wide range of applications and ability to handle diverse data types make it an essential tool in various scientific and engineering disciplines.

3.3.2. SVR

Support Vector Regression (SVR) is a robust and versatile machine learning algorithm used for regression tasks. Originating from the principles of Support Vector Machines, SVR aims to predict continuous outcomes based on input features. While SVM is widely recognized for its effectiveness in classification problems, SVR extends these capabilities to the realm of regression, making it a valuable tool in various applications, such as financial forecasting, time series prediction, and biological data analysis.

At its core, SVR is designed to find a function that approximates the relationship between input features and a continuous target variable. The primary objective of SVR is to minimize the error of predictions while maintaining a degree of robustness by considering only a subset of the training data, known as support vectors. These support vectors are the critical elements of the training set that define the model.

SVR operates by introducing a margin of tolerance (epsilon), within which errors are ignored. This margin is referred to as the epsilon-insensitive zone. The essence of SVR is to ensure that the predictions lie within this margin for as many data points as possible, thereby balancing the complexity of the model with its predictive accuracy.

The SVR problem can be formulated as follows:

Given a training dataset ${\left\{(x_i,y_i)\right\}}_{i=1}^n$, where x_i represents the input features and y_i the target variable, the goal is to find a function f(x) that deviates from the actual targets y_i by at most ϵ. The function f(x) is typically expressed as:

f(x)=⟨w,x⟩+b where w is the weight vector and b is the bias term. The optimization problem is defined as (Vapnik [19]):

$${\mathrm{m}\mathrm{i}\mathrm{n}}_{w,b}{{\displaystyle \frac{1}{2}}‖w‖}^2+C\sum _{i=1}^n\left({\xi }_i+{\xi }_i^{*}\right)$$

subject to

$$\begin{gathered} y_i-⟨\mathrm{w},{\mathrm{x}}_i⟩-\mathrm{b}\le \mathsf{ϵ}+{\mathsf{\xi }}_i \\ ⟨\mathrm{w},{\mathrm{x}}_i⟩+\mathrm{b}-y_i\le \mathsf{ϵ}+{\xi }_i^{*} \\ {\xi }_i,{\xi }_i^{*}\ge 0 \end{gathered}$$

Here, ${\xi }_i,\mathrm{and}{\xi }_i^{\mathrm{*}}\mathrm{are}$slack variables that allow for deviations greater than ϵ, and C is a regularization parameter that determines the trade-off between the flatness of f(x) and the amount by which deviations larger than ϵ are tolerated.

One of the strengths of SVR, inherited from SVM, is the ability to handle non-linear relationships through the kernel trick. By mapping the input features into a higher-dimensional space using a kernel function, SVR can fit complex, non-linear functions. Commonly used kernels include the linear, polynomial, and radial basis function (RBF) kernels. In this research, the radial basis function is adopted.

3.3.3. GWO

The Grey Wolf Optimizer (GWO) is a metaheuristic algorithm inspired by the leadership structure and cooperative hunting strategies of grey wolves. These animals typically live in packs of 5–12 members and exhibit a well-defined social hierarchy. This hierarchy is the foundation of GWO’s design.

(1)

Social Hierarchy of Grey Wolves

In a grey wolf pack, leadership and responsibilities are distributed among four main roles:

Alpha (α): The alpha wolves, usually a dominant male and female, lead the pack by making critical decisions regarding hunting, resting, and movement. Interestingly, alpha wolves are not always the strongest physically; their primary role is managing and maintaining order. While their decisions are generally followed by the rest of the pack, occasional democratic behaviors are observed where alpha wolves heed suggestions from others. Submission to alpha decisions is typically shown by the other wolves lowering their tails.

Beta (β): Positioned just below alpha wolves, the beta wolves act as advisors and enforcers. They help implement alpha’s decisions, maintain discipline, and provide feedback. Betas are also the likely successors if an alpha dies or becomes incapable of leading.

Omega (ω): At the lowest rank, omega wolves serve as scapegoats and stress relievers for the group, which prevents internal conflicts. Although they eat last and hold minimal authority, omegas are essential for pack stability. Occasionally, they also serve as caretakers for pups.

Delta (δ): Wolves that do not belong to the above categories are called subordinates (delta wolves). They serve various roles such as scouts (patrolling territory), sentinels (guarding), hunters (assisting in hunting), and caregivers for weak or injured members.

Group Hunting Behavior

Grey wolves hunt cooperatively, and their hunting process can be divided into three stages:

i.

Tracking and approaching the prey.

ii.

Pursuing, encircling, and harassing until the prey is immobilized.

iii.

Attacking once the prey is exhausted.

GWO incorporates these two behavioral aspects—social hierarchy and hunting strategy—into its optimization mechanism (Muro et al. [20]).

(2)

Mathematical Modeling in GWO

• Hierarchy Representation

The algorithm assigns the best candidate solution as alpha (α), the second-best as beta (β), and the third-best as delta (δ). All other solutions are considered omega (ω). These three top solutions guide the remaining wolves during the optimization process.

2.

Encircling Prey

Grey wolves encircle their prey during hunting, which is modeled as follows:

$$\overrightarrow{D}=\left|\overrightarrow{C}·\overrightarrow{X_p}\left(t\right)-\overrightarrow{X}\left(t\right)\right|$$

$$\overrightarrow{X}\left(t+1\right)=\overrightarrow{X_p}\left(t\right)-\overrightarrow{A}·\overrightarrow{D}$$

where:

t = current iteration,

$\overrightarrow{X}$= position vector of the wolf,

$\overrightarrow{X_p}$ = position vector of the prey,

$\overrightarrow{A}$ and $\overrightarrow{C}$ = coefficient vectors,

$\overrightarrow{D}$ = vector used to specify a new position of the grey wolf.

The coefficient vectors are defined as (Mirjalili et al., [21]):

$$\overrightarrow{A}=2\overrightarrow{a}·\overrightarrow{r_1}-\overrightarrow{a}$$

$$\overrightarrow{C}=2\overrightarrow{r_2}$$

Here:

$\overrightarrow{a}$ a decreases linearly from 2 to 0 over iterations,

r1 and r2 are random vectors in [0, 1].

3.

Hunting (Exploitation Phase)

To simulate hunting, GWO assumes α, β, and δ wolves have the most accurate information about the prey’s position. The algorithm stores these three best solutions and updates others accordingly:

$${\overrightarrow{D}}_{\alpha }=\left|{\overrightarrow{C}}_1·\overrightarrow{X_{\alpha }}-\overrightarrow{X}\right|,{\overrightarrow{D}}_{\beta }=\left|{\overrightarrow{C}}_2·\overrightarrow{X_{\beta }}-\overrightarrow{X}\right|,{\overrightarrow{D}}_{\delta }=\left|{\overrightarrow{C}}_3·\overrightarrow{X_{\delta }}-\overrightarrow{X}\right|$$

$$\overrightarrow{X_1}=\overrightarrow{X_{\alpha }}-{\overrightarrow{A}}_1·\left({\overrightarrow{D}}_{\alpha }\right),\overrightarrow{X_2}=\overrightarrow{X_{\beta }}-{\overrightarrow{A}}_2·\left({\overrightarrow{D}}_{\beta }\right),\overrightarrow{X_3}=\overrightarrow{X_{\delta }}-{\overrightarrow{A}}_3·\left({\overrightarrow{D}}_{\delta }\right)$$

$$\overrightarrow{X}\left(t+1\right)={\displaystyle \frac{\overrightarrow{X_1}+\overrightarrow{X_2}+\overrightarrow{X_3}}{3}}$$

This mechanism ensures wolves move toward the prey’s estimated location.

4.

Attacking the Prey

The attack occurs as $\overrightarrow{a}$ decreases from 2 to 0. When $\overrightarrow{a}$ lies within [-1, 1], wolves close in on the prey, leading to exploitation and convergence.

5.

Searching for Prey (Exploration)

During exploration, wolves diverge to search for prey. This is controlled by $\overrightarrow{A}$; if ∣A∣> 1, wolves move away from the prey to explore new areas. Another component, $\overrightarrow{C}$, introduces randomness by adjusting the emphasis on the prey’s position (C > 1 amplifies attraction; C < 1 reduces it). Unlike $\overrightarrow{A}$, $\overrightarrow{C}$does not decrease linearly, which helps maintain global search capability and prevents premature convergence. Overall, these mechanisms allow GWO to balance exploration and exploitation, reducing the risk of getting trapped in local optima.

Figure 1 presents the flowchart of the GWO algorithm.

3.3.4. Model Construction of EMD-SVR-GWO

The EMD-SVR-GWO model is a hybrid forecasting framework designed to improve prediction accuracy by leveraging the strengths of three advanced techniques: Empirical Mode Decomposition, Support Vector Regression, and the Grey Wolf Optimizer. This model is particularly effective for complex, nonlinear, and non-stationary time series data, such as financial indices, economic indicators, or environmental data. The construction of the EMD-SVR-GWO model involves three primary stages:

1. Empirical Mode Decomposition

EMD is a data-driven decomposition technique that breaks down the original data into a set of components called Intrinsic Mode Functions and a residual trend. Each IMF represents oscillatory modes embedded within the data, capturing different frequency characteristics, while the residual captures the long-term trend.

The decomposition enables the isolation of distinct patterns (high-frequency fluctuations, medium-term cycles, and long-term trends), which are modeled separately in the next stage.

2. Support Vector Regression

After decomposition, each IMF and the residual are treated as independent sub-series. Support Vector Regression is applied to predict each component individually. SVR is chosen for its robustness in handling nonlinear relationships and its strong generalization capabilities.

Each IMF captures different dynamics of the original time series, and SVR models are trained separately to forecast these dynamics effectively.

3. Grey Wolf Optimizer

The final stage involves the reconstruction of the original time series by combining the predicted IMFs and residual. To achieve the most accurate forecast, the Grey Wolf Optimizer is employed to determine the optimal combination of these forecasts. GWO is a metaheuristic optimization algorithm inspired by the hunting behavior and social hierarchy of grey wolves in nature.

The EMD-SVR-GWO hybrid model offers a powerful and flexible framework for time series forecasting. EMD enhances the model by decomposing complex data into simpler components, SVR captures nonlinear relationships within each component, and GWO optimizes the final forecast by adjusting the weights for the most accurate reconstruction. This approach is highly effective for datasets with mixed trends, seasonality, and noise, providing superior forecasting performance compared to traditional single-model methods.

4. Data and Results

This section begins with a description of the dataset used in the analysis, followed by the presentation of results from all forecasting models. Subsequently, the models’ performance is assessed using three evaluation criteria, and their forecasting accuracy is compared.

4.1. BDI Time Series Data

The dataset consists of monthly average daily BDI values spanning January 2019 to August 2024. For the purpose of analysis, the data is divided into two segments: the in-sample set, covering January 2019 to December 2023, used for model estimation, and the out-of-sample set, covering January to August 2024, used for prediction. A summary of the monthly BDI data is provided in Table 1.

4.2. Model Results

This section summarizes the outcomes produced by the four models, while their forecasting accuracy is compared in the following section.

4.2.1. Grey Forecasting

The accuracy of forecasts generated by the grey model depends on the length of the initial data segment selected from the time series. To address this, we executed multiple grey forecasts using different initial sequence sizes. The best results, indicated by the lowest prediction errors, occurred when the initial sequence length was set to four observations. The procedure for producing forecasts with the grey model involves the following steps:

• Class Ratio Test:

Before applying the GM(1,1) model, the data must satisfy the class ratio condition. As described in Section 3.1, the ratio values $\sigma (k)$ should lie between 0 and 1 for the series to qualify for grey modeling. The class ratio test results, summarized in Table 2, confirm that the grey model is suitable for this dataset.

2.

Accumulated Generating Operation (AGO):

Following Equation (2), the AGO was applied to obtain the new sequence. The resulting values are presented in Table 3.

3.

Mean value generating sequence

We calculated the mean value generating sequence, as shown in Table 4.

4.

Time series prediction model

The coefficients a and b were estimated using the least squares method, yielding the following values:

a=−0.1749295723 and b=1313.9202714553.

These estimates are used to get

$${{\displaystyle \widehat{X}}}^{(0)}\left(k+1\right)=\left(1-{{\displaystyle e}}^{\widehat{a}}\right)\left[{{\displaystyle X}}^{(0)}\left(1\right)-{\displaystyle \frac{\widehat{b}}{\widehat{a}}}\right]{{\displaystyle e}}^{-\widehat{a}k}$$

(8)

These parameters are then applied to Equation (8) to compute the predicted values of the time series. The forecasts for the period from January 2024 to August 2024 are presented in Table 8.

4.2.2. ARIMA Model

An examination of the historical data in Figure 2, along with the SAS output shown in Figure 3, indicates the presence of a trend in the series. To confirm stationarity, we performed an ADF test, which produced a statistic of -2.377 and a p-value of 0.1483. Since the p-value exceeds 0.05, the series is non-stationary. Therefore, differencing is required to achieve stationarity. After applying the first difference, the data became stationary, as illustrated in Figure 4 and Figure 5.

After converting the original data into a stationary series, the ACF and PACF plots were utilized to identify potential Box–Jenkins models. Three candidate models are summarized in Table 5. To evaluate model adequacy, the residuals from each model were analyzed using the Ljung–Box test, based on the SAS output.

Table 5 reports the Ljung–Box statistic and corresponding p-values for each candidate model. According to these results, Models 1 and 2 fail to meet adequacy criteria because not all p-values exceed 0.05. In contrast, Model 3 satisfies the adequacy requirement, as its Ljung–Box statistics (Chi-Square values) are relatively small and all associated p-values are greater than 0.05.

The best model identified is ARIMA(0,1,1), which is used to predict BDI. However, it achieves only good performance rather than excellent.

To further improve the model selection process, we use the auto_arima function to automatically identify the optimal parameters for BDI forecasting. There is an auto_arima function from the pmdarima Python package that automates the process of fitting ARIMA models to time series data. It systematically evaluates a range of models by varying the autoregressive (p), differencing (d), and moving average (q) parameters—and even seasonal components if applicable. To select the best model, auto_arima uses information criteria such as the Akaike Information Criterion (AIC), and the Bayesian Information Criterion (BIC), which balance the model’s goodness-of-fit against its complexity. By employing a stepwise search strategy, it efficiently narrows down the candidate models and ultimately selects the one with the lowest AIC (or the chosen metric), ensuring an optimal trade-off between fit and parsimony for forecasting purposes. The results produced by auto_arima are presented in Table 6.

The auto_arima procedure selected an ARIMA(1,0,0) model, which is essentially an AR(1) model. Despite its simplicity, this model delivered better forecast performance, as indicated by its lower AIC compared to more complex alternatives. The model’s parsimony helps prevent overfitting—a common concern when working with limited data—while still capturing the essential dynamics of the series.

In summary, although manual diagnostics like ACF/PACF plots and the ADF test suggest nonstationarity, auto_arima selects a simpler model (e.g., ARIMA(1,0,0)) that yields better forecasting performance. This happens because auto_arima optimizes forecast accuracy using criteria like AIC, which in small samples can favor avoiding unnecessary differencing.

Thus, while manual analysis is crucial for understanding the data’s properties, the empirical success of auto_arima shows that a simpler model can be more effective in practice. This underscores the need to balance theoretical diagnostics with practical performance to choose a model that not only fits the data well but also delivers robust predictions. In this research, the model suggested by auto_arima is adopted and the predicted values of the time series from January 2024 to August 2024 are calculated and presented in Table 8.

4.2.3. SVR

Support Vector Regression is a powerful machine learning technique for predicting the Baltic Dry Index. In this study, a total of 80 monthly BDI data points are used. The first 72 consecutive months serve as input features for the model, while the last 8 months are used as output targets.

To prevent training errors caused by high-dimensional data or large variations in feature values, the dataset must be normalized before training the SVR model. The standardization process follows the formula:

$$X^{\prime }=\frac{X-\mu }{\sigma }$$

where

X is the original feature value,

μ is the mean of the feature in the training set,

σ is the standard deviation of the feature in the training set,

X^′ is the standardized value.

Next, the SVR model is trained using an appropriate kernel function, which determines how the model captures non-linear relationships in the data. The three most common kernels are linear, polynomial, and radial basis function. For BDI prediction, the RBF kernel is selected due to its effectiveness in handling complex, non-linear patterns.

Finally, the trained SVR model is deployed for forecasting future BDI values. It is used to generate one-step predictions for the last 8 months of the dataset. The predicted BDI values for the time series from January 2024 to August 2024 are presented in Table 8.

4.2.4. EMD-SVR-GWO

In the proposed EMD-SVR-GWO hybrid ensemble model, the first stage is to apply EMD to decompose the data of BDI into several IMFs. Too many IMFs may result in a poor final result due to the accumulation error of each IMF. This paper decomposes the original data into four IMFs to avoid aforementioned problems. Figure 6 present the decomposition results of the BDI via EMD, where all IMFs are listed from the highest frequency to the lowest frequency.

The SVR is applied to each IMFs in the second stage. The RBF kernel is selected for SVR. The trained SVR model is deployed for forecasting the next 8 periods of each IMF.

The third stage is to optimize the weight of each IMF using GWO based on least squares. It is worth noting that when using the least squares method as an evaluation criterion, actual data is required. However, in forecasting, actual data is unknown. In this study the predicted values obtained from simple regression are used as the actual data. The optimal weight for each IMF is shown in Table 7. After obtaining the weight of each IMF, multiply this weight by the predicted value of each IMF and sum them up to obtain the predicted value of BDI. The predicted BDI are shown in Table 8.

4.3. Comparison of Forecasting Results

The BDI forecasts for the out-of-sample period (January to August 2024) were generated using each forecasting approach. The predicted values, along with the actual BDI figures for comparison, are presented in Table 8.

Table 8. Actual and predicted BDI.

	Actual BDI	Grey forecast	ARIMA (1,0,0)	SVR	EMD-SVR-GWO
1	1673.61	2871.9526	2305.03	1402.1657	1903.5586
2	1658.16	1915.8460	2213.99	1441.9440	1902.6450
3	2232.90	1179.6161	2136.49	1484.2511	1901.7250
4	1731.33	2517.5462	2070.52	1527.7972	1900.8048
5	1890.36	1941.7500	2014.37	1571.3694	1899.8906
6	1922.00	1614.6354	1966.57	1613.8601	1898.9877
7	1925.30	2044.0939	1925.88	1654.2896	1898.1015
8	1716.24	1947.6592	1891.24	1691.8222	1897.2360

Table 8. Actual and predicted BDI.

	Actual BDI	Grey forecast	ARIMA (1,0,0)	SVR	EMD-SVR-GWO
1	1673.61	2871.9526	2305.03	1402.1657	1903.5586
2	1658.16	1915.8460	2213.99	1441.9440	1902.6450
3	2232.90	1179.6161	2136.49	1484.2511	1901.7250
4	1731.33	2517.5462	2070.52	1527.7972	1900.8048
5	1890.36	1941.7500	2014.37	1571.3694	1899.8906
6	1922.00	1614.6354	1966.57	1613.8601	1898.9877
7	1925.30	2044.0939	1925.88	1654.2896	1898.1015
8	1716.24	1947.6592	1891.24	1691.8222	1897.2360

Yokum and Armstrong [22] carried out two studies examining experts’ views on the criteria for selecting forecasting methods. Their findings indicated that accuracy is regarded as the most important factor by the majority of researchers. However, because no single accuracy metric is universally applicable to all forecasting scenarios, multiple measures are often employed to provide a comprehensive evaluation of model performance. As noted by Makridakis et al. [23], the ranking of models can vary depending on the metric applied.

In this study, three widely used accuracy measures are adopted: root mean squared error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). These measures are defined as follows:

$$RMSE=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^n{\left(Y_i-{\stackrel{\wedge }{Y}}_i\right)}^2}}{n}}}$$

(9)

$$MAE={\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left|Y_i-\stackrel{\wedge }{Y_i}\right|}}{n}}$$

(10)

$$MAPE={\displaystyle \frac{100{\displaystyle \sum _{i=1}^n\left|\frac{Y_i-\stackrel{\wedge }{Y_i}}{Y_i}\right|}}{n}}$$

(11)

where Y_i and ${\stackrel{\wedge }{Y}}_i$ denote the actual and predicted values of the time series for period i, respectively. All three performance metrics yield positive values, and a smaller value indicates a more accurate forecasting model.

The comparative evaluation of forecasting accuracy across all methods is summarized in Table 9. According to the results, the EMD-SVR-GWO model delivers the highest accuracy, as it consistently records the lowest values across all three metrics. ARIMA ranks as the second most effective model, regardless of which measure is applied. The difference in accuracy between SVR and ARIMA is relatively minor. In contrast, the Grey forecasting approach performs the poorest in predicting the monthly BDI.

A commonly cited reference for the MAPE-based classification is Lewis [24]. In this work, Lewis outlines that a MAPE below 10% indicates excellent forecast accuracy, 10%-20% is good, 20% -50% is acceptable, and values above 50% are considered poor.

According to Lewis’s classification for forecast accuracy, our proposed EMD-SVR-GWO model achieves a MAPE of 8.3455, placing it in the excellent category. In contrast, the ARIMA and SVR models record MAPEs of 14.2829 and 15.3686 respectively, which are classified as good. Although the Grey Forecast model has a MAPE of 27.26, falling into the acceptable range, its performance is noticeably inferior compared to the other models. These results highlight the superior predictive accuracy of our novel EMD-SVR-GWO model in forecasting BDI.

5. Conclusions

The Baltic Dry Index is a key indicator of global shipping costs and economic activity, influencing decisions in trade, investment, and policymaking. However, its high volatility makes accurate forecasting essential. Reliable predictions help shipping companies optimize operations, investors anticipate market trends, and policymakers make informed decisions, reducing financial uncertainty and improving strategic planning in the shipping and trade industries.

This paper compares various univariate forecasting methods to develop a more precise short-term BDI prediction model, providing valuable insights for decision-makers. Four forecasting techniques are examined: Grey Forecast, ARIMA, Support Vector Regression, and EMD-SVR-GWO. Our results indicate that EMD-SVR-GWO outperforms ARIMA and other two methods (SVR and Grey Forecast). Our proposed approach goes beyond simply combining EMD with SVR by introducing an additional composition step. Notably, this step plays a significant role in enhancing the overall forecasting accuracy. The EMD-SVR-GWO model achieves a MAPE of 8.3455, classified into the excellent category.

From an industry perspective, the novel EMD-SVR-GWO model provided by this study serves as a valuable reference for ship-owners and charterers in making chartering decisions. For instance, if the Baltic Dry Index is projected to rise, ship owners should consider purchasing new or second-hand vessels or securing time charter contracts. If they already hold long-term transportation contracts, sub-chartering their vessels may be a strategic move. On the other hand, charterers should act promptly to secure time charter contracts or long-term transportation agreements with ship-owners. Conversely, if the Baltic Dry Index is expected to decline, the opposite strategies should be adopted.

Future research should integrate the Grey Wolf Optimizer or other metaheuristic algorithms—such as Particle Swarm Optimization—to optimize the hyperparameters of the Support Vector Regression model, thereby further enhancing the robustness and accuracy of the BDI model. Furthermore, the current univariate approach considers only BDI data and overlooks variables like macroeconomic policies across countries and the role of investor sentiment. Future model enhancements will aim to integrate these factors for improved forecasting accuracy.