Forecasting inter-urban transport demand for a logistics company: A combined grey–periodic extension model with remnant correction

Abstract

Accurately predicting short-term transport demand for an individual logistics company involved in a competitive market is critical to make short-term operation decisions. This article proposes a combined grey–periodic extension model with remnant correction to forecast the short-term inter-urban transport demand of a logistics company involved in a nationwide competitive market, showing changes in trend and seasonal fluctuations with irregular periods different to the macroeconomic cycle. A basic grey–periodic extension model of an additive pattern, namely, the main combination model, is first constructed to fit the changing trends and the featured seasonal fluctuation periods. In order to improve prediction accuracy and model adaptability, the grey model is repeatedly modelled to fit the remnant tail time series of the main combination model until prediction accuracy is satisfied. The modelling approach is applied to a logistics company engaged in a nationwide less-than-truckload road transportation business in China. The results demonstrate that the proposed modelling approach produces good forecasting results and goodness of fit, also showing good model adaptability to the analysed object in a changing macro environment. This fact makes this modelling approach an option to analyse the short-term transportation demand of an individual logistics company.

Introduction

The main business of transport-oriented logistics companies is freight transport. Accurately predicting transport demand is very important for these logistics companies to arrange transportation tasks, optimise resource allocation and reduce the operating cost of production. Especially for large-scale transportation logistics companies engaged in nationwide business, accurately predicting short-term inter-urban transport demand can help to optimise transport plans and vehicle scheduling and improve profitability. ¹ At present, the logistics industry in China is in an upgrading period of integration; a few large logistics enterprises have been growing. However, the efficiency of resource integration and utilisation needs to be improved. This article mainly focuses on short-term inter-urban transport demand forecasting for large transport-oriented logistics companies, providing support to short-term operation decision-making.

West et al. ² acknowledge that the uncertainty, randomness and fuzziness of freight transport demand leads to complexity when forecasting. Freight transport demand forecasting methods can be divided into qualitative and quantitative forecasting, according to Zhou and Dai. ³ Qualitative forecasting mainly uses the personal experience and subjective judgement ability of professionals or decision-makers to find out the law and to make judgements about the future. This forecasting method is simple, but subjective and one sided. Its prediction accuracy is also not high. The qualitative forecasting methods used in freight transport demand forecasting mainly include personal judgement, the Delphi technique and the subjective probability method, which are usually used as supplementary decision-making tools. Quantitative forecasting methods use a variety of mathematical models to predict future development based on historical statistical data and relevant information. Commonly used methods for quantitatively predicting transport demand include time series analysis, regression analysis, the transport coefficient method, the velocity ratio method and the neural network prediction method.^4–9

The quantitative prediction process includes model selection, model calculation and result adjustment. ¹⁰ The first step is to select a model or a model combination according to the characteristics of the object to be forecast and the decision-making purpose. The second step is to carry out data collection and model calculation. The third step is to analyse the reliability and validity of the results and to carry out model adjustment if needed. In the literature, forecasting macro national transportation demand or regional long-term comprehensive transportation demand has been mostly studied.^11–13 Single models are mostly used, and a few combined models are seen. Most of the prediction results show that the factors influencing macro and long-term freight transport demand are relatively clear, the trend and periods of which are closely related to the trend and periods of national and regional economic development. This kind of prediction is not sufficient to support the microscopic operation management of a firm or company, however, because the short-term freight transportation demand of a firm or company suffers from more uncertain influencing factors. The freight transport demand of a firm or company is related not only to the influencing factors of macro transportation demand but also to the micro influencing factors such as market position, marketing strategy and customer relationships. McCarthy ¹⁴ and Lee et al. ¹⁵ acknowledge that these two kinds of factors are continually changing, and the connections between them are difficult to accurately express, showing a kind of grey system relation. Fang et al. ¹⁶ argued that to some extent, the micro transportation demand of a firm or company will rise and fall along with macroeconomic fluctuations, but its fluctuation period and development trend are often unclear and fuzzy.

The grey model GM(1,1) and its improved models have been successfully applied to forecast long-term and short-term transportation demands and have shown promising results. In fact, grey models do not require data as strictly as other regression models and have good adaptability to short-term and data-limited predictions. In terms of processing the seasonal fluctuations and trends of a transportation demand time series with a fixed period, it has been suggested to separate the time series by smoothing the data into a trend component, fluctuation component and error component and then using different models to fit the period and the trend, respectively, and combine them, according to Barrow. ¹⁷ This method was successfully applied to forecast the transport demand of a small logistics company (an export container transportation business near the port area) developed by Fang et al. ¹⁶ Seasonal Autoregressive Integrated Moving Average (SARIMA) model is another forecasting method usually used to analyse non-stationary time series with trend and seasonality. ¹⁸ However, when the seasonal fluctuation period is not very stable or clear, this method may confuse the subsequent analysis. At this time, the data series needs to be further modelled to find the trend and periods. A Fourier series indicates that any time series function can be seen as the superposition of infinite vibration waves with different frequencies. ¹⁹ Brockwell and Davis ²⁰ stated that power spectrum analysis and periodic extension are the two most effective methods to identify and extract the dominant periods of a time series.

