# Grain Consumption Forecasting: One Modified MLR Model Combined with Time Series Forecasting Theory

Authors: Chunhua Zhu; Jiaojiao Wang; Jiake Tian
Journal: International Journal of Mathematics and Mathematical Sciences
 (2020)

---

## Abstract

In the classical multivariate prediction model, most research studies focused on the selection of relevant behaviour factors and the stability of historical data for improving the predicting accuracy of the main behaviour factor, and the historical data of the main behaviour factor have never been considered as one relevant behaviour factor, which in fact can be the first key impact factor; besides, the historical data can directly predict the main behaviour in the time series forecasting model, such as the ARIMA model. In this paper, one modified MLR model combined with time series forecasting theory is presented and applied in grain consumption forecasting. In the proposed model, to improve the current grain consumption forecasting, how to select impact factors is also discussed by combining the grey relational degree and Pearson correlation coefficient with given weights, and the optimal preprocessing parameter by the moving average filtering is computed for eliminating the abnormal points and stabilizing the data. Finally, the selected main impact factors are inputted to the proposed modified MLR model for forecasting grain consumption. Simulation results have shown that the five-year mean absolute percentage error of ration and feed grain is 2.34% and 3.27%, respectively, and the prediction accuracy has improved up to 2 times compared with the BP model and LSTM model. Moreover, the robustness of the model is verified by prediction analysis at different time intervals of historical data.


---

## Body

## 1. Introduction

With the change of social, economic, and environmental factors, the grain consumption structure appears to have new features. There have been various prediction methods by research methods, focuses, and perspectives, which can be divided into qualitative analysis and quantitative analysis. Qualitative prediction is the analysis of the specialists based on the characteristics of historical data and intuitive materials, and the corresponding results rely on the experience and analysis ability of researchers. The quantitative prediction method is to build one mathematical model to predict the changing laws and the development trend in future, which has obtained wider acceptance and application. In general, the quantitative analysis methods include the time series prediction model, econometric equation (single equation and simultaneous equation), regression analysis prediction model, etc [1]. Cheng et al. considered the actual level of China’s current grain consumption and the future trend of China’s population development and predicted that the total grain consumption demand in China will not exceed 640 million tons in 2030 [2]. Gao used the time series method to calculate per capita grain consumption by selecting a national sample of residents and predicted total food consumption. He said that in 2020, China’s total grain consumption will reach to 595 million tons [3]. Chen et al. used the single equation econometric model to complete the prediction of per capita rural consumption and population changes in China; here the two variables were multiplied to obtain the total grain consumption, which will reach 76 million tons in 2020 and 65 million tons in 2030 [4]. Liao et al. proposed the modified CAPSIM-PODIUM model to forecast and analyze the grain demand of the whole country and nine major watershed slices, including rice, wheat, and corn; the results have shown that China’s total grain demand will reach 507 million tons in 2020 [5]. Pan et al. first established a regression analysis; here only two factors of food price and per capita disposable income are regarded as significant impacts on food consumption. Thus, the direct and indirect grain consumption models of rural residents per capita were obtained by the linear fitting model, and the predicted value of total grain consumption was obtained by combining the total population. The results have shown that direct and indirect grain consumption in China’s rural areas will reach 96.8658 million tons and 87.3173 million tons by 2020 and 58.1396 tons and 91.1647 million tons by 2030 [6]. From the existing research, the impact factors and historical data features are more important for the grain consumption prediction model. However, so far, it is not enough for the more comprehensive analysis. In this paper, one new combined model is proposed, including the selection of impact factors and the consideration of optimal preprocessing.The main contributions of this paper are listed as follows:(1)
To smooth original data and optimize the selection of impact factors, firstly, the abnormal points in the original data are removed by moving average filtering with the optimal window width, which makes the data trend to be smoothed and stabilized, and the relational degree of grain consumption and its impact factors are calculated by the proposed combined method, and thus the key impact factors are selected.(2)
To minimize the prediction error of grain consumption, we combine the advantages of two models, and one new modified MLR model combined with time series forecasting theory based on data barycenter is proposed.


## 2. Research of Consumption Forecasting Method

China’s grain consumption can be divided into two categories: food consumption and nonfood consumption. Food consumption includes ration and feed grain consumption, while nonfood consumption includes industrial grain, seed grain, and loss grain. As shown in Figure1, food consumption is the most important mode of grain consumption in China, so the research focus of this paper is the changing trend of ration and feed grain.Figure 1
The proportion of China’s grain consumption by category in 2017.


### 2.1. Moving Average Filtering with Optimal Window Span

Moving average filtering has been applied widely in numerical analysis, also known as linear filtering, which belongs to one-dimensional filtering. The algorithm is realized by recursion, that is, replacing the original value of the corresponding position with the average value of several neighborhood points, and moving average value is regarded as the trend representative value of the moving average term in this period. In this paper, all the historical data will be preprocessed by average moving filter for eliminating abnormal point and random fluctuations. Generally, when the filtering window widthN is a positive odd number, the general formula of moving average filtering is defined as(1)fn=xn+xn−1+xn−2+⋯+xn−N+1N,where xn is the input sequence, fn is the output sequence, and N is the window width.Under different window widths, the smoothness and stabilization of the filtering original data will have a large difference. Here, coefficient of variation is adopted to evaluate the discreteness of observed values, which can be expressed by the ratio of standard deviation to average value of each observed value, and the symbol is CV [7]:(2)CV=σμ,where σis the standard deviation, σ=1/n∑i=1nxi−μ2, and μ is the average value, μ=1/n∑i=1nxi. The smaller the coefficient of variation is, the better stability the group of data has. In fact, for obtaining higher prediction accuracy, the largest window width with the smallest coefficient of variation may not bring expected results; therefore, the optimal filtering window width will be chosen by the coefficient of variation and the prediction error.


### 2.2. Analysis of Impact Factors of Grain Consumption

#### 2.2.1. Grey Relational Degree

The method of grey correlation analysis, also called “grey relational degree,” is used widely to calculate the correlation between variables [8], which can analyze the degree of correlation between two systems or two factors within the system under the condition of the change of time and speed and judge whether two data series are highly related by comparing their geometric similarity. In this paper, assuming the grain consumption data are a reference sequence and their impact factors are comparative sequence, the grey relational degree between the reference sequence and comparative sequence is calculated as follows [9]:(1)
The reference sequence isx0k=x01,x02,…,x0n; comparative sequence is(3)xik=xi1,xi2,…,xin,i=1,2,…m.(2)
Normalization is adopted to carry on a dimensionless preprocessing, and the normalized reference and comparative series are expressed asz0k and zik, respectively.(4)z0k=x0kmaxx0k,zik=xikmaxxik.(3)
Calculating the absolute difference sequenceΔ0ik between reference sequence z0k and comparative sequence zik:(5)Δ0ik=z0k−zik=Δi1,Δi2,…,Δin.(4)
Searching the maximum absolute difference sequenceΔmax and minimum absolute difference sequence Δmin: the i-th grey relational coefficient is defined as(6)L0ik=ρΔmax+ΔminΔ0ik+ρΔmax,whereρ is the resolution coefficient, which is generally between 0.5 and 1. The smaller the parameter p is, the stronger the distinguishing ability for the difference among the relational coefficients is, usually taking 0.5.(5)
Averaging the relational coefficients: the relational degree can be obtained as(7)R0ik=1n∑k=1nL0ik=1nL0i1+L0i2+⋯+L0in.


#### 2.2.2. Pearson Correlation Coefficient

Pearson correlation coefficient is a classical statistic proposed by Karl Pearson. In statistics, it is widely used to measure the degree of correlation between two variables, which can reflect the linear correlation between two variables. It is defined as the quotient of covariance and standard deviation between two variables [10]. Assuming two series X=x1,x2,…,xn and Y=y1,y2,…,yn, the Pearson correlation coefficient between them is(8)r=n∑i=1nxiyi−∑i=1nxi∑i=1nyin∑i=1nxi2−∑i=1nxi2n∑i=1nyi2−∑i=1nyi2,where r is between −1 and 1. The larger the absolute value of the correlation coefficient is, the higher the correlation between X and Y is, and the positive correlation coefficient means positive correlation, and the negative correlation coefficient represents negative correlation.


#### 2.2.3. Combined Model of Selecting Main Impact Factors

