An Empirical Analysis of Sino-Russia Foreign Trade Turnover Time Series: Based on EMD-LSTM Model

Abstract

The time series of foreign trade turnover is complex and variable and contains linear and nonlinear information. This paper proposes preprocessing the dataset by the EMD algorithm and combining the linear prediction advantage of the SARIMA model with the nonlinear prediction advantage of the EMD-LSTM model to construct the SARIMA-EMD-LSTM hybrid model by the weight assignment method. The forecast performance of the single models is compared with that of the hybrid models by using MAPE and RMSE metrics. Furthermore, it is confirmed that the weight assignment approach can benefit from the hybrid models. The results show that the SARIMA model can capture the fluctuation pattern of the time series, but it cannot effectively predict the sudden drop in foreign trade turnover caused by special reasons and has the lowest accuracy in long-term forecasting. The EMD-LSTM model successfully resolves the hysteresis phenomenon and has the highest forecast accuracy of all models, with a MAPE of 7.4304%. Therefore, it can be effectively used to forecast the Sino-Russia foreign trade turnover time series post-epidemic. Hybrid models cannot take advantage of SARIMA linear and LSTM nonlinear forecasting, so weight assignment is not the best method to construct hybrid models.

1. Introduction

China and Russia have maintained close strategic collaboration in international and regional affairs, which has strongly maintained peace and stability in the region and the world. China has held on to its position as Russia’s top trading partner for 12 consecutive years, and Russia ranks 10^th among China’s major trading partners. Russia is a major energy exporter in the world, while China is a leading industrial power in the world, importing large amounts of energy from Russia. There is no doubt that energy is a key part of the bilateral trade and economic relations between China and Russia. Over the last 10 years, China’s GDP has ranked approximately 2^nd and Russia’s at 11^th all over the world. The economic sanction of European and American countries is the most direct factor limiting the growth of Russia’s GDP and foreign trade turnover.

To ease the stress on the Russian economy from the economic sanctions of Europe and American countries, the “Look East” strategy has been proposed and implemented, which has provided an opportunity to strengthen economic and trade cooperation with East Asian countries. So, Russia strengthened its economic and trade cooperation with China. Sino-Russia foreign trade turnover from 2013 to 2020 was $714.9 billion. The main export and import commodities are minerals (33% of total trade), machinery, and equipment (29% of total trade) (RU-STAT, 2021). Energy cooperation is always the most important, fruitful, and wide-ranging area of practical cooperation between China and Russia. Among Russian major trade partners, China is the most essential one, and the “Look East” strategy shows remarkable results. The Sino-Russia foreign trade turnover time series, which is the main external expression of the import and export markets, presents a very complex operating mechanism. At the same time, the time series of foreign trade turnover implicitly shows the current situation and future change trend of Sino-Russia foreign trade development. Therefore, it is particularly important to understand its fluctuation pattern to optimize the structure of import and export trade, increase foreign exchange reserves and create more employment in economic activities.

2. Literature Review

Traditional time series models include the grey prediction model (GM(1,1)), holt-winters (HW), auto-regressive (AR), moving average (MA), vector auto-regression (VAR), autoregressive moving average (ARMA) model, seasonal auto-regressive integrated moving average (SARIMA) models, etc. Scholars use traditional time series models to conduct in-depth studies of foreign trade turnover forecasts. For instance, Zhang et al. (2013) improved on the traditional GM(1,1) model to increase the forecast accuracy and finds that there is still much potential for China-Russia trade cooperation. The ARIMA model is a classical time series analysis method with high forecast accuracy in the short term and has been widely used in the field of foreign trade markets. The ARIMA model is used by Farooqi (2014) to forecast Pakistan’s total exports and imports from 1947 to 2013, and he finds that the ARIMA (2, 2, 2) model is suitable for forecasting annual imports and the ARIMA (1, 2, 2) model is suitable for forecasting annual exports.