For a time series with trend and periodic waves, it can neither reveal the periodic feature very well by only using grey system models nor can it reveal the trend by only using periodic extension models (PEMs). The approach proposed by Fang et al. ¹⁶ can cope with the effect of macroeconomic fluctuation on small-sized enterprises. However, data pre-processing eliminates the irregular data, which helps to improve the stability of the observed data series. In this article, the sample logistics company’s business data have the dual characteristics of trend and non-fixed periods, which result from the additive fluctuations of macro and micro factors. Liu and Lin ²¹ indicated that trying a combination model of the grey model and PEM gives a way to solve this problem.

The objective of this study is to construct a combined grey-PEM with remnant correction to cope with short-term and data-limited prediction problem with trend and non-stationary periods, which is characterised by the sample business data from a real logistics company. The context of the model is characterised by one distinct element: the short time series with trend and non-stationary periods. In the modelling process, the main combination model is constructed by superposing a grey model GM(1,1) with a PEM, in which GM(1,1) is chosen to fit the trend, and the PEM is constructed to capture the featured seasonal fluctuation periods. In order to improve prediction accuracy and model adaptability, remnant correction is carried out by fitting the remnant tail time series of the main combination model until the prediction accuracy is satisfied.

The rest of this article is organised as follows: In the following section, a detailed modelling process and the model architecture are presented. This is followed by model application, including data description, modelling results, forecasting performance and model comparison. Conclusions and implications are then drawn.

Modelling

The basic idea of building a GM(1,1)-periodic extension model with remnant correction (RC-GMPEM) is to make full use of the trend characteristics of exponential growth in the GM(1,1) model to improve the predicting accuracy. ²¹ While taking into account the seasonal fluctuations in the time series, the PEM is used to identify precisely the seasonal fluctuation period of the observed time series and carry out the forecasting revision in the observed fluctuant points; then, the trend component and seasonal component are superposed using the additive modelling of the seasonal time series in order to obtain the main combination model. ¹⁷

The modelling process of RC-GMPEM is illustrated in Figure 1. The modelling process can be divided into five steps. The first step is to test the smoothness and quasi-exponential of the original time series. Passing the test of the quasi-exponential is the premise to establish the GM(1,1) model. The second step is to build the GM(1,1) model of the time series according to the modelling method of grey system theory, if the time series has passed the test and fits the trend components of the time series. In the third step, the PEM is built to identify the fluctuation period of the time series; this needs to correct the peak points of the large fitting errors caused by seasonal fluctuations. The fourth step is to combine the GM(1,1) model with the PEM to obtain the main combination model (GMPEM) according to the additive modelling idea of the seasonal time series and then test the prediction accuracy of the main combination model. In the fifth step, GM(1,1) is repeatedly modelled to fit the remnant tail series of the main combination model until the prediction accuracy is satisfied in the fourth step and then we have the final RC-GMPEM. We describe the detailed modelling steps in the following subsections.

Figure 1.

Modelling process of the RC-GMPEM.

Test the smoothness and quasi-exponential of the time series

Assuming that $X^{(0)}=(x^{(0)}(1),x^{(0)}(2),\dots ,x^{(0)}(n))$ is the inter-urban transport demand time series of a logistics company. Obviously, $x^{(0)}(k)\ge 0,k=1,2,\dots ,n$ . Its 1-accumulated generating operation (AGO) $X^{(1)}$ is $X^{(0)}$ , and $X^{(1)}=(x^{(1)}(1),x^{(1)}(2),\dots ,x^{(1)}(n))$ , where $x^{(1)}(k)=\sum _{i=1}^k{x}^{(0)}(i)$ , $k=1,2,\dots ,n$ . Similarly, $X^{(r)}=(x^{(r)}(1),x^{(r)}(2),\dots ,x^{(r)}(n))$ is called r-AGO series of $X^{(0)}$ , where $x^{(r)}(k)=\sum _{i=1}^k{x}^{(r-1)}(i),k=1,2,\dots ,n$ .

Test the smoothness of ${\mathrm{X}}^{(0)}$

The smoothness of the series can be measured by the smoothness ratio $\rho (k)$ . Smaller $\rho (k)$ means better smoothness of the time series

$$\rho \left(k\right)=\frac{x^{(0)}\left(k\right)}{\sum _{i=1}^{k-1}{x}^{(0)}\left(0\right)}=\frac{x^{(0)}\left(k\right)}{x^{(1)}(k-1)},k=2,3,\dots ,n$$

(1)

When the $\rho (k)$ of the time series $X^{(0)}$ meets the following three conditions, we say the time series $X^{(0)}$ meets the quasi-smoothness conditions for grey forecasting. ²¹

1. $\rho (k+1)/\rho (k)<1,k=2,3,\dots ,n-1$

2. $\rho (k)\in [0,\epsilon ],k=3,4,\dots ,n$

3. $\epsilon <0.5$

where $\epsilon$ is a term of shocking disturbance, representing a disturbed sequence. ²² We have to seek an appropriate model or combined models by eliminating shock wave disturbances in order to recover the true state of the systems’ behavioural data, so that the accuracy of the consequent predictions can be improved. Therefore, if the original time series does not pass the quasi-smoothness test, we have to adopt some appropriate treatment to the original data series as follows until the time series passes the test according to Liu and Lin. ²¹