Grain consumption includes ration consumption, feed-grain consumption, and industrial/seed/loss grain consumption. The changing social, economic, and environmental factors will produce a different effect on each component of the grain consumption; however, some uncorrelated impact factors and those with lower correlation will reduce the prediction accuracy and increase the complexity. Therefore, it is necessary to select main impact factors for the following prediction of the grain consumption.From Sections2.2.1 and 2.2.2, grey correlation analysis can embody the geometric similarity between any two factors, and Pearson correlation coefficient can reflect the linear correlation between two variables. How to combine the above two correlation degrees and select the real key impact factors will be one of innovations in this paper. In the existing references, only one correlation analysis usually is used to rank the correlation degree between one main behavioral factor and the reference series. Though for the same impact factors, there can produce the different ranking order of correlation degree by the above two correlation analysis methods [11]. One new combined model of selecting main impact factors is proposed in Figure 2. Here, grain consumption is the main behavioral factor, and its more possible impact factors are the reference series. From Figure 2, the correlation degree between them will firstly be computed separately by grey correlation analysis and Pearson correlation coefficient, and then one weighted operator will be adopted to obtain the final relational degree. Thus, the key impact factors can be chosen by the ranking order and one threshold of the final relational degree. Here, the threshold can be determined by the feedback of the following prediction error.Figure 2
Selection of impact factors.


#### 2.2.4. Prediction Model of Impact Factors

The ARIMA model is one time series prediction method proposed by Box and Jenkins in the early 1970s, which is also called the Box–Jenkins model [12, 13]. In the ARIMA (p, d, q) model, p is called autoregressive order and q is called moving average order. For obtaining better prediction, the original data should be stationary. The unit root test is often used to check the stationarity of time series. If the time series is nonstationary, it can be transformed into stationary time series by several times of difference [14].Suppose thatωt is the predictive value in t year and ωt−1,ωt−2,…,ωt−p are the original values of various impact factors in past p years and set ωt=1−Ldyt; here yt is a single integer sequence with d order and ωt is the stationary series; the general form of the ARMA model can be written as(9)ωt=φ1ωt−1+φ2ωt−2+⋯+φpωt−p+εt+θ1εt−1+⋯+θqεt−q.SupposeL is the lag operator, then(10)Lωt=ωt−1,Lpωt=ωt−p,and equation (9) can be rewritten as(11)φLωt=ΘLεt,where(12)φL=1−φ1L−φ2L2−⋯−φpLp,ΘL=1+θ1L+θ2L2+⋯+θqLq.Thed order difference transformation of ARMA(p, q) as equation (10) is called ARIMA(p, d, q):(13)φL1−Ldyt=ΘLεt,where εt is a white noise process with its mean value being 0 and variance being σ2.


### 2.3. MLR Prediction Model Based on Data Barycenter

In the classical multivariate linear regression (MLR) model, the least-square method is usually used to compute the prediction parameters. Compared with the least-square method, the barycenter operators [15] can provide higher stable parameters for a multivariate regression model [16].Assuming the dataset isxi,i=1,2,…,n, its first-order barycenter can be written as(14)x=x1+x2+x3+⋯+xnn=1n∑i=1nxi.If recorded∑i=1nxi=∑i=1nxi1, called 1/n∑i=1nxi1 a first-order barycenter operator. Similarly, its second-order barycenter can be written as(15)x=x1+x1+x2+⋅⋅⋅x1+x2+x3+⋅⋅⋅+xn1+2+3+⋅⋅⋅+n=∑i=1n∑i′=1ixi′∑i=1ni,where ∑i=1n∑i′=1ixi′=∑i=1n∑i′=1ix1i′=∑i=1ixi2, called 1/1+2+⋅⋅⋅+n∑i=1nxi2, a second-order barycenter operator.And so forth, thek-order barycenter operator is(16)∑i=1nxik1/k−1!∑i=1nn−i+1n−i+2⋅⋅⋅n−i+k−1,where ∑i=1nxik=1/k−1!∑i=1nn−i+1n−i+2⋅⋅⋅n−i+k−1xik≥2.In the classical multivariate linear regression model [17, 18], the prediction value can be expressed as(17)y1=β0+β1x11+β2x12+⋯+βp−1x1p−1,y2=β0+β1x21+β2x22+⋯+βp−1x2p−1,⋯yn=β0+β1xn1+β2xn2+⋯+βp−1xnp−1,here, β0,β1,…,βp−1 are the prediction parameters; x is the independent variables, and the number are p−1 of them, n sets of observations of variables are expressed asxi1,xi2,…,xip−1,yi,i=1,2,…,n.In this paper, thek-order barycenter operators are adopted to compute the parameters of multivariate prediction model. The order of k-order barycenter operator is equal to the number of the prediction model parameters. The new proposed model based on data barycenter is called DB-based MLR.In the DB-based MLR, equation (17) can be transformed as(18)∑i=1p−1YiK=β0∑i=1p−11K+β1∑i=1p−1xi1k+β2∑i=1p−1xi1k+⋅⋅⋅+βp−1∑i=1p−1xip−1k.Letk = 1, 2, …, p in equation (18); we can get one linear equation set. The corresponding parameters can be solved by Cramer's rule:(19)β^0=D1D,β^1=D2D,…,β^p−1=DpD,where(20)D=∑i=1n11∑i=1nxi11⋯∑i=1nxip−11∑i=1n12∑i=1nxi12⋯∑i=1nxip−12⋮⋮⋮⋮∑i=1n1p∑i=1nxi1p⋯∑i=1nxip−1p.Replacing thej-th column of D with ∑i=1nyi1,∑i=1nyi2,∑i=1nyip (that is, the left term in equation (18)), we can get the coefficients Dj, j=1,2,…,p. Furthermore, taking equation (19) into equation (17), the optimal predictive model parameters can be obtained by several iterations.


### 2.4. Time Series Prediction Model

The time series prediction model is one of the quantitative prediction methods, which arranges the historical data of the predicted object into a time series in time order, analyzes its changing trend with time, and establishes a mathematical model to extrapolate the future value. [19]. The general structure of time series prediction model is expressed as follows:(21)yt=ϕ1yt−1+ϕ2yt−2+ϕ3yt−3+⋯+ϕnyt−n,where ϕ1,ϕ2,…,ϕn are the prediction parameters and yt,yt−1,…,yt−n are the historical values of grain consumption. By the modeling idea and parameter calculating algorithms, the time series models can be divided into the moving average (MA) model, trend extrapolation model, exponential regression (ES) model, autoregressive (AR) model, ARMA model, ARIMA model, etc [20–22]. It is same for the above time series to predict the future value by fitting its historical data, only using the different computing methods of model parameters.


### 2.5. Modified MLR Model Combined with Time Series Forecasting Theory

Summarizing the analysis of the time series model and classical MLR model, we can find that the two models cannot make full use of the historical data and the impact factors at the same time. In this paper, the historical data of prediction variable will be considered as one key impact factor, which will be imputed to the MLR model with the other chosen impact factors. The proposed modified MLR model can be expressed as(22)yt=β0+β1x1t+β2x2t+⋯βnxnt+ϕ1yt−1+ϕ2yt−2+ϕ3yt−3+⋯+ϕnyt−n,where yt,yt−1,…,yt−n are the historical data of the prediction value yt; they are the “internal factors” of the MLR model; besides, x1t,x2t,…,xnt are the chosen key impact factors, which are “external factors.” Therefore, the modified MLR model can take advantages from adding the internal factors. In this model, the model parameters β0,β1,…,βn and ϕ1,ϕ2,…,ϕn can be computed at one time; therefore, the proposed model is different with the simple combining the two models, respectively. The prediction performance will be discussed in the following sections.


### 2.6. The Evaluation Index of Prediction Error

In this paper, the absolute percentage error (APE), the mean absolute percentage error (MAPE), and the Theil inequality coefficient of prediction residuals are used as the criteria for the model prediction. MAPE is the average ratio of all absolute percentage errors (APEs) to real values, and it is a percentage. In general, absolute error can avoid the problem of errors cancelling out each other; therefore, MAPE can accurately reflect the actual prediction error and can better reflect the reliability of prediction model. Besides, Theil inequality coefficient is a correlation index to measure the prediction accuracy. The closer the Theil inequality coefficient is to 0, the smaller the difference between the predicted value and the real value will be, which indicates the better fitting degree of the prediction model. The mean absolute percentage error and the Theil inequality coefficient are defined as(23)APE=yi−y^iyi×100%,MAPE=1n∑i=1nyi−y^iyi×100%,Theil’sU=1/n∑i=1nyi−y^i21/n∑i=1nyi2+1/n∑i=1ny^i2.