The application of the ARIMA model greatly contributed to forecasting research in financial markets and itself became a classic control model for various new forecasting methods (Kim &Won, 2018; Suh et al., 2014). Mladenovic et al. (2016), used the seasonal Holt-Winters and ARIMA methods to analyze export trend changes for the Republic of Serbia. It is confirmed that the two methods are consistent in predicting export trend changes, which provides a valid reference method for the analysis of other macroeconomic indicators. In another example, Urrutia et al. (2019) used the ARIMA model and the Bayesian Artificial Neural Network (BANN) model to forecast the turnover of Philippine imports and exports, as well as to compare the predictive performance of the two models. They conclude that there is no significant difference between the actual and predicted values of Philippine imports and exports, so the BANN model is the most appropriate model for forecasting the turnover of Philippine imports and exports.

With the research of deep learning, LSTM neural networks alleviate the gradient disappearance and explosion phenomenon that easily occurs during the training process of traditional recurrent neural networks (RNN) (Gers et al., 2000), and the forecasting advantage in foreign trade time series is gradually highlighted. Many scholars have studied foreign trade turnover time series by using LSTM models. Qu et al. (2019) proposed forecasting the trend of trade turnover using a recurrent neural network (RNN) based on LSTM. The MSE exceeds that of the traditional time series model at 12.2% and proves the effectiveness of the LSTM-RNN model. The LSTM neural network is also applied in fields such as agricultural crop yield forecasting (Haider et al., 2019) and stock price forecasting (Nelson et al., 2017; Gyamerah & Korda, 2021), which shows superior predictive performance in comparison with other machine learning methods.

Because of the complexity and cyclical characteristics of foreign trade turnover time series, single models have gradually revealed shortcomings in long-term forecasting. To compensate for the shortcomings of single models, various combined models have been developed. Yu et al. (2022) combined the SARIMA model with the LSTM model and innovatively constructed the SALSTM model to forecast economic indicators related to the competitiveness of export products and proved its effectivity and superiority by the mean square error (MSE), which provides a reference basis for national trade policy formulation. In another study, the ARIMA-LSTM hybrid model is used by Dave et al. (2021) to forecast the turnover of Indonesian exports. In addition, the ARIMA-LSTM hybrid model is also applied to forecast stock prices (Alzheev & Kochkarov, 2020; Zainuri et al., 2021). Zhang and Sun (2019) examined the combined method of hybrid models and constructed the ARIMA-SVM hybrid model by the weight assignment method to forecast and analyze Chinese foreign trade turnover from 1980 to 2016. The result shows that the forecast accuracy of the hybrid model is higher than that of the single model.

Although the LSTM model demonstrates superiority in the forecasting of time series, because the series is not transformed from non-stationary to stationary, the forecast accuracy still needs to be improved. The empirical mode decomposition (EMD) method can decompose nonlinear and non-smooth data into smooth data, which is a good selection for data preprocessing. The original time series is decomposed by EMD into a fixed number of intrinsic mode functions (IMFs) at different frequencies, and then the decomposed signal is simulated using other models. In the end, all the simulated results are superimposed to improve the forecast accuracy of the hybrid models (Li et al, 2018; Solanki et al., 2020). As an example, Qin et al. (2018) use the hybrid model of integrated empirical modal decomposition (EEMD) and support vector regression (SVR) to forecast the Chinese consumer price index (CPI) and verify that the hybrid model has higher forecast accuracy than the EEMD single model. Shu and Gao (2020) use EMD and a convolutional neural network (CNN) to preprocess the original stock price series, which is modeled by the LSTM neural network and obtain the final predictive price after linear transformation. The hybrid network can achieve better performance by modeling different frequencies than other models. However, there are fewer studies that combine the SARIMA model with the EMD-LSTM model to forecast trade turnover time series.

To summarize the above studies, it can be found that the SARIMA model has advantages in short-term forecasting and that the LSTM neural network algorithm has more potential for application in the economic field. The hybrid model has a higher forecast accuracy than the single model. Therefore, this paper proposes the following hypotheses:

H1: The EMD algorithm can effectively resolve the hysteresis phenomenon generated by the LSTM model and improve the forecast accuracy;

H2: The weight assignment hybrid model combines the linear forecasting advantage of the SARIMA model with the nonlinear forecasting advantage of the EMD-LSTM model, so it has the highest forecast accuracy of all models;

H3: Weight assignment is an efficient method to construct the hybrid model of foreign trade turnover time series.

3. Data and Model Specification

3.1. Data Source and Assessment Method