Taking a proper constant c, and making $y^{(0)}(k)=x^{(0)}(k)+c,k=1,2,\dots ,n$ , until the smoothness ratio of the new time series $Y^{(0)}=(y^{(0)}(1),y^{(0)}(2),\dots ,y^{(0)}(n))$ passes the quasi-smoothness test.

Test the quasi-exponential of ${\mathrm{X}}^{(r)}$

Generally, the randomness of the original time series may be reduced after repeated AGO, and the generated series will show the law of exponential growth according to Benítez et al. ²³ At the same time, Kayacan et al. ²⁴ stated that, if the r-AGO series of $X^{(0)}$ has had obvious exponential growth, additional AGO may destroy its regularity and make the exponential law grey. Hence, AGO should be stopped in the right scope.

For $X^{(r)}=(x^{(r)}(1),x^{(r)}(2),\dots ,x^{(r)}(n))$ , if ${\rho }^{(r)}(k)\in [a,b],b-a=\delta \mathrm{and}\epsilon <0.5$ for any k, we can say that $X^{(r)}$ is quasi-exponential. ${\rho }^{(r)}(k)$ is the class ratio of $X^{(r)}$ , ${\rho }^{(r)}(k)=x^{(r)}(k)/x^{(r)}(k-1);k=2,3,\dots ,n$ .

Because ${\rho }^{(1)}(k)=x^{(1)}(k)/x^{(1)}(k-1)=x^{(0)}(k)+x^{(1)}(k-1)/x^{(1)}(k-1)=1+\rho (k)$ , and the original time series $X^{(0)}$ or $Y^{(0)}$ that has passed the smoothness test meets $\rho (k)<0.5$ , and ${\rho }^{(1)}(k)\in [1,1.5),\epsilon <0.5$ . Therefore, a time series that has passed the smoothness test just needs one-time AGO, whose 1-AGO series $X^{(1)}$ could be quasi-exponential.

Build the GM(1,1) model for the generated transport demand series

By building the differential equation of the GM(1,1) model of the generated series $X^{(1)}$ , we can obtain the time response function by estimating the parameters.

The GM(1,1) differential equation of $X^{(1)}$ is

$$\frac{d{x}^{(1)}}{\mathrm{dt}}+a{x}^{(1)}=b$$

(2)

where a and b are the waiting-estimated parameters.

At the moment k, when the time interval $\Delta t$ is 1 unit, the discrete form of $d{x}^{(1)}/\mathrm{dt}$ can be approximately expressed as

$$\frac{\Delta x^{(1)}}{\Delta t}=x^{(1)}(k+1)-x^{(1)}\left(k\right)=x^{(0)}(k+1)$$

(3)

When $\Delta t$ is 1 unit, the discrete value of the variable $x^{(r)}(t)$ can be approximately expressed as

$$x^{(1)}\left(t\right)\approx \frac{1}{2}[x^{(1)}\left(k\right)+x^{(1)}(k+1)]$$

(4)

Assuming that $Z^{(1)}=(z^{(1)}(2),z^{(1)}(3),\dots ,z^{(1)}(n))$ is the nearest-neighbour mean generating series of $X^{(1)}$ , where $z^{(1)}(k)=0.5(X^{(1)}(k)+X^{(1)}(k-1)),k=2,3,\dots ,n$ . Therefore, the discrete form of the first-order one-variable differential equation of GM(1,1), $X^{(1)}$ , can be expressed as

$$x^{(0)}(k+1)=-a{z}^{(1)}(k+1)+b$$

(5)

Let $k=1,2,\dots ,n-1$ , and substituting into equation (5), we obtain the equations