## 2.1. Moving Average Filtering with Optimal Window Span

Moving average filtering has been applied widely in numerical analysis, also known as linear filtering, which belongs to one-dimensional filtering. The algorithm is realized by recursion, that is, replacing the original value of the corresponding position with the average value of several neighborhood points, and moving average value is regarded as the trend representative value of the moving average term in this period. In this paper, all the historical data will be preprocessed by average moving filter for eliminating abnormal point and random fluctuations. Generally, when the filtering window widthN is a positive odd number, the general formula of moving average filtering is defined as(1)fn=xn+xn−1+xn−2+⋯+xn−N+1N,where xn is the input sequence, fn is the output sequence, and N is the window width.Under different window widths, the smoothness and stabilization of the filtering original data will have a large difference. Here, coefficient of variation is adopted to evaluate the discreteness of observed values, which can be expressed by the ratio of standard deviation to average value of each observed value, and the symbol is CV [7]:(2)CV=σμ,where σis the standard deviation, σ=1/n∑i=1nxi−μ2, and μ is the average value, μ=1/n∑i=1nxi. The smaller the coefficient of variation is, the better stability the group of data has. In fact, for obtaining higher prediction accuracy, the largest window width with the smallest coefficient of variation may not bring expected results; therefore, the optimal filtering window width will be chosen by the coefficient of variation and the prediction error.


## 2.2. Analysis of Impact Factors of Grain Consumption

### 2.2.1. Grey Relational Degree

The method of grey correlation analysis, also called “grey relational degree,” is used widely to calculate the correlation between variables [8], which can analyze the degree of correlation between two systems or two factors within the system under the condition of the change of time and speed and judge whether two data series are highly related by comparing their geometric similarity. In this paper, assuming the grain consumption data are a reference sequence and their impact factors are comparative sequence, the grey relational degree between the reference sequence and comparative sequence is calculated as follows [9]:(1)
The reference sequence isx0k=x01,x02,…,x0n; comparative sequence is(3)xik=xi1,xi2,…,xin,i=1,2,…m.(2)
Normalization is adopted to carry on a dimensionless preprocessing, and the normalized reference and comparative series are expressed asz0k and zik, respectively.(4)z0k=x0kmaxx0k,zik=xikmaxxik.(3)
Calculating the absolute difference sequenceΔ0ik between reference sequence z0k and comparative sequence zik:(5)Δ0ik=z0k−zik=Δi1,Δi2,…,Δin.(4)
Searching the maximum absolute difference sequenceΔmax and minimum absolute difference sequence Δmin: the i-th grey relational coefficient is defined as(6)L0ik=ρΔmax+ΔminΔ0ik+ρΔmax,whereρ is the resolution coefficient, which is generally between 0.5 and 1. The smaller the parameter p is, the stronger the distinguishing ability for the difference among the relational coefficients is, usually taking 0.5.(5)
Averaging the relational coefficients: the relational degree can be obtained as(7)R0ik=1n∑k=1nL0ik=1nL0i1+L0i2+⋯+L0in.


### 2.2.2. Pearson Correlation Coefficient

Pearson correlation coefficient is a classical statistic proposed by Karl Pearson. In statistics, it is widely used to measure the degree of correlation between two variables, which can reflect the linear correlation between two variables. It is defined as the quotient of covariance and standard deviation between two variables [10]. Assuming two series X=x1,x2,…,xn and Y=y1,y2,…,yn, the Pearson correlation coefficient between them is(8)r=n∑i=1nxiyi−∑i=1nxi∑i=1nyin∑i=1nxi2−∑i=1nxi2n∑i=1nyi2−∑i=1nyi2,where r is between −1 and 1. The larger the absolute value of the correlation coefficient is, the higher the correlation between X and Y is, and the positive correlation coefficient means positive correlation, and the negative correlation coefficient represents negative correlation.


### 2.2.3. Combined Model of Selecting Main Impact Factors

Grain consumption includes ration consumption, feed-grain consumption, and industrial/seed/loss grain consumption. The changing social, economic, and environmental factors will produce a different effect on each component of the grain consumption; however, some uncorrelated impact factors and those with lower correlation will reduce the prediction accuracy and increase the complexity. Therefore, it is necessary to select main impact factors for the following prediction of the grain consumption.From Sections2.2.1 and 2.2.2, grey correlation analysis can embody the geometric similarity between any two factors, and Pearson correlation coefficient can reflect the linear correlation between two variables. How to combine the above two correlation degrees and select the real key impact factors will be one of innovations in this paper. In the existing references, only one correlation analysis usually is used to rank the correlation degree between one main behavioral factor and the reference series. Though for the same impact factors, there can produce the different ranking order of correlation degree by the above two correlation analysis methods [11]. One new combined model of selecting main impact factors is proposed in Figure 2. Here, grain consumption is the main behavioral factor, and its more possible impact factors are the reference series. From Figure 2, the correlation degree between them will firstly be computed separately by grey correlation analysis and Pearson correlation coefficient, and then one weighted operator will be adopted to obtain the final relational degree. Thus, the key impact factors can be chosen by the ranking order and one threshold of the final relational degree. Here, the threshold can be determined by the feedback of the following prediction error.Figure 2
Selection of impact factors.


### 2.2.4. Prediction Model of Impact Factors

The ARIMA model is one time series prediction method proposed by Box and Jenkins in the early 1970s, which is also called the Box–Jenkins model [12, 13]. In the ARIMA (p, d, q) model, p is called autoregressive order and q is called moving average order. For obtaining better prediction, the original data should be stationary. The unit root test is often used to check the stationarity of time series. If the time series is nonstationary, it can be transformed into stationary time series by several times of difference [14].Suppose thatωt is the predictive value in t year and ωt−1,ωt−2,…,ωt−p are the original values of various impact factors in past p years and set ωt=1−Ldyt; here yt is a single integer sequence with d order and ωt is the stationary series; the general form of the ARMA model can be written as(9)ωt=φ1ωt−1+φ2ωt−2+⋯+φpωt−p+εt+θ1εt−1+⋯+θqεt−q.SupposeL is the lag operator, then(10)Lωt=ωt−1,Lpωt=ωt−p,and equation (9) can be rewritten as(11)φLωt=ΘLεt,where(12)φL=1−φ1L−φ2L2−⋯−φpLp,ΘL=1+θ1L+θ2L2+⋯+θqLq.Thed order difference transformation of ARMA(p, q) as equation (10) is called ARIMA(p, d, q):(13)φL1−Ldyt=ΘLεt,where εt is a white noise process with its mean value being 0 and variance being σ2.


## 2.2.1. Grey Relational Degree

The method of grey correlation analysis, also called “grey relational degree,” is used widely to calculate the correlation between variables [8], which can analyze the degree of correlation between two systems or two factors within the system under the condition of the change of time and speed and judge whether two data series are highly related by comparing their geometric similarity. In this paper, assuming the grain consumption data are a reference sequence and their impact factors are comparative sequence, the grey relational degree between the reference sequence and comparative sequence is calculated as follows [9]:(1)
The reference sequence isx0k=x01,x02,…,x0n; comparative sequence is(3)xik=xi1,xi2,…,xin,i=1,2,…m.(2)
Normalization is adopted to carry on a dimensionless preprocessing, and the normalized reference and comparative series are expressed asz0k and zik, respectively.(4)z0k=x0kmaxx0k,zik=xikmaxxik.(3)
Calculating the absolute difference sequenceΔ0ik between reference sequence z0k and comparative sequence zik:(5)Δ0ik=z0k−zik=Δi1,Δi2,…,Δin.(4)
Searching the maximum absolute difference sequenceΔmax and minimum absolute difference sequence Δmin: the i-th grey relational coefficient is defined as(6)L0ik=ρΔmax+ΔminΔ0ik+ρΔmax,whereρ is the resolution coefficient, which is generally between 0.5 and 1. The smaller the parameter p is, the stronger the distinguishing ability for the difference among the relational coefficients is, usually taking 0.5.(5)
Averaging the relational coefficients: the relational degree can be obtained as(7)R0ik=1n∑k=1nL0ik=1nL0i1+L0i2+⋯+L0in.