This paper uses monthly historical data on Sino-Russia foreign trade turnover from 2013 to 2021 as a research sample; the December 2021 trade turnover is not yet available, so the 2021 average turnover is used instead. The foreign trade turnover from 2013–2019 is the training set, and the rest is the test set. The ratio of the training set to the test set is 7:2, and the predicted values of the test set are compared with the actual values for validation. The mean square error (MSE), mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE) are used as indicators for the evaluation of the model. The mathematical formula is:

\(\begin{aligned}\mathrm{MSE}=\frac{1}{n} \sum_{i=1}^{n}\left(\hat{y}_{i}-y_{i}\right)^{2}\\\end{aligned}\)       (1)

\(\begin{aligned}\mathrm{MAE}=\frac{1}{n} \sum_{i=1}^{n}\left|y_{i}-\hat{y}_{i}\right|\\\end{aligned}\)       (2)

\(\begin{aligned}\mathrm{RMSE}=\sqrt{\frac{1}{n} \sum_{i=1}^{n}\left(\hat{y}_{i}-y_{i}\right)^{2}}\\\end{aligned}\)       (3)

\(\begin{aligned}\mathrm{RMSE}=\frac{1}{n} \sum_{i=1}^{n}\left|\frac{y_{i}-\hat{y}_{i}}{y_{i}}\right|\\\end{aligned}\)       (4)

where y_i is the actual value, ŷ_i is the predicted value, and n is the sample size.

3.2. Relevant Models

3.2.1. The SARIMA Model

The SARIMA model is a combination of seasonal differencing and an ARIMA model for forecasting time series data with cyclical characteristics (John & Patrick, 2016). The model expression is:

SARIMA = (p, d, q) × (P, D, Q)s       (5)

where p is the auto-regressive orders, P is the seasonal auto-regressive orders, d is the differential counts, D is the seasonal differential counts, q is the moving average orders, Q is the seasonal moving-average orders, and s is the number of periods.

The minimum Akaike information criterion (AIC) (Akaike, 1974) can effectively compensate for the subjectivity of the fixed order of auto-correlation function (ACF) and partial auto-correlation function (PACF) plots and can find the optimal hyper-parameters more quickly within a limited range of orders. In this paper, the model parameters are determined by the AIC values, and the mathematical formula is:

\(\begin{aligned}\mathrm{AIC}=2 k-2 \ln (\hat{L})\\\end{aligned}\)       (6)

where k is the number of parameters and \(\begin{aligned}\hat{L}\end{aligned}\) is the maximum likelihood function.

3.2.2. The EMD Algorithm

The methods for transforming non-stationary time series into stationary time series include log transformation, exponential smoothing, difference, and decomposition. The EMD algorithm is a signal processing method that converts non-stationary and nonlinear data into stationary and linear data and is an important tool for exploring the hidden time series relationships in the data. It does not require a predetermined basis function and relies on the signal itself to adaptively decompose the IMF and residual terms for the time series (Li et al., 2019). The formula is:

X(t) = Σⁿ_i=0C_i(t) + R_n(t)       (7)

where X(t) is the time series, C_i is the IMF, and R_n and is the residual term, which usually represents the overall trend of the time series.

3.2.3. The LSTM Model

The LSTM model adds an input gate, output gate, and forgetting gate to the neurons in the RNN. The forgetting gate is mainly used to decide which information needs to be retained and which information needs to be updated, so the training results are reflected in the cell state and output information. This structure is a better solution to the problem of gradient disappearance or gradient explosion in RNNs during training (Figure 1).

Figure 1: The Structure of the LSTM Neural Network

The data first pass through the forget gate, which determines the forgotten and retained information in the cell state C_t−1 at moment t−1. The activation function σ maps [h_t−1, x_t] to 0 and 1, where 0 means “completely discard” and 1 means “completely retain” (Ni et al., 2017). The output value f_t of the activation function, σ determines whether [h_t−1, x_t] is forgotten or not. The mathematical expression for f_t is:

f_t = σ (W_f · [h_t−1, x_t] + b_f)       (8)

where h_t−1 is the output information at moment t−1, x_t is the input information at moment t, W_f is the weight of [h_t−1, x_t], b_f, is the bias amount, and σ is the activation function.