$$\{\begin{array}{c}\hfill x^{(0)}\left(2\right)=-a{z}^{(1)}\left(2\right)+b\hfill \\ \multicolumn{1}{c}{x^{(0)}\left(3\right)=-a{z}^{(1)}\left(3\right)+b}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{x^{(0)}\left(n\right)=-a{z}^{(1)}\left(n\right)+b}\end{array}$$

(6)

where a and b are the waiting-estimated parameters.

Assuming that

$$\hat{a}={[a,b]}^T,Y=\left[\begin{array}{c}\hfill x^{(0)}\left(2\right)\hfill \\ \multicolumn{1}{c}{x^{(0)}\left(3\right)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{x^{(0)}\left(n\right)}\end{array}\right],B=\left[\begin{array}{c}\hfill -z^{(1)}\left(2\right)1\hfill \\ \multicolumn{1}{c}{-z^{(1)}\left(3\right)1}\\ \multicolumn{1}{c}{⋮⋮}\\ \multicolumn{1}{c}{-z^{(1)}\left(n\right)1}\end{array}\right]$$

then equation (6) is equivalent to

$$Y=B\hat{a}$$

(7)

Then, we can obtain the expression of the waiting-estimated parameters a and b

$$\widehat{a}=\left[\begin{array}{l}a\\ b\end{array}\right]={\left(B^{T}B\right)}^{-1}{B}^{T}Y=\left[\begin{array}{l}\frac{\frac{1}{n-1}\left({\displaystyle \sum _{k=2}^n{x}^{(0)}\left(k\right)}\right)\cdot \left({\displaystyle \sum _{k=2}^n{z}^{(1)}\left(k\right)}\right)-{\displaystyle \sum _{k=2}^n{x}^{(0)}\left(k\right)\cdot z^{(1)}\left(k\right)}}{{\displaystyle \sum _{k=2}^n{\left(z^{(1)}\left(k\right)\right)}^2}-\frac{1}{n-1}{\left({\displaystyle \sum _{k=2}^n{z}^{(1)}\left(k\right)}\right)}^2}\\ \frac{1}{n-1}\left[{\displaystyle \sum _{k=2}^n{x}^{(0)}\left(k\right)}+a{\displaystyle \sum _{k=2}^n{z}^{(1)}\left(k\right)}\right]\end{array}\right]$$

According to equation (4), by solving the differential equation, we obtain

$$x^{(1)}\left(t\right)=c{e}^{-\mathrm{at}}+\frac{b}{a}$$

(8)

The discrete form is

$$x^{(1)}\left(k\right)=c{e}^{-\mathrm{ak}}+\frac{b}{a}$$

(9)

Let $k=1$ . We can obtain the parameter $c=(x^{(1)}(1)-(b/a))e^a=(x^{(0)}(1)-(b/a))e^a$ . Therefore, $x^{(1)}(k)$ can be expressed as

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

(10)

Let $k=k+1$ , and substituting into equation (10), we obtain

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

(11)

Equation (11) is the response function of the 1-AGO series of the GM(1,1) of the inter-urban transport demand time series.

Thus, according to $x^{(1)}(k+1)=\sum _{i=1}^{k+1}{x}^{(0)}(i)=\sum _{i=1}^k{x}^{(0)}(i)+x^{(0)}(k+1)$ , we can obtain the following GM(1,1) of the transport demand time series by applying the inverse AGO approach

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

(12)

Build the PEM model for the generated transport demand series

In order to overcome the shortcoming that the GM(1,1) model is not excellent at fitting seasonal fluctuations, we can build the PEM of the remnant series, which is constituted by the predicted values and observed values of the GM(1,1) model. Through this, we can extract the dominant period by discerning the periodic fluctuation phenomena of the remnant series and then modify the predicted values on the inflection point of periodic variation. The modelling steps are as follows:

1. Build the remnant series $X\prime$ of the transport demand model GM(1,1). Assuming that $x\prime (k)=x^{(0)}(k)-{\hat{x}}^{(0)}(k)$ , the remnant series is $X\prime =(x\prime (1),x\prime (2),\dots ,x\prime (n))$ .

2. Build the mean generating function of the remnant series $X\prime$ as follows

   $${\overline{x}}_m\left(i\right)=\frac{\sum _{j=0}^{n_m-1}x\prime (i+\mathrm{jm})}{n_m},i=1,2,\dots ,m,1\le m\le M$$
   (13)

In equation (13), n is the length of the remnant series $X\prime$ , $n_m$ is the largest integer smaller than $n/m$ and M is the largest integer smaller than $n/2$ .

Through these calculations, we can obtain the mean generating function matrix of $X\prime$

$$\left[\begin{array}{c}\hfill {\overline{x}}_1\left(1\right){\overline{x}}_2\left(1\right){\overline{x}}_3\left(1\right)\cdots {\overline{x}}_M\left(1\right)\hfill \\ \multicolumn{1}{c}{{\overline{x}}_2\left(2\right){\overline{x}}_3\left(2\right)\cdots {\overline{x}}_M\left(2\right)}\\ \multicolumn{1}{c}{{\overline{x}}_3\left(3\right)\cdots {\overline{x}}_M\left(3\right)}\\ \multicolumn{1}{c}{\ddots ⋮}\\ \multicolumn{1}{c}{{\overline{x}}_M\left(M\right)}\end{array}\right]$$

(14)

We now make the periodic extension the mean generating function ${\overline{x}}_m(i)$

$$f_m\left(k\right)={\overline{x}}_m\left(i\right),k=i[mod\left(m\right)],k=1,2,\dots ,n$$

(15)

In equation (15), $f_m(k)$ is called the extensional function of the mean generating function.

• 3. Extract the dominant period m. According to variance analysis theory, we can use equation (16) to test whether the series $X\prime$ has a dominant hidden period m

  $$F^{(m)}=\frac{(n-m)S^{(m)}}{(m-1)S}$$
  (16)

where $F^{(m)}$ is an F distribution of degrees of freedom $(m-1,n-m)$ , $S^{(m)}=\sum _{i=1}^m{n_i({\overline{x}}_m(i)-\overline{x})}^2$ , $n_i=n/i$ , $\overline{x}=(\sum _{i=1}^{n}x\prime (i))/n$ and $S=\sum _{i=1}^m{\sum _{j=1}^n[x\prime (i+(j-1)m)-{\overline{x}}_m(i)]}^2$ .

For a given confidence level $\alpha$ , if $F^{(m)}>F_{\alpha }(m-1,n-m)$ , we say that there exists a dominant period m in the series $X\prime$ .

• 4. Extract the other dominant periods. We generate a new series $X″$ by subtracting the element of the remnant series $X\prime$ from the corresponding extensional function, that is

  $$x″\left(k\right)=x\prime \left(k\right)-f_m\left(k\right)$$
  (17)

By repeating steps (2) and (3) for the new series $X″$ , we can extract the other dominant periods.

• 5. Superposition. Superpose the values of different periods at the same time and write it as $f(k)$

  $$f\left(k\right)=\sum _{i=1}^m{f}_i\left(k\right)$$
  (18)

$f(k)$ is the PEM built by superposing the period.

Build the main combination model (GMPEM)

According to the additive modelling idea of seasonal time series, we make $x\prime (k)$ approximately equal to $f(k)$ , and we can then obtain the combination expression of the GMPEM

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

(19)

${\widehat{\widehat{x}}}^{(0)}(k)$ is the predictive value of the transport demand calculated by the GMPEM.

Test the prediction accuracy of the GMPEM

The prediction accuracy of the GMPEM can be measured by variance ratio c and small error probability p. ²¹ The calculations are as follows

$$c=\frac{S_2}{S_1}$$

(20)

$$p=P\{|\tilde{x}\prime \left(k\right)-\tilde{x}\prime |<0.6745S_1\}=\frac{n\prime }{n}$$

(21)

where $x\prime (k)=x^{(0)}(k)-{\widehat{\widehat{x}}}^{(0)}(k)$ is the remnant of the combination model, $\tilde{x}\prime =(1/n)\sum _{k=1}^n\tilde{x}\prime (k)$ is the remnant mean, $S_1=\sqrt{(1/(n-1)){\sum _{k=1}^n[x^{(0)}(k)-(1/n)\sum _{k=1}^n{x}^{(0)}(k)]}^2}$ is the standard deviation of original freight time series, $S_2=\sqrt{(1/(n-1)){\sum _{k=1}^n[\tilde{x}\prime (k)-\tilde{x}\prime ]}^2}$ is the standard deviation of the remnants of the combination model, $n\prime$ is the number of basic events that meets $|\tilde{x}\prime (k)-\tilde{x}\prime |<0.6745S_1$ and n is the length of the original time series.

The prediction accuracy can be quantified into different levels according to the values of c and p (see Table 1). If the values of c and p are not at the same level, we take the greater value to be the prediction accuracy level of the model.

Table 1.

Prediction accuracy level measured by c and p.

c	<0.35	<0.50	<0.65	≥0.65
p	>0.95	>0.80	>0.70	≤0.70
Prediction	A	B	C	D
Accuracy level	Good	Qualified	Weak	Unqualified

Source: Liu and Lin. ²¹

Correct the random error of the GMPEM by remnant correction

In order to reduce the influence of random error of the GMPEM on the prediction accuracy, we build the GM(1,1) model of the remnants of the combination model to correct the random error, so that we can improve the precision of the prediction accuracy further.

Assuming that ${\widehat{\widehat{x}}}^{(0)}(k)$ is the predicted series of the combination model (GMPEM), and ${\widehat{\widehat{X}}}^{(0)}(k)=({\widehat{\widehat{x}}}^{(0)}(1),{\widehat{\widehat{x}}}^{(0)}(2),\dots ,{\widehat{\widehat{x}}}^{(0)}(n))$ , the remnant of the GMPEM can be defined as

$$\tilde{x}\prime \left(k\right)=x^{(0)}\left(k\right)-{\widehat{\widehat{x}}}^{(0)}\left(k\right)$$

(22)

Therefore, we can obtain the remnant series of the GMPEM, namely, ${\tilde{X}}^{\prime }=(\tilde{x}\prime (1),\tilde{x}\prime (2),\dots ,\tilde{x}\prime (n))$ . If there is a $k_0$ satisfying (1) $\forall k>k_0$ , $\tilde{x}\prime (k)$ with consistent symbols, and (2) $n-k_0\ge 4$ , we say the remnant tail time series $(\tilde{x}\prime (k_0),\tilde{x}\prime (k_0+1),\dots ,\tilde{x}\prime (n))$ can be modelled. Hence, we mark it as

$${\widetilde{X\prime }}^{(0)}=(\tilde{x}\prime \left(k_0\right),\tilde{x}\prime (k_0+1),\dots ,\tilde{x}\prime \left(n\right))$$

Based on this method of modelling GM(1,1), we can build a GM(1,1) model of the 1-AGO series ${\widetilde{X\prime }}^{(1)}=(\tilde{x}{\prime }^{(1)}(k_0),\tilde{x}{\prime }^{(1)}(k_0+1),\dots ,\tilde{x}{\prime }^{(1)}(n))$ of ${\widetilde{X\prime }}^{(0)}$

$${\widehat{{\tilde{x}}^{\prime }}}^{(1)}\left(k+1\right)=\left({\tilde{x}}^{\prime }\left(k_0\right)-\frac{b^{\prime }}{a^{\prime }}\right)e^{-a^{\prime }(k-k_0+1)}+\frac{b^{\prime }}{a^{\prime }},k\ge k_0-1$$

(23)

Therefore, the GM(1,1) form of the fitted series ${\widehat{\tilde{X\prime }}}^{(0)}=(\widehat{\tilde{x}\prime }(k_0),\widehat{\tilde{x}\prime }(k_0+1),\dots ,\widehat{\tilde{x}\prime }(n))$ of ${\widetilde{X\prime }}^{(0)}$ is

$$\widehat{{\tilde{x}}^{\prime }}\left(k+1\right)=\left(1-e^{a^{\prime }}\right)\left({\tilde{x}}^{\prime }\left(k_0\right)-\frac{b^{\prime }}{a^{\prime }}\right)e^{-a^{\prime }(k-k_0+1)},k\ge k_0-1$$

(24)

Based on the additive modelling idea of seasonal time series, by combining equations (19) and (24), we can obtain the final RC-GMPEM model, which has been corrected using remnants

$$\begin{array}{cc}\multicolumn{1}{c}{{\widehat{\widehat{\widehat{x}}}}^{(0)}\left(k\right)}& ={\widehat{\widehat{x}}}^{(0)}\left(k\right)+\widehat{\tilde{x}\prime }\left(k\right)={\widehat{x}}^{(0)}\left(k\right)+f\left(k\right)+\widehat{\tilde{x}\prime }\left(k\right)\hfill \\ \multicolumn{1}{c}{}& =(1-e^a)(x^{(0)}\left(1\right)-\frac{b}{a})e^{-a(k-1)}+\delta \left(k\right)(1-e^{a\prime })\hfill \\ \multicolumn{1}{c}{}& (\tilde{x}\prime \left(k_0\right)-\frac{b\prime }{a\prime })e^{-a\prime (k-k_0)}+f\left(k\right)\hfill \end{array}$$

(25)

where $\delta (k)=\{\begin{array}{c}\hfill 1,k\ge k_0\hfill \\ \multicolumn{1}{c}{0,k<k_0}\end{array}$ .

Model applications

Data

The proposed model is applied to predict the short-term transport demand in between three hub cities for a large five-star transportation logistics company in China, whose major business is domestic road less-than-truckload (LTL) freight transportation between big cities. The data cover all the transport volumes of Shanghai–Guangzhou, Shanghai–Shenzhen and Guangzhou–Shenzhen during a 36-month period (January 2009 to December 2011). The real raw data are shown in Table 5 in Appendix 1. These three time series show a tendency towards irregular periods. Our proposed model is used to analyse the time series.

Modelling results

Let $X_1^{(0)}=(x_1^{(0)}(1),x_1^{(0)}(2),\dots ,x_1^{(0)}(36))$ , $X_2^{(0)}=(x_2^{(0)}(1),x_2^{(0)}(2),\dots ,x_2^{(0)}(36))$ and $X_3^{(0)}=(x_3^{(0)}(1),x_3^{(0)}(2),\dots ,x_3^{(0)}(36))$ be the original time series of the transport volume of Shanghai–Guangzhou, Shanghai–Shenzhen and Guangzhou–Shenzhen, respectively, and ${\hat{x}}_1^{(0)}(k)$ , ${\hat{x}}_2^{(0)}(k)$ , ${\hat{x}}_3^{(0)}(k)$ be the form of GM(1,1). We obtain the three GM(1,1) models by carrying out the calculations in section ‘Build the GM(1,1) model for the generated transport demand series’ in MATLAB

$${\hat{x}}_1^{(0)}\left(k\right)=1257521.78e^{0.0288759(k-1)}$$

$${\hat{x}}_2^{(0)}\left(k\right)=1550808.88e^{0.0100999(k-1)}$$

$${\hat{x}}_3^{(0)}\left(k\right)=1637979.74e^{0.0100999(k-1)}$$

Based on a confidence level $\alpha =0.05$ , we obtain the three PEMs and their dominant periods through the method presented in section ‘Build the PEM model for the generated transport demand series’

$$f^1\left(k\right)=\sum _m{f}_m^1\left(k\right),m=6,8,10,11,14,12$$

$$f^2\left(k\right)=\sum _m{f}_m^2\left(k\right),m=4,6,10,8,11,13$$

$$f^3\left(k\right)=\sum _m{f}_m^3\left(k\right),m=4,6,11,13$$

The main combination models (GMPEMs) are

$${\widehat{\widehat{x}}}_1^{(0)}\left(k\right)=1257521.78e^{0.0288759(k-1)}+f^1\left(k\right)$$

$${\widehat{\widehat{x}}}_2^{(0)}\left(k\right)=1550808.88e^{0.0100999(k-1)}+f^2\left(k\right)$$

$${\widehat{\widehat{x}}}_3^{(0)}\left(k\right)=1637979.74e^{0.0100999(k-1)}+f^3\left(k\right)$$

The prediction accuracy of the three GMPEMs is shown in Table 2.

Table 2.

Prediction accuracy of the three GMPEMs.

Sample	c	p	Accuracy level
Shanghai–Guangzhou	0.21	1	Good
Shanghai–Shenzhen	0.23	1	Good
Guangzhou–Shenzhen	0.52	0.72	Weak

GMPEM: GM(1,1)-periodic extension model.

Table 2 shows that the prediction accuracy of the GMPEMs fits the observation series of Shanghai–Guangzhou and Shanghai–Shenzhen well but is weak for Guangzhou–Shenzhen. There is a possibility that the prediction accuracy differences among the different samples could be due to random errors in the time series, and therefore, if feasible, it is necessary to eliminate the random error disturbance adversely affecting the prediction accuracy as much as possible, especially for the time series of Guangzhou–Shenzhen.

We use the method in section ‘Correct the random error of the GMPEM by remnant correction’ to correct the random errors for the three sample time series. The results are as follows:

1. For Shanghai–Guangzhou, there is the presence of $k_0=27$ , making the remnant tail time series $\{\widetilde{{x\prime }_3}(27),\widetilde{{x\prime }_3}(28),\dots ,\widetilde{{x\prime }_3}(36)\}$ able to be modelled as a GM(1,1) with $\widetilde{{x\prime }_1}(k)=x_1^{(0)}(k)-{\tilde{x}}_1^{\prime }(k)=x_1^{(0)}(k)-{\overline{x}}_1^{(0)}(k)$ . The estimated parameters are $\left[\begin{array}{c}\hfill a\prime \hfill \\ \multicolumn{1}{c}{b\prime }\end{array}\right]=\left[\begin{array}{c}\hfill 0.138257\hfill \\ \multicolumn{1}{c}{-261528}\end{array}\right]$ , and the GM(1,1) of the remnant tail time series are

   $$\widehat{\tilde{{x\prime }_1}}\left(k\right)=-273586.36e^{-0.138257(k-27)}k\ge 27$$

2. For Shanghai–Shenzhen, we detect the presence of $k_0=36$ , indicating that it does not meet the requirements of building a GM(1,1) for the remnant tail time series. Therefore, the GMPEM cannot be corrected.

3. For Guangzhou–Shenzhen, there is the presence of $k_0=28$ ; therefore, the remnant tail time series $\{\widetilde{{x\prime }_3}(28),\widetilde{{x\prime }_3}(29),\dots ,\widetilde{{x\prime }_3}(36)\}$ can be modelled as a GM(1,1) with $\widetilde{{x\prime }_3}(k)=x_3^{(0)}(k)-{\overline{x}}_3^{(0)}(k)$ . The estimated parameters $\left[\begin{array}{c}\hfill a\prime \hfill \\ \multicolumn{1}{c}{b\prime }\end{array}\right]=\left[\begin{array}{c}\hfill 0.052402\hfill \\ \multicolumn{1}{c}{-403970}\end{array}\right]$ , and the GM(1,1) of the remnant tail time series are

   $$\widehat{\tilde{{x\prime }_3}}\left(k\right)=-392552.66e^{-0.052402(k-28)}k\ge 28$$

After random error correction for the remnant tail time series of the main combination forecasting model, we can obtain the final combination forecasting models (RC-GMPEMs)

$${\widehat{\widehat{\widehat{x}}}}_1^{(0)}\left(k\right)=1257521.78e^{0.0288759(k-1)}-\delta \left(k\right)273586.36e^{-0.138257(k-27)}+f^1\left(k\right)$$

(26)

where if $k<27$ , then $\delta (k)=0$ ; if $k\ge 27$ , then $\delta (k)=1$

$${\widehat{\widehat{\widehat{x}}}}_2^{(0)}\left(k\right)={\widehat{\widehat{x}}}_2^{(0)}\left(k\right)=1550808.88e^{0.0100999(k-1)}+f^2\left(k\right)$$

(27)

$${\widehat{\widehat{\widehat{x}}}}_3^{(0)}\left(k\right)=1637979.74e^{0.0100999(k-1)}-\delta \left(k\right)392552.66e^{-0.052402(k-28)}+f^3\left(k\right)$$

(28)

where if $k<27$ , then $\delta (k)=0$ ; if $k\ge 27$ , then $\delta (k)=1$ .

Forecasting results and model comparison

As a comparison analysis of forecasting accuracy in the sample, we use the additive polynomial trend model ¹⁶ in equation (29), SARIMA model and RC-GMPEMs to forecast the same time series

$$\hat{y}=(\mu +{\beta }_{1}x+{\beta }_2{x}^2+{\beta }_3{x}^3+s_i)(1+{\overline{S}}_c)$$

(29)

where $\hat{y}$ denotes the dependent variable; x is the predictor variable; $\mu ,{\beta }_1,{\beta }_2$ and ${\beta }_3$ are constant parameters; $s_i$ is the random error; and ${\overline{S}}_c$ denotes the periodicity factor.

The forecast results are shown in Figures 2–4. The dominant period of SARIMA models is 13; the predicted values are obtained from forecasting point 14.

Figure 2.

Fitted transport demand and actual demand of Shanghai–Guangzhou.

Figure 3.

Fitted transport demand and actual demand of Shanghai–Shenzhen.

Figure 4.

Fitted transport demand and actual demand of Guangzhou–Shenzhen.

Figures 2–4 show that the RC-GMPEMs have better prediction performance than the other two. Specific criteria, stated by Nelder and Wedderburn, ²⁵ such as the coefficient of determination $R^2$ and F-test, are used to measure the goodness of fit of a proposed forecasting model. However, this is not the case with nonlinear models. According to Gujarati, ²⁶ the ${\overline{R}}^2$ can be used to measure the goodness of fit of the proposed RC-GMPEMs

$${\overline{R}}^2=1-\frac{\sum _{i=1}^n{{\hat{u}}_i}^2}{\sum _{i=1}^n{(x_i-\overline{x})}^2}$$

(30)

where n is the observation number and ${\hat{u}}_i=x_i-{\hat{x}}_i$ , where ${\hat{x}}_i$ are the estimated values from the RC-GMPEMs. The results are shown in Table 3.

Table 3.

Goodness of fit.

Sample	Goodness of fit ${\overline{R}}^2$
	RC-GMPEMs	SARIMA	Polynomial
Shanghai–Guangzhou	0.9603	0.4451	0.8882
Shanghai–Shenzhen	0.9392	0.2612	0.8394
Guangzhou–Shenzhen	0.8296	0.3246	0.7190

RC-GMPEM: GM(1,1)-periodic extension model with remnant correction; SARIMA: Seasonal Autoregressive Integrated Moving Average.

In terms of ${\overline{R}}^2$ in Table 3, RC-GMPEMs have the best performance, and SARIMA models have the worst performance. The results denote that the polynomial trend models and the RC-GMPEMs have advantages over the SARIMA model in coping with short time series with trend and seasonality.

For further comparative analysis, the mean absolute relative error is used as another dimensionless index to measure the forecasting accuracy of the polynomial trend model and the proposed RC-GMPEMs. We computed the relative errors of the polynomial trend models and the RC-GMPEMs. Figures 5–7 show that the relative errors of the RC-GMPEMs are obviously smaller than those of the polynomial trend models.

Figure 5.

Relative errors of the forecasted demand of Shanghai–Guangzhou.

Figure 6.

Relative errors of the forecasted demand of Shanghai–Shenzhen.

Figure 7.

Relative errors of the forecasted demand of Guangzhou–Shenzhen.

The mean absolute relative errors listed in Table 4 show that RC-GMPEMs perform better than the polynomial trend models, which demonstrate that the RC-GMPEMs have advantages over the polynomial trend models in coping with short time series with trend and non-stationary periods by carrying out remnant corrections on the main combination model.

Table 4.

Mean absolute relative errors for the polynomial trend models and RC-GMPEMs.

Forecasting method	Shanghai–Guangzhou	Shanghai–Shenzhen	Guangzhou–Shenzhen
Polynomial trend model	11.92	10.43	12.23
RC-GMPEM	7.53	6.68	9.94

RC-GMPEM: GM(1,1)-periodic extension model with remnant correction.

Conclusion

Classical transportation economics usually focuses on macro comprehensive transportation demand in a state, a region or in between two cities, which offers a research approach to revealing the relations of transportation demand. These relations are helpful for logistics companies but are insufficient. For an individual logistics company, the market denotes the macroeconomic environment and the transportation industry environment, as well as its own market operation strategies (such as service level, pricing strategy, etc.). Therefore, forecasting the transportation demand of an individual company involved in the market will support the logistics company in making its market operation decision-making. In this article, the grey model GM(1,1) combined with a PEM is used to forecast a company’s short-term transportation demand for two reasons. First, the grey model is used to capture the trend component because the factors affecting a company’s transportation demand are more complicated and fuzzier than those influencing the macro transportation demand. Second, the PEM is chosen to find the implied periods from an observed data series because seasonality is related to the economic cycle but does not coincide with the cycle.

The remnant correction model is used to improve the reliability and validation of the main model. When the same model is applied to different objects (i.e. different observed data series), remnants will be generated, because each market has different influencing factors and changing laws, including trends, period length and fluctuation range. The remnants, which are usually included as the last component of a model, as random perturbation, are often referred to as white noise. The fact is that remnants are related to the reliability and validity of the model. When remnants are large and cannot be ignored, the model should be adjusted until satisfactory reliability and validity are achieved. In this article, prediction accuracy analysis, by checking the remnants of the GM(1,1)-PEM, is performed to determine whether the main model needs adjustment. If the main model needs fixing, GM(1,1) will be used to fit the remnant series, fixing the main model.

The proposed model shows good adaptability in the case study in terms of the choice of the sample company and sample data series. This sample company is a large-scale transportation-oriented logistics company, whose business is an LTL road transportation service in a nationwide hub-and-spoke network in China. The LTL transportation market in China is now clearly divided into two levels. The upper level is a nationwide business, which is taken up by a few large-scale companies similar to the sample one, while numerous Small and Medium-sized Enterprises (SMEs) compete in the remaining market. Competition is fierce, and concentration is very low. As the level of logistics management has increased in recent years, market concentration has slowly grown. The sample company has a good representation of the upper market segmentation. The observed time series cover all the business data of the three origin-destination hub cities of the sample company during a 36-month period (January 2009 to December 2011), during which the national economy had only just recovered from the financial crisis. The forecasting results show that the model has good adaptability to the macroeconomic environment.