## 2.2.2. Pearson Correlation Coefficient

Pearson correlation coefficient is a classical statistic proposed by Karl Pearson. In statistics, it is widely used to measure the degree of correlation between two variables, which can reflect the linear correlation between two variables. It is defined as the quotient of covariance and standard deviation between two variables [10]. Assuming two series X=x1,x2,…,xn and Y=y1,y2,…,yn, the Pearson correlation coefficient between them is(8)r=n∑i=1nxiyi−∑i=1nxi∑i=1nyin∑i=1nxi2−∑i=1nxi2n∑i=1nyi2−∑i=1nyi2,where r is between −1 and 1. The larger the absolute value of the correlation coefficient is, the higher the correlation between X and Y is, and the positive correlation coefficient means positive correlation, and the negative correlation coefficient represents negative correlation.


## 2.2.3. Combined Model of Selecting Main Impact Factors

Grain consumption includes ration consumption, feed-grain consumption, and industrial/seed/loss grain consumption. The changing social, economic, and environmental factors will produce a different effect on each component of the grain consumption; however, some uncorrelated impact factors and those with lower correlation will reduce the prediction accuracy and increase the complexity. Therefore, it is necessary to select main impact factors for the following prediction of the grain consumption.From Sections2.2.1 and 2.2.2, grey correlation analysis can embody the geometric similarity between any two factors, and Pearson correlation coefficient can reflect the linear correlation between two variables. How to combine the above two correlation degrees and select the real key impact factors will be one of innovations in this paper. In the existing references, only one correlation analysis usually is used to rank the correlation degree between one main behavioral factor and the reference series. Though for the same impact factors, there can produce the different ranking order of correlation degree by the above two correlation analysis methods [11]. One new combined model of selecting main impact factors is proposed in Figure 2. Here, grain consumption is the main behavioral factor, and its more possible impact factors are the reference series. From Figure 2, the correlation degree between them will firstly be computed separately by grey correlation analysis and Pearson correlation coefficient, and then one weighted operator will be adopted to obtain the final relational degree. Thus, the key impact factors can be chosen by the ranking order and one threshold of the final relational degree. Here, the threshold can be determined by the feedback of the following prediction error.Figure 2
Selection of impact factors.


## 2.2.4. Prediction Model of Impact Factors

The ARIMA model is one time series prediction method proposed by Box and Jenkins in the early 1970s, which is also called the Box–Jenkins model [12, 13]. In the ARIMA (p, d, q) model, p is called autoregressive order and q is called moving average order. For obtaining better prediction, the original data should be stationary. The unit root test is often used to check the stationarity of time series. If the time series is nonstationary, it can be transformed into stationary time series by several times of difference [14].Suppose thatωt is the predictive value in t year and ωt−1,ωt−2,…,ωt−p are the original values of various impact factors in past p years and set ωt=1−Ldyt; here yt is a single integer sequence with d order and ωt is the stationary series; the general form of the ARMA model can be written as(9)ωt=φ1ωt−1+φ2ωt−2+⋯+φpωt−p+εt+θ1εt−1+⋯+θqεt−q.SupposeL is the lag operator, then(10)Lωt=ωt−1,Lpωt=ωt−p,and equation (9) can be rewritten as(11)φLωt=ΘLεt,where(12)φL=1−φ1L−φ2L2−⋯−φpLp,ΘL=1+θ1L+θ2L2+⋯+θqLq.Thed order difference transformation of ARMA(p, q) as equation (10) is called ARIMA(p, d, q):(13)φL1−Ldyt=ΘLεt,where εt is a white noise process with its mean value being 0 and variance being σ2.


## 2.3. MLR Prediction Model Based on Data Barycenter

In the classical multivariate linear regression (MLR) model, the least-square method is usually used to compute the prediction parameters. Compared with the least-square method, the barycenter operators [15] can provide higher stable parameters for a multivariate regression model [16].Assuming the dataset isxi,i=1,2,…,n, its first-order barycenter can be written as(14)x=x1+x2+x3+⋯+xnn=1n∑i=1nxi.If recorded∑i=1nxi=∑i=1nxi1, called 1/n∑i=1nxi1 a first-order barycenter operator. Similarly, its second-order barycenter can be written as(15)x=x1+x1+x2+⋅⋅⋅x1+x2+x3+⋅⋅⋅+xn1+2+3+⋅⋅⋅+n=∑i=1n∑i′=1ixi′∑i=1ni,where ∑i=1n∑i′=1ixi′=∑i=1n∑i′=1ix1i′=∑i=1ixi2, called 1/1+2+⋅⋅⋅+n∑i=1nxi2, a second-order barycenter operator.And so forth, thek-order barycenter operator is(16)∑i=1nxik1/k−1!∑i=1nn−i+1n−i+2⋅⋅⋅n−i+k−1,where ∑i=1nxik=1/k−1!∑i=1nn−i+1n−i+2⋅⋅⋅n−i+k−1xik≥2.In the classical multivariate linear regression model [17, 18], the prediction value can be expressed as(17)y1=β0+β1x11+β2x12+⋯+βp−1x1p−1,y2=β0+β1x21+β2x22+⋯+βp−1x2p−1,⋯yn=β0+β1xn1+β2xn2+⋯+βp−1xnp−1,here, β0,β1,…,βp−1 are the prediction parameters; x is the independent variables, and the number are p−1 of them, n sets of observations of variables are expressed asxi1,xi2,…,xip−1,yi,i=1,2,…,n.In this paper, thek-order barycenter operators are adopted to compute the parameters of multivariate prediction model. The order of k-order barycenter operator is equal to the number of the prediction model parameters. The new proposed model based on data barycenter is called DB-based MLR.In the DB-based MLR, equation (17) can be transformed as(18)∑i=1p−1YiK=β0∑i=1p−11K+β1∑i=1p−1xi1k+β2∑i=1p−1xi1k+⋅⋅⋅+βp−1∑i=1p−1xip−1k.Letk = 1, 2, …, p in equation (18); we can get one linear equation set. The corresponding parameters can be solved by Cramer's rule:(19)β^0=D1D,β^1=D2D,…,β^p−1=DpD,where(20)D=∑i=1n11∑i=1nxi11⋯∑i=1nxip−11∑i=1n12∑i=1nxi12⋯∑i=1nxip−12⋮⋮⋮⋮∑i=1n1p∑i=1nxi1p⋯∑i=1nxip−1p.Replacing thej-th column of D with ∑i=1nyi1,∑i=1nyi2,∑i=1nyip (that is, the left term in equation (18)), we can get the coefficients Dj, j=1,2,…,p. Furthermore, taking equation (19) into equation (17), the optimal predictive model parameters can be obtained by several iterations.


## 2.4. Time Series Prediction Model

The time series prediction model is one of the quantitative prediction methods, which arranges the historical data of the predicted object into a time series in time order, analyzes its changing trend with time, and establishes a mathematical model to extrapolate the future value. [19]. The general structure of time series prediction model is expressed as follows:(21)yt=ϕ1yt−1+ϕ2yt−2+ϕ3yt−3+⋯+ϕnyt−n,where ϕ1,ϕ2,…,ϕn are the prediction parameters and yt,yt−1,…,yt−n are the historical values of grain consumption. By the modeling idea and parameter calculating algorithms, the time series models can be divided into the moving average (MA) model, trend extrapolation model, exponential regression (ES) model, autoregressive (AR) model, ARMA model, ARIMA model, etc [20–22]. It is same for the above time series to predict the future value by fitting its historical data, only using the different computing methods of model parameters.


## 2.5. Modified MLR Model Combined with Time Series Forecasting Theory

Summarizing the analysis of the time series model and classical MLR model, we can find that the two models cannot make full use of the historical data and the impact factors at the same time. In this paper, the historical data of prediction variable will be considered as one key impact factor, which will be imputed to the MLR model with the other chosen impact factors. The proposed modified MLR model can be expressed as(22)yt=β0+β1x1t+β2x2t+⋯βnxnt+ϕ1yt−1+ϕ2yt−2+ϕ3yt−3+⋯+ϕnyt−n,where yt,yt−1,…,yt−n are the historical data of the prediction value yt; they are the “internal factors” of the MLR model; besides, x1t,x2t,…,xnt are the chosen key impact factors, which are “external factors.” Therefore, the modified MLR model can take advantages from adding the internal factors. In this model, the model parameters β0,β1,…,βn and ϕ1,ϕ2,…,ϕn can be computed at one time; therefore, the proposed model is different with the simple combining the two models, respectively. The prediction performance will be discussed in the following sections.