The data pass through the input gate, and the valid information contained in [h_t−1, x_t] is determined. The output value i_t of the activation function, σ determines the proportion of valid information. The amount of valid information in [h_t−1, x_t] is calculated using the function to generate a new vector of candidate values \(\begin{aligned}\tilde{c}_{t}\\\end{aligned}\) (Song et al., 2020). The mathematical expressions for i_t and \(\begin{aligned}\tilde{c}_{t}\\\end{aligned}\) are:

i_t = σ (W_t · [h_t−1, x_t] + b_i)       (9)

\(\begin{aligned}\tilde{c}_{t}=\tanh \left(W_{c} \cdot\left[h_{t-1}, x_{t}\right]+b_{c}\right)\\\end{aligned}\)       (10)

Update the cell state C_t−1 at moment t−1, and acquire the cell state C_t at moment t. The mathematical expression for C_t is:

\(\begin{aligned}C_{t}=f_{t} C_{t-1}+i_{t} \tilde{c}_{t}\\\end{aligned}\)       (11)

where f_tC_t−1 is the retained information, and \(\begin{aligned}i_{t} \tilde{c}_{t}\\\end{aligned}\) is the added new information.

The data pass through the output gate, and the output value o_t of the activation function σ a determines the proportion of the output information, and C_t passes through the function to calculate the output information. The product of tanh(C_t) and o_t is the output information h_t at moment t. The mathematical equations for o_t and h_t are as follows:

o_t = σ (W_o · [h_t−1, x_t] + b_o)       (12)

h_t = o_t tanh(C_t)       (13)

3.2.4. Hybrid Model

The optimal hyperparameters of the SARIMA model were determined by grid search with cross-validation (GridSearchCV), and the predicted values of the SARIMA model were obtained. The EMD algorithm is used as a preprocessing method for signal feature extraction, and each IMF component is used as input information for the LSTM model to obtain the predicted values of the EMD-LSTM model. The SARIMA-EMD-LSTM hybrid model is construed by weight assignment (Figure 2). The mathematical formulation of the hybrid model is:

Y = a . Y₁ + b . Y₂       (14)

Figure 2: Illustration of the SARIMA-EMD-LSTM Hybrid Model Structure

where Y is the forecasting value of the SARIMA-EMD-LSTM hybrid model, Y₁ is the forecasting value of the SARIMA model, Y₂ is the forecasting value of the EMD-LSTM model, and a, b is the weighting factor.

The weight factors a and b are determined based on the MAE values of the EMD-LSTM model and SARIMA model, which take values in the range [0,1], a + b = 1, and the formula is:

\(\begin{aligned}a=\frac{\text { MAE }_{\text {SARIMA }}}{\text { MAE }_{\text {SARIMA }}+\text { MAE }_{\text {EMD-LSTM }}}\\\end{aligned}\)       (15)

\(\begin{aligned}b=\frac{\mathrm{MAE}_{\mathrm{EMD}-\mathrm{LSTM}}}{\mathrm{MAE}_{\mathrm{SARIMA}}+\mathrm{MAE}_{\mathrm{EMD}-\mathrm{LSTM}}}\\\end{aligned}\)       (16)

3.3. Statistical Treatment

This paper uses Python coding, and the IDE uses Python 2021 software. The SARIMA model is constructed by the arima_model package, and the EMD-LSTM model is constructed by the keras package and the PyEMD package.

4. Results and Discussion

4.1. Results of the SARIMA Model

The Seasonal and Trend decomposition using Loess (STL) method is used to decompose the Sino-Russia foreign trade turnover time series to obtain the decomposition values of the seasonal effect, trend effect, and residual volatility effect. The decomposition result shows that there is a clear trend in the time series, with a downward trend from 2014–2016 and an upwards trend from 2016–2020, which requires a differencing operation. Meanwhile, there is a stable cyclicality in trade turnover, with a cycle of 12 months. The Sino-Russia foreign trade volume time series has both seasonal and nonseasonal information, so it needs to be fitted and forecasted using the SARIMA model. After the second-order differences are applied to the time series, the results of the Augmented Dickey-Fuller test (ADF) and Ljung-Box test (LB) show that they are smooth and have no white noise.

Optimal hyper-parameters are determined by using GridSearchCV in machine learning. The search for parameters is performed in two stages, first determining the optimal parameters for the ARIMA model and then determining the seasonal parameters.