## 2.6. The Evaluation Index of Prediction Error

In this paper, the absolute percentage error (APE), the mean absolute percentage error (MAPE), and the Theil inequality coefficient of prediction residuals are used as the criteria for the model prediction. MAPE is the average ratio of all absolute percentage errors (APEs) to real values, and it is a percentage. In general, absolute error can avoid the problem of errors cancelling out each other; therefore, MAPE can accurately reflect the actual prediction error and can better reflect the reliability of prediction model. Besides, Theil inequality coefficient is a correlation index to measure the prediction accuracy. The closer the Theil inequality coefficient is to 0, the smaller the difference between the predicted value and the real value will be, which indicates the better fitting degree of the prediction model. The mean absolute percentage error and the Theil inequality coefficient are defined as(23)APE=yi−y^iyi×100%,MAPE=1n∑i=1nyi−y^iyi×100%,Theil’sU=1/n∑i=1nyi−y^i21/n∑i=1nyi2+1/n∑i=1ny^i2.


## 3. Simulation Analysis

### 3.1. Moving Average Filtering with Optimal Window Width

In experiment, the data are from the website of the National Statistical Office and “China Statistical Yearbook 2018” [23]. All impact factors and different kinds of grain consumptions will be preprocessed by the moving average filtering under different window widths, as shown in Figures 3 and 4. Here, urban ration (1981∼2017) is one example of the original data in the experiment.Figure 3
Filtering result of urban ration.Figure 4
The coefficient of variation test.From Figures3 and 4, for different window width N, the larger the window width is, the smaller the coefficient of variation is and the smoother the filtered sequence is. Considering the following prediction accuracy, for given training data, the optimal window width can be determined by two indexes: coefficient of variation test and prediction error.


### 3.2. Selection of Impact Factors

As shown in Sections2.2.1 and 2.2.2, the relational degree and order between each kind of grain consumption and the corresponding several impact factors are analyzed. The selection of key impact factors is shown as Table 1. Here, the grain consumption includes urban ration, rural ration, urban feed grain, and rural feed grain. The impact factors include population, urbanization level, Engel’s coefficient, per capita income, price index of agricultural products, etc.Table 1
Selection of key impact factors by the proposed combined model.
Classification of grain consumption
Key impact factors
x
11
x
22
Urban ration
Urban population
Urban per capita income
Rural ration
Rural population
Rural Engel’s coefficient
Urban feed grain
Urban population
Urbanization level
Rural feed grain
Urbanization level
Rural Engel’s coefficient


### 3.3. Analysis of Model Prediction

In the modified prediction model proposed above, by considering the complexity of the algorithm and the comparison of experimental results, the final model determined by the project is as follows:(24)yt=β0+β1x1t+β2x2t+ϕ1yt−1,where model parameters β0,β1,β2, and ϕ1 can be calculated by the data barycenter method.The original data include grain consumption and several impact factors. In simulation, the first 33 years (1981–2013) of grain consumption data are trained in the model, and then the next five years (2014–2017) of grain consumption data are predicted. Defining the year of 1981 ast = 1 and 2013 as t = 33, the urban ration consumption, urban feed-grain consumption, rural ration consumption, and rural feed-grain consumption are, respectively, expressed as y1,y2,y3, and y4; the five factors are, respectively, defined as x1,x2,x3,x4, and x5. The forecasting process is shown as follows:(1)
Whent = 1∼37, the relational degree between y1 and x1,x2,x3,x4, and x5 is calculated; the two key impact factors will be selected by the proposed combined method, expressed as x11t and x22t.(2)
Whent = 1∼33, y1t, x11t, x22t, and y1t−1 are preprocessed by average moving filtering and are expressed as y1′t, x11′t, x22′t, and y1′t−1, respectively, which are used as training data into y1′t=β0+β1x11′t+β2x22′t+ϕ1y1′t−1. The model parameters β0,β1,β2, and β3 will calculated by data barycenter, and then the prediction model can be obtained.(3)
x11t+1 and x22t+1 are predicted by the ARIMA model.(4)
Assumingt = t + 1 and y1t+1=β0+β1x11t+1+β2x22t+1+ϕ1y0t, the urban ration consumption next year is obtained by the proposed joint model.(5)
Repeating (3)∼(4), the prediction of the urban ration consumption can be completed.Similar to the forecasting process of rural ration consumption, the urban feed-grain consumptiony2, rural ration consumption y3, and rural feed-grain consumption y4 can be obtained, respectively. For obtaining the total grain consumption, we can sum the above prediction results.


#### 3.3.1. Fitting Results

By the proposed method as equation (24), for every kind of grain consumption, two key impact factors and the historical grain consumption are chosen as the input of the time series-MLR joint model. However, for the classical MLR model, the input is only the two chosen key impact factors. The fitting results of the two models are shown in Figure 5(a) for feed grain prediction and Figure 5(b) for ration grain prediction.Figure 5
The fitting results of different models (1981–2013). (a) Feed grain consumption. (b) Ration consumption.
(a)
(b)From Table2, it can be seen that the fitting goodness of ration and feed grain all can be improved. The fitting results of the proposed model can provide better prediction reliability.Table 2
The fitting goodness of grain consumption trend.
Grain consumption
Ration
Feed grain
Predicted model
Multivariate regression
Joint model
Multivariate regression
Joint model
R
2
0.9656
0.9863
0.9917
0.9965


#### 3.3.2. Analysis of Prediction Results

For the experiments of predicting the future trend of grain consumption, we choose the filter window width equal to 3, 5, or 7 and establish the prediction model. Different filtering window widths will affect the model prediction performance; therefore, if we feed the prediction error back to adjust the filter parameters, the optimal filtering window width can be chosen. Here, the chosen optimal window width is 5 and 3, respectively, for ration consumption prediction and feed grain consumption prediction. In this paper, the MLR model and modified MLR model are established based on these optimal preprocessing methods.In addition, in our simulation, a classical time series ARIMA model is adopted to predict the ration and feed grain consumption in China. Besides, deep learning prediction is also considered for comparison, which is the emerging research hotspot in recent years [24–30]. Five kinds of prediction models are considered for comparing APE, MAPE, and Theil’s U performance.(1)
The classical time series ARIMA model:(25)yt=φ1yt−1+φ2yt−2+⋯+φpyt−p+εt+θ1εt−1+⋯+θqεt−q,whereyt−1,yt−2,…,yt−p is the actual value of the past p years; the order of p and q can be realized by model recognition and φ1,φ2,…,φp and θ1,θ2,…,θq can be calculated by the least-square method.(2)
The MLR model:The two impact factors with the greatest relational degree as the key impact factors in the modeling were input into the MLR prediction model:(26)yt=α0+α1x1t+α2x2t,wherex1t,x2t are key impact factors; model parameters α0,α1,α2 can be calculated by the data barycenter method.(3)
The modified MLR model:Considering the complexity of the algorithm and the prediction performance, the proposed MLR model is(27)yt=β0+β1x1t+β2x2t+ϕ1yt−1,where model parametersβ0,β1,β2, and ϕ1 can be calculated by the data barycenter method; yt is the future prediction value; yt−1 is the historical data of the prediction variable; and x1t,x2t are the chosen optimal impact factors.(4)
The BP model:Backpropagation (BP) is one of the artificial neural networks (ANNs) applied widely, which can train the model parameters by reversing the error propagation algorithm [24–28]. The prediction frame of the BP model is shown in Figure 6.There are three layers to predict the future value from the three inputted history data; here the hidden layer is used to analyze the characteristics of input data. In Figure6, xt−2,xt−1,xt is the inputted history data of grain consumption, yt is the future prediction value, and W1 and W2 are the weighting values of the trained model.(5)
The LSTM model:LSTM is one special recurrent neural network (RNN), which substitutes the original neural units of hidden layers for memory units [29, 30]. Each memory unit has an input door, output door, forgotten door, and memory cells, as shown in Figure 7.Figure 6
The prediction frame of the BP model.Figure 7
The prediction frame of the LSTM model.Herext−2,xt−1,xt is the input history data of grain consumption, ct−2,ct−2,ct−1,ct is the output of different memory unit, yt−3,yt−2,yt−1 is the prediction value of each layer, and so on; we can get the final future value yt.The comparison of the prediction results of these five models is shown in Tables3 and 4.Table 3
Comparison of ration consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Year
2014
2015
2016
2017
2018
Original ration consumption
25532
24411
24192
23784
23361
MLR model
Predicted value
23081
23437
23436
23314
23243
APE
9.60%
3.99%
3.13%
1.98%
0.51%
MAPE
3.84%
Theil’sU
0.0262
ARIMA model
Predicted value
25410
25647
25678
25831
25831
APE
0.48%
5.06%
6.14%
8.61%
10.57%
MAPE
6.17%
Theil’sU
0.0336
BP model
Predicted value
25760
28837
26097
28026
26489
APE
0.89%
18.13%
7.87%
17.84%
13.39%
MAPE
11.62%
Theil’sU
0.0623
LSTM model
Predicted value
22424
22861
24273
22158
22928
APE
12.17%
6.35%
0.33%
6.84%
1.85%
MAPE
5.51%
Theil’sU
0.0366
Modified MLR model
Predicted value
24029
24031
24205
24067
24059
APE
5.89%
1.56%
0.05%
1.19%
2.99%
MAPE
2.34%
Theil’sU
0.0159Table 4
Comparison of feed grain consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Year
2014
2015
2016
2017
2018
Original feed grain consumption
22247
23047
23317
23920
25300
MLR model
Predicted value
23921
24155
24286
24393
24483
APE
7.52%
4.81%
4.16%
1.98%
3.23%
MAPE
4.34%
Theil’sU
0.0226
ARIMA model
Predicted value
21619
21703
21730
21730
21730
APE
2.82%
5.83%
6.81%
9.16%
14.11%
MAPE
7.75%
Theil’sU
0.0466
BP model
Predicted value
23145
23233
26432
25253
25977
APE
4.04%
0.81%
13.36%
5.57%
2.68%
MAPE
5.29%
Theil’sU
0.0330
LSTM model
Predicted value
24378
24876
24956
27131
27481
APE
9.58%
7.94%
7.03%
13.42%
8.62%
MAPE
9.32%
Theil’sU
0.0459
Modified MLR model
Predicted value
23056
23100
23125
23146
23163
APE
3.64%
0.23%
0.82%
3.24%
8.45%
MAPE
3.27%
Theil’sU
0.0232From Tables3 and 4, for the ARIMA model, the closer the forecast year is, the higher the prediction accuracy of the model is. For the MLR model, the further away the predicted year is, the higher the prediction accuracy of the model is. This is because the time series prediction model is more suitable for short-term prediction and the multiple regression model is suitable for long-term prediction. However, it can be found that the prediction accuracy of the proposed modified MLR model is higher than that of the single time series model and the classical MLR model. It is mainly because the combined model can not only track the change of impact factors but also grasp the internal change trend of grain consumption, thereby leading to a small error range. And the proposed model is superior to the BP model and LSTM model for the two indexes of MAPE and Theil’s U.Meanwhile, in order to verify the stability of the modified MLR model, two experiments are carried out: (1) the training time interval is 1981–2014, and the predicting time interval is 2015–2018; (2) the training time interval is 1981–2013, and the predicting time interval is 2014–2017. The prediction performance of ration consumption and feed grain consumption is shown in Table5.Table 5
Comparison of feed grain consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Training time interval (predicting time interval)
1981–2014 (2015–2018)
1981–2013 (2014–2017)
Year
2015
2016
2017
2018
2014
2015
2016
2017
Ration consumption
Actual value
24411
24192
23784
23361
25532
24411
24192
23784
Predicted value
24142
23977
23948
23985
24029
24031
24205
24067
APE
1.10%
0.89%
0.69%
2.67%
5.89%
1.56%
0.05%
1.19%
RMSE
365.67
787.79
MAPE
1.34%
2.17%
Theil’sU
0.0076
0.0162
Feed grain consumption
Actual value
23047
23317
23920
25300
22247
23047
23317
23920
Predicted value
22982
22994
23006
22999
23056
23100
23125
23146
APE
0.28%
1.39%
3.82%
9.09%
3.64%
0.23%
0.82%
3.24%
RMSE
1248.85
568.38
MAPE
3.64%
1.98%
Theil’sU
0.0266
0.0123From Table5, it can be seen that the prediction error limit is kept within 2% for ration and feed grain consumption under different time intervals, which indicates that the proposed model can realize stable prediction with satisfying accuracy for the medium-short-term prediction.


## 3.1. Moving Average Filtering with Optimal Window Width

In experiment, the data are from the website of the National Statistical Office and “China Statistical Yearbook 2018” [23]. All impact factors and different kinds of grain consumptions will be preprocessed by the moving average filtering under different window widths, as shown in Figures 3 and 4. Here, urban ration (1981∼2017) is one example of the original data in the experiment.Figure 3
Filtering result of urban ration.Figure 4
The coefficient of variation test.From Figures3 and 4, for different window width N, the larger the window width is, the smaller the coefficient of variation is and the smoother the filtered sequence is. Considering the following prediction accuracy, for given training data, the optimal window width can be determined by two indexes: coefficient of variation test and prediction error.


## 3.2. Selection of Impact Factors

As shown in Sections2.2.1 and 2.2.2, the relational degree and order between each kind of grain consumption and the corresponding several impact factors are analyzed. The selection of key impact factors is shown as Table 1. Here, the grain consumption includes urban ration, rural ration, urban feed grain, and rural feed grain. The impact factors include population, urbanization level, Engel’s coefficient, per capita income, price index of agricultural products, etc.Table 1
Selection of key impact factors by the proposed combined model.
Classification of grain consumption
Key impact factors
x
11
x
22
Urban ration
Urban population
Urban per capita income
Rural ration
Rural population
Rural Engel’s coefficient
Urban feed grain
Urban population
Urbanization level
Rural feed grain
Urbanization level
Rural Engel’s coefficient


## 3.3. Analysis of Model Prediction

In the modified prediction model proposed above, by considering the complexity of the algorithm and the comparison of experimental results, the final model determined by the project is as follows:(24)yt=β0+β1x1t+β2x2t+ϕ1yt−1,where model parameters β0,β1,β2, and ϕ1 can be calculated by the data barycenter method.The original data include grain consumption and several impact factors. In simulation, the first 33 years (1981–2013) of grain consumption data are trained in the model, and then the next five years (2014–2017) of grain consumption data are predicted. Defining the year of 1981 ast = 1 and 2013 as t = 33, the urban ration consumption, urban feed-grain consumption, rural ration consumption, and rural feed-grain consumption are, respectively, expressed as y1,y2,y3, and y4; the five factors are, respectively, defined as x1,x2,x3,x4, and x5. The forecasting process is shown as follows:(1)
Whent = 1∼37, the relational degree between y1 and x1,x2,x3,x4, and x5 is calculated; the two key impact factors will be selected by the proposed combined method, expressed as x11t and x22t.(2)
Whent = 1∼33, y1t, x11t, x22t, and y1t−1 are preprocessed by average moving filtering and are expressed as y1′t, x11′t, x22′t, and y1′t−1, respectively, which are used as training data into y1′t=β0+β1x11′t+β2x22′t+ϕ1y1′t−1. The model parameters β0,β1,β2, and β3 will calculated by data barycenter, and then the prediction model can be obtained.(3)
x11t+1 and x22t+1 are predicted by the ARIMA model.(4)
Assumingt = t + 1 and y1t+1=β0+β1x11t+1+β2x22t+1+ϕ1y0t, the urban ration consumption next year is obtained by the proposed joint model.(5)
Repeating (3)∼(4), the prediction of the urban ration consumption can be completed.Similar to the forecasting process of rural ration consumption, the urban feed-grain consumptiony2, rural ration consumption y3, and rural feed-grain consumption y4 can be obtained, respectively. For obtaining the total grain consumption, we can sum the above prediction results.


### 3.3.1. Fitting Results