Figure 3 shows that the ACF is trailing, with the auto-correlation falling into the 95% confidence interval at p = 2, and the PACF is also trailing, with the auto-correlation falling into the 95% confidence interval at q = 5. The range of values for the parameters of the ARIMA model is set as follows: p ∈ [1,2], d = 2, q ∈ [1, 5]. When the parameters are the following options: (1,2,5), (2,2,5), (2,2,4), and (2,2,1), the value of AIC is less than 530. Only the ARIMA(2,2,1) model has all coefficients with p values significant at the 0.05 level. The range of values for the seasonality parameter in the grid search is set as P ∈ [0, 4], D = 2, Q ∈ [0, 6]. When the ARIMA model is added to the seasonal parameters P = 1, D = 2, Q = 0, the model is significantly optimized, and the AIC value is significantly reduced from 529.45 to 377.68. Therefore, the optimal model obtained by the two-stage search for parameters is SARIMA(2,2,1) × (1,2,0,12).

Figure 3: ACF and PACF Plots

Figure 4 (Normal Q-Q) shows that the residuals are approximately on a straight line, which indicates that the residuals are normally distributed; Figure 4 (Correlogram) shows that all auto-correlations are within the 95% confidence interval. Therefore, the residual series of the SARIMA model is white noise.

Figure 4: The Residual Information of the SARIMA Model

As shown in Table 1, all coefficients of the model are statistically significant, and therefore, the SARIMA(2,2,1) × (1,2,0,12) model can be used to predict the Russian trade turnover time series.

Table 1: Statistical Significance Test of SARIMA Model Coefficients

Note: Standard errors in parentheses, *p < 0.1, **p < 0.05, and ***p < 0.01. Source: Composed by authors.

4.2. The EMD-LSTM Model Results

Each IMF represents the signal functions of each layer after the original signal has been decomposed by EMD. As shown in Figure 6, each IMF satisfies the constraints (Huang et al., 1998).

As shown in Figure 5, IMF1 exhibits strong nonlinear characteristics and contains a large amount of irregular noise. Moreover, it has a high instantaneous frequency and an insignificant periodicity, so it can be regarded as a high-frequency function of the load series. IMF2 demonstrates significant cyclical characteristics, so it can be seen as the cyclical function of the load series. IMF4 has a low instantaneous frequency and an insignificant periodicity, so it can be seen as a low-frequency function of the load series. The residual shows a significant change trend, so it is considered to be the trend function of the load series.

Figure 5: Illustration of EMD Decomposition

Each IFM and residual is used as input information to the LSTM model for prediction. The dataset is normalized in the range [0,1] using MinMaxScaler from the scikit-learn library. The equation is:

\(\begin{aligned}X_{\mathrm{norm}}^{(i)}=\frac{X^{(i)}-X_{\min }}{X_{\max }-X_{\min }}\\\end{aligned}\)       (17)

The dataset is also transformed into a supervised learning dataset, and the value at moment t is predicted by the value at moment t−1. By stacking LSTM hidden layers to increase the depth of the neural network, it can effectively reduce neurons and improve training efficiency and forecast accuracy. In this paper, we use the Keras deep learning framework to quickly build a stacked LSTM neural network model. The sequential model is built by adding an LSTM layer (with 256 units), a dropout layer (Krizhevsky et al., 2017), three LSTM hidden layers (with 128, 64, and 32 units), and finally, a dense layer to aggregate its dimensionality to 1. Moreover, when dropout is 0.5, it generates the most network structures randomly, which effectively enhances the generalization ability of the model and prevents overfitting (Song et al., 2021; Ding, 2022).

The activation function Sigmoid can easily cause problems such as gradient explosion and gradient disappearance, so Relu is used for the activation function. The optimization training is performed using the Adam optimizer, which can adaptively adjust the learning rate (Yi et al., 2020). The time step is set to 1, the feature vector is also 1, and the pitch size is 40. After 300 iterations, the loss rate tends to a stable value close to 0, and convergence is reached. The time series is reconstructed to obtain the predicted values of the EMD-LSTM model. Finally, the predicted values of the SARIMA-EMD-LSTM hybrid model are obtained based on the weighting factors.

As shown in Table 2, we can make the following conclusions:

Table 2: Predicted Performance for Each Model

Source: Composed by authors.