By the proposed method as equation (24), for every kind of grain consumption, two key impact factors and the historical grain consumption are chosen as the input of the time series-MLR joint model. However, for the classical MLR model, the input is only the two chosen key impact factors. The fitting results of the two models are shown in Figure 5(a) for feed grain prediction and Figure 5(b) for ration grain prediction.Figure 5
The fitting results of different models (1981–2013). (a) Feed grain consumption. (b) Ration consumption.
(a)
(b)From Table2, it can be seen that the fitting goodness of ration and feed grain all can be improved. The fitting results of the proposed model can provide better prediction reliability.Table 2
The fitting goodness of grain consumption trend.
Grain consumption
Ration
Feed grain
Predicted model
Multivariate regression
Joint model
Multivariate regression
Joint model
R
2
0.9656
0.9863
0.9917
0.9965


### 3.3.2. Analysis of Prediction Results

For the experiments of predicting the future trend of grain consumption, we choose the filter window width equal to 3, 5, or 7 and establish the prediction model. Different filtering window widths will affect the model prediction performance; therefore, if we feed the prediction error back to adjust the filter parameters, the optimal filtering window width can be chosen. Here, the chosen optimal window width is 5 and 3, respectively, for ration consumption prediction and feed grain consumption prediction. In this paper, the MLR model and modified MLR model are established based on these optimal preprocessing methods.In addition, in our simulation, a classical time series ARIMA model is adopted to predict the ration and feed grain consumption in China. Besides, deep learning prediction is also considered for comparison, which is the emerging research hotspot in recent years [24–30]. Five kinds of prediction models are considered for comparing APE, MAPE, and Theil’s U performance.(1)
The classical time series ARIMA model:(25)yt=φ1yt−1+φ2yt−2+⋯+φpyt−p+εt+θ1εt−1+⋯+θqεt−q,whereyt−1,yt−2,…,yt−p is the actual value of the past p years; the order of p and q can be realized by model recognition and φ1,φ2,…,φp and θ1,θ2,…,θq can be calculated by the least-square method.(2)
The MLR model:The two impact factors with the greatest relational degree as the key impact factors in the modeling were input into the MLR prediction model:(26)yt=α0+α1x1t+α2x2t,wherex1t,x2t are key impact factors; model parameters α0,α1,α2 can be calculated by the data barycenter method.(3)
The modified MLR model:Considering the complexity of the algorithm and the prediction performance, the proposed MLR model is(27)yt=β0+β1x1t+β2x2t+ϕ1yt−1,where model parametersβ0,β1,β2, and ϕ1 can be calculated by the data barycenter method; yt is the future prediction value; yt−1 is the historical data of the prediction variable; and x1t,x2t are the chosen optimal impact factors.(4)
The BP model:Backpropagation (BP) is one of the artificial neural networks (ANNs) applied widely, which can train the model parameters by reversing the error propagation algorithm [24–28]. The prediction frame of the BP model is shown in Figure 6.There are three layers to predict the future value from the three inputted history data; here the hidden layer is used to analyze the characteristics of input data. In Figure6, xt−2,xt−1,xt is the inputted history data of grain consumption, yt is the future prediction value, and W1 and W2 are the weighting values of the trained model.(5)
The LSTM model:LSTM is one special recurrent neural network (RNN), which substitutes the original neural units of hidden layers for memory units [29, 30]. Each memory unit has an input door, output door, forgotten door, and memory cells, as shown in Figure 7.Figure 6
The prediction frame of the BP model.Figure 7
The prediction frame of the LSTM model.Herext−2,xt−1,xt is the input history data of grain consumption, ct−2,ct−2,ct−1,ct is the output of different memory unit, yt−3,yt−2,yt−1 is the prediction value of each layer, and so on; we can get the final future value yt.The comparison of the prediction results of these five models is shown in Tables3 and 4.Table 3
Comparison of ration consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Year
2014
2015
2016
2017
2018
Original ration consumption
25532
24411
24192
23784
23361
MLR model
Predicted value
23081
23437
23436
23314
23243
APE
9.60%
3.99%
3.13%
1.98%
0.51%
MAPE
3.84%
Theil’sU
0.0262
ARIMA model
Predicted value
25410
25647
25678
25831
25831
APE
0.48%
5.06%
6.14%
8.61%
10.57%
MAPE
6.17%
Theil’sU
0.0336
BP model
Predicted value
25760
28837
26097
28026
26489
APE
0.89%
18.13%
7.87%
17.84%
13.39%
MAPE
11.62%
Theil’sU
0.0623
LSTM model
Predicted value
22424
22861
24273
22158
22928
APE
12.17%
6.35%
0.33%
6.84%
1.85%
MAPE
5.51%
Theil’sU
0.0366
Modified MLR model
Predicted value
24029
24031
24205
24067
24059
APE
5.89%
1.56%
0.05%
1.19%
2.99%
MAPE
2.34%
Theil’sU
0.0159Table 4
Comparison of feed grain consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Year
2014
2015
2016
2017
2018
Original feed grain consumption
22247
23047
23317
23920
25300
MLR model
Predicted value
23921
24155
24286
24393
24483
APE
7.52%
4.81%
4.16%
1.98%
3.23%
MAPE
4.34%
Theil’sU
0.0226
ARIMA model
Predicted value
21619
21703
21730
21730
21730
APE
2.82%
5.83%
6.81%
9.16%
14.11%
MAPE
7.75%
Theil’sU
0.0466
BP model
Predicted value
23145
23233
26432
25253
25977
APE
4.04%
0.81%
13.36%
5.57%
2.68%
MAPE
5.29%
Theil’sU
0.0330
LSTM model
Predicted value
24378
24876
24956
27131
27481
APE
9.58%
7.94%
7.03%
13.42%
8.62%
MAPE
9.32%
Theil’sU
0.0459
Modified MLR model
Predicted value
23056
23100
23125
23146
23163
APE
3.64%
0.23%
0.82%
3.24%
8.45%
MAPE
3.27%
Theil’sU
0.0232From Tables3 and 4, for the ARIMA model, the closer the forecast year is, the higher the prediction accuracy of the model is. For the MLR model, the further away the predicted year is, the higher the prediction accuracy of the model is. This is because the time series prediction model is more suitable for short-term prediction and the multiple regression model is suitable for long-term prediction. However, it can be found that the prediction accuracy of the proposed modified MLR model is higher than that of the single time series model and the classical MLR model. It is mainly because the combined model can not only track the change of impact factors but also grasp the internal change trend of grain consumption, thereby leading to a small error range. And the proposed model is superior to the BP model and LSTM model for the two indexes of MAPE and Theil’s U.Meanwhile, in order to verify the stability of the modified MLR model, two experiments are carried out: (1) the training time interval is 1981–2014, and the predicting time interval is 2015–2018; (2) the training time interval is 1981–2013, and the predicting time interval is 2014–2017. The prediction performance of ration consumption and feed grain consumption is shown in Table5.Table 5
Comparison of feed grain consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Training time interval (predicting time interval)
1981–2014 (2015–2018)
1981–2013 (2014–2017)
Year
2015
2016
2017
2018
2014
2015
2016
2017
Ration consumption
Actual value
24411
24192
23784
23361
25532
24411
24192
23784
Predicted value
24142
23977
23948
23985
24029
24031
24205
24067
APE
1.10%
0.89%
0.69%
2.67%
5.89%
1.56%
0.05%
1.19%
RMSE
365.67
787.79
MAPE
1.34%
2.17%
Theil’sU
0.0076
0.0162
Feed grain consumption
Actual value
23047
23317
23920
25300
22247
23047
23317
23920
Predicted value
22982
22994
23006
22999
23056
23100
23125
23146
APE
0.28%
1.39%
3.82%
9.09%
3.64%
0.23%
0.82%
3.24%
RMSE
1248.85
568.38
MAPE
3.64%
1.98%
Theil’sU
0.0266
0.0123From Table5, it can be seen that the prediction error limit is kept within 2% for ration and feed grain consumption under different time intervals, which indicates that the proposed model can realize stable prediction with satisfying accuracy for the medium-short-term prediction.


## 3.3.1. Fitting Results

By the proposed method as equation (24), for every kind of grain consumption, two key impact factors and the historical grain consumption are chosen as the input of the time series-MLR joint model. However, for the classical MLR model, the input is only the two chosen key impact factors. The fitting results of the two models are shown in Figure 5(a) for feed grain prediction and Figure 5(b) for ration grain prediction.Figure 5
The fitting results of different models (1981–2013). (a) Feed grain consumption. (b) Ration consumption.
(a)
(b)From Table2, it can be seen that the fitting goodness of ration and feed grain all can be improved. The fitting results of the proposed model can provide better prediction reliability.Table 2
The fitting goodness of grain consumption trend.
Grain consumption
Ration
Feed grain
Predicted model
Multivariate regression
Joint model
Multivariate regression
Joint model
R
2
0.9656
0.9863
0.9917
0.9965


## 3.3.2. Analysis of Prediction Results

For the experiments of predicting the future trend of grain consumption, we choose the filter window width equal to 3, 5, or 7 and establish the prediction model. Different filtering window widths will affect the model prediction performance; therefore, if we feed the prediction error back to adjust the filter parameters, the optimal filtering window width can be chosen. Here, the chosen optimal window width is 5 and 3, respectively, for ration consumption prediction and feed grain consumption prediction. In this paper, the MLR model and modified MLR model are established based on these optimal preprocessing methods.In addition, in our simulation, a classical time series ARIMA model is adopted to predict the ration and feed grain consumption in China. Besides, deep learning prediction is also considered for comparison, which is the emerging research hotspot in recent years [24–30]. Five kinds of prediction models are considered for comparing APE, MAPE, and Theil’s U performance.(1)
The classical time series ARIMA model:(25)yt=φ1yt−1+φ2yt−2+⋯+φpyt−p+εt+θ1εt−1+⋯+θqεt−q,whereyt−1,yt−2,…,yt−p is the actual value of the past p years; the order of p and q can be realized by model recognition and φ1,φ2,…,φp and θ1,θ2,…,θq can be calculated by the least-square method.(2)
The MLR model:The two impact factors with the greatest relational degree as the key impact factors in the modeling were input into the MLR prediction model:(26)yt=α0+α1x1t+α2x2t,wherex1t,x2t are key impact factors; model parameters α0,α1,α2 can be calculated by the data barycenter method.(3)
The modified MLR model:Considering the complexity of the algorithm and the prediction performance, the proposed MLR model is(27)yt=β0+β1x1t+β2x2t+ϕ1yt−1,where model parametersβ0,β1,β2, and ϕ1 can be calculated by the data barycenter method; yt is the future prediction value; yt−1 is the historical data of the prediction variable; and x1t,x2t are the chosen optimal impact factors.(4)
The BP model:Backpropagation (BP) is one of the artificial neural networks (ANNs) applied widely, which can train the model parameters by reversing the error propagation algorithm [24–28]. The prediction frame of the BP model is shown in Figure 6.There are three layers to predict the future value from the three inputted history data; here the hidden layer is used to analyze the characteristics of input data. In Figure6, xt−2,xt−1,xt is the inputted history data of grain consumption, yt is the future prediction value, and W1 and W2 are the weighting values of the trained model.(5)
The LSTM model:LSTM is one special recurrent neural network (RNN), which substitutes the original neural units of hidden layers for memory units [29, 30]. Each memory unit has an input door, output door, forgotten door, and memory cells, as shown in Figure 7.Figure 6
The prediction frame of the BP model.Figure 7
The prediction frame of the LSTM model.Herext−2,xt−1,xt is the input history data of grain consumption, ct−2,ct−2,ct−1,ct is the output of different memory unit, yt−3,yt−2,yt−1 is the prediction value of each layer, and so on; we can get the final future value yt.The comparison of the prediction results of these five models is shown in Tables3 and 4.Table 3
Comparison of ration consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Year
2014
2015
2016
2017
2018
Original ration consumption
25532
24411
24192
23784
23361
MLR model
Predicted value
23081
23437
23436
23314
23243
APE
9.60%
3.99%
3.13%
1.98%
0.51%
MAPE
3.84%
Theil’sU
0.0262
ARIMA model
Predicted value
25410
25647
25678
25831
25831
APE
0.48%
5.06%
6.14%
8.61%
10.57%
MAPE
6.17%
Theil’sU
0.0336
BP model
Predicted value
25760
28837
26097
28026
26489
APE
0.89%
18.13%
7.87%
17.84%
13.39%
MAPE
11.62%
Theil’sU
0.0623
LSTM model
Predicted value
22424
22861
24273
22158
22928
APE
12.17%
6.35%
0.33%
6.84%
1.85%
MAPE
5.51%
Theil’sU
0.0366
Modified MLR model
Predicted value
24029
24031
24205
24067
24059
APE
5.89%
1.56%
0.05%
1.19%
2.99%
MAPE
2.34%
Theil’sU
0.0159Table 4
Comparison of feed grain consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Year
2014
2015
2016
2017
2018
Original feed grain consumption
22247
23047
23317
23920
25300
MLR model
Predicted value
23921
24155
24286
24393
24483
APE
7.52%
4.81%
4.16%
1.98%
3.23%
MAPE
4.34%
Theil’sU
0.0226
ARIMA model
Predicted value
21619
21703
21730
21730
21730
APE
2.82%
5.83%
6.81%
9.16%
14.11%
MAPE
7.75%
Theil’sU
0.0466
BP model
Predicted value
23145
23233
26432
25253
25977
APE
4.04%
0.81%
13.36%
5.57%
2.68%
MAPE
5.29%
Theil’sU
0.0330
LSTM model
Predicted value
24378
24876
24956
27131
27481
APE
9.58%
7.94%
7.03%
13.42%
8.62%
MAPE
9.32%
Theil’sU
0.0459
Modified MLR model
Predicted value
23056
23100
23125
23146
23163
APE
3.64%
0.23%
0.82%
3.24%
8.45%
MAPE
3.27%
Theil’sU
0.0232From Tables3 and 4, for the ARIMA model, the closer the forecast year is, the higher the prediction accuracy of the model is. For the MLR model, the further away the predicted year is, the higher the prediction accuracy of the model is. This is because the time series prediction model is more suitable for short-term prediction and the multiple regression model is suitable for long-term prediction. However, it can be found that the prediction accuracy of the proposed modified MLR model is higher than that of the single time series model and the classical MLR model. It is mainly because the combined model can not only track the change of impact factors but also grasp the internal change trend of grain consumption, thereby leading to a small error range. And the proposed model is superior to the BP model and LSTM model for the two indexes of MAPE and Theil’s U.Meanwhile, in order to verify the stability of the modified MLR model, two experiments are carried out: (1) the training time interval is 1981–2014, and the predicting time interval is 2015–2018; (2) the training time interval is 1981–2013, and the predicting time interval is 2014–2017. The prediction performance of ration consumption and feed grain consumption is shown in Table5.Table 5
Comparison of feed grain consumption prediction results in China (2014–2018) (unit: ten thousand tons).
Training time interval (predicting time interval)
1981–2014 (2015–2018)
1981–2013 (2014–2017)
Year
2015
2016
2017
2018
2014
2015
2016
2017
Ration consumption
Actual value
24411
24192
23784
23361
25532
24411
24192
23784
Predicted value
24142
23977
23948
23985
24029
24031
24205
24067
APE
1.10%
0.89%
0.69%
2.67%
5.89%
1.56%
0.05%
1.19%
RMSE
365.67
787.79
MAPE
1.34%
2.17%
Theil’sU
0.0076
0.0162
Feed grain consumption
Actual value
23047
23317
23920
25300
22247
23047
23317
23920
Predicted value
22982
22994
23006
22999
23056
23100
23125
23146
APE
0.28%
1.39%
3.82%
9.09%
3.64%
0.23%
0.82%
3.24%
RMSE
1248.85
568.38
MAPE
3.64%
1.98%
Theil’sU
0.0266
0.0123From Table5, it can be seen that the prediction error limit is kept within 2% for ration and feed grain consumption under different time intervals, which indicates that the proposed model can realize stable prediction with satisfying accuracy for the medium-short-term prediction.


## 4. Conclusions

In this paper, one new modified MLR model combined with time series forecasting theory is proposed. In the proposed model, the historical data of prediction variable will be considered as one key impact factor, which will be imputed to the MLR model with the other chosen impact factors. Thus, it can make full use of the advantages of the time series model and the MLR model, and the prediction performance has been improved strikingly, compared with the MLR model, ARIMA model, and even the BP model and LSTM model. Besides, there are other innovative works; the original data are preprocessed by moving average filtering under optimal window, and the data barycenter method is also adapted to compute the model parameters of MLR for improving the robustness. The proposed prediction model can be applied widely in short-medium-term prediction of grain prediction or the other fields.



---

*Source: 3956717-2020-05-13.xml*