1. The SARIMA model can only capture linear information but not nonlinear information. The forecast result shows that the RMSE is 8.933681, and the MAPE is 14.4404%, which is the worst forecast accuracy of all models, indicating that the long-term forecasting effect of the SARIMA model is not significant.

2. The forecast accuracy of the EMD-LSTM model improved by 2.862% over that of the LSTM model, and this conclusion is consistent with H1.

3. The EMD-LSTM model has the highest forecast accuracy of all models: the MAPE is 7.4304%, which is 2.6398% higher than that of the SARIMA-EMD-LSTM model. This result indicates that the hybrid model does not exploit the linear predictive advantage of the SARIMA model with the nonlinear predictive advantage of the LSTM model. Therefore, this conclusion is inconsistent with H2.

4. The forecast accuracy of the SARIMA-EMD-LSTM hybrid model is lower than that of the EMD-LSTM model, and the forecast accuracy of the SARIMA-LSTM model is lower than that of the LSTM model, which indicates that the hybrid model constructed by the method of weight assignment cannot take advantage of the SARIMA linear and LSTM nonlinear forecasting, so the weight assignment is not the best-combined method for constructing the hybrid model. This conclusion is inconsistent with H3. The hybrid model is generated to effectively compensate for the limitations of the single model, but it does not mean that the single model forecasting is certainly worse than the hybrid model. The reason is that the forecast accuracy depends on the combined method of the hybrid model.

5. The SARIMA model accurately captures the fluctuation pattern of the time series from November 2020 to October 2021. Furthermore, the predicted value for January 2021 is the same as the true value. However, the abnormal decline in trade turnover due to COVID-19 is not reflected in the SARIMA model, which is a shortcoming of the SARIMA model fitting and forecasting.

6. The EMD algorithm effectively resolves the phenomenon of hysteresis generated by the LSTM model. This conclusion is consistent with Hypothesis 1, and it can accurately capture the fluctuation pattern of the Russian trade turnover time series. The ARIMA-LSTM and SARIMA-EMD-LSTM models are significantly influenced by the SARIMA model, as seen in Figure 6(c), with almost the same predicted fluctuation pattern of the time series.

5. Conclusion and Limitations

The Sino-Russia foreign trade turnover time series is complex and has significant seasonal fluctuations, which are transformed into a smooth series by differencing methods. In this paper, the two-stage grid search is used to find the optimal hyper-parameters, and according to the change in AIC values, it is necessary to add seasonal parameters. The EMD algorithm decomposes the time series into a finite number of IMFs and uses them as input information for the LSTM model. Then, the time series were reconstructed to obtain the predicted values of the EMD-LSTM model. The predicted values of the SARIMA-EMD-LSTM hybrid model are obtained by weighting factors. The EMD-LSTM model has the highest forecast accuracy of all models and can be applied to the forecasting of the Sino-Russia foreign trade turnover time series post-epidemic. However, the following limitations remain in this paper:

1. The EMD decomposition algorithm effectively solves the hysteresis phenomenon while improving the forecast accuracy of the time series, as a result in the paper has confirmed. Although during the research, the problem of information leakage did not arise, it is an undisputed fact that when using the EMD algorithm, the problem of information leakage is prone to arise, and no effective solution is proposed in this paper.

2. In this paper, the optimal hyper-parameters are selected after several validations in the construction of the LSTM model, with a strong subjectivity. The Bayesian optimization algorithm is an efficient global optimization algorithm (Torun et al., 2018). It is widely used for hyperparameter optimization in machine learning and deep learning. Using the Bayesian optimization algorithm to determine the optimal hyper-parameters of the LSTM model can improve the forecast accuracy, which is a further research direction.

3. The paper validates that the weight assignment is not the best-combined method for constructing the hybrid model. However, other combined methods, such as the linear-nonlinear hybrid model and correction hybrid model, have not been validated.

4. Foreign trade turnover is affected by various factors, such as politics, geography, war, and epidemics. Only trade turnover is used as the feature variable of the LSTM model in this paper, and other factor are not considered. Russian GDP, international oil prices, and related trade products such as timber, metals, chemicals, machinery, and electricity are the main factors that impact Sino-Russia foreign trade turnover (Zhou & Kang 2021). Therefore, in later research, more macro-political and macroeconomic factors will be taken into account to build a complete system of trade turnover indicators and to improve the forecast accuracy of the model.