Optimal solution for novel grey polynomial prediction model

Accepted Manuscript

Optimal solution for novel grey polynomial prediction model Baolei Wei, Naiming Xie, Aqin Hu PII: DOI: Reference:

S0307-904X(18)30288-9 10.1016/j.apm.2018.06.035 APM 12334

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

1 March 2018 12 June 2018 18 June 2018

Please cite this article as: Baolei Wei, Naiming Xie, Aqin Hu, Optimal solution for novel grey polynomial prediction model, Applied Mathematical Modelling (2018), doi: 10.1016/j.apm.2018.06.035

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • The novel grey polynomial prediction model with nonlinear parameters is proposed. • The parameter optimization criteria in the existing grey prediction models are unified by weighted least square method.

• The algorithmic framework for selecting the polynomial order and estimating the nonlinear parameters

CR IP T

is proposed.

• The affine transformation of the accumulating sequence proves to be independent of modeling perfor-

AC

CE

PT

ED

M

AN US

mance.

1

ACCEPTED MANUSCRIPT

Optimal solution for novel grey polynomial prediction model Baolei Weia,∗, Naiming Xiea,b , Aqin Hua a College

of Economics and Management, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China of Grey System Studies, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China

CR IP T

b Institute

Abstract

The grey prediction model, as a time-series analysis tool, has been used in various fields only with partly known distribution information. The grey polynomial model is a novel method to solve the problem that the original sequence is in accord with a more general trend rather than the special homogeneous or nonhomogeneous trend, but how to select the polynomial order still needs further study. In this paper the tuned background coefficient is introduced into the grey polynomial model and then the algorithmic frame-

AN US

work for polynomial order selection, background coefficient search and parameter estimation is proposed. The quantitative relations between the affine transformation of accumulating sequence and the parameter estimates are deduced. The modeling performance proves to be independent of the affine transformation. The numerical example and application are carried out to assess the modeling efficiency in comparison with other conventional models.

1. Introduction

ED

1

M

Keywords: grey prediction model, background coefficient, weighted least square method, line search

Time series prediction has always been an important issue in economic, finance, marketing, as well as

3

social problems. Nowadays, there exist hundreds of sequential prediction techniques, such as LR (Linear

4

Regression), ES (Exponential Smoothing), ARIMA (Auto Regressive Integrated Moving Average), SVR

5

(Support Vector Regression), HMM (Hidden Markov Model) and ANN (Artificial Neural Network). The

6

results of these methods are always ensured just under a basic assumption of known distribution or large

7

sample. However, it is sometimes impossible to collect time series with large sample size especially for a

8

new system, such as the cold-start problem in recommendation system. In Deng’s seminal study [1], grey

9

system theory focusing on the small sample problems was first proposed. Grey prediction is an important

10

branch of grey system theory and generally written as GM(p, q) which is a p-th order model with q variables.

12 13

CE

AC

11

PT

2

Over the last 35 years, it has been widely utilized in various fields such as energy, economics, environment, transportation and so on [2–4]. The GM(1,1) is the foundation of grey prediction models and has been substantiated in real world

14

applications. In the meantime, much importance has been consistently placed on improving the modeling

15

accuracy, and most of the effort is invested in the following aspects: background coefficient optimization [5–8],

16

initial condition selection [9–12], hybrid approach combined with background coefficient and initial condition ∗ Corresponding

author. Email addresses: [email protected] (Baolei Wei), [email protected] (Naiming Xie)

Preprint submitted to Applied Mathematical Modelling

June 22, 2018

ACCEPTED MANUSCRIPT

17

optimization [13], fractional order accumulation [14], time-varying weighted accumulating generation [15],

18

properties analysis including the necessary and sufficient condition for modeling [16], the relative error upper

19

bound estimation [17], the independence of model’s accuracy and affine transformation [18], the influence of

20

sample size on modeling performance [19–21]. The modeling procedures of univariate grey prediction models are represented by GM(1,1) model, and

22

many extension models have already been reported in recent years. Guo and Guo [22] probed into the nature

23

of grey prediction model and looked ahead to the extension principle but short of details. Subsequently,

24

Cui et al. [23] put forward NGM(1,1) model by replacing the grey action constant with time linear function

25

so that it is applicable to the sequence in accord with quasi non-homogeneous exponential law. Similarly,

26

Qian, Dang and Liu [24] developed GM(1,1,tα ) model by adding a time power function to the grey action

CR IP T

21

term, and NGM(1,1) model proves to be a special form of this model with α = 1. In a general way, Ma, Hu and Liu [25] constructed a kernel regularized grey model KGM(1,1) by using a function of time as the

29

grey action term, and Luo and Wei [26] proposed a grey polynomial model GPM(1,1,N ) by using a time

30

polynomial function as the grey action term. These two models have larger application ranges, but each of

31

them suffers from some limitations, respectively, such as the selection of regularized parameter and kernel

32

function in KGM(1,1) model, and the selection of polynomial order in GPM(1,1,N ) model.

AN US

28

27

The discrete grey prediction models have long been considered novel methodologies, and they are essen-

34

tially the extensions of the existing models, such as DGM(1,1) model [27] corresponding to the GM(1,1)

35

model, and NDGM(1,1) model [28] corresponding to the NGM(1,1) model, etc. The multivariate grey predic-

36

tion models can be also regarded as extensions, such as GM(1,N ) model [29] corresponding to the GM(1,1)

37

model, NGM(1,N ) model [30] corresponding to the NGM(1,1) model, and KGM(1,N ) model corresponding

38

to the KGM(1,1) model [31], etc. In addition, the optimization strategies in the GM(1,1) model have been

39

gradually employed to improve the accuracies of the discrete and multivariate grey models. For example,

40

Zeng and Li [32] improved the NGM(1,N ) model by using the particle swarm optimization algorithm to

41

search the optimal background coefficient.

ED

M

33

In summary, the univariate model is the basis of the discrete and multivariate model construction. Luo

43

and Wei [26] proved that the above GM(1,1), NGM(1,1) and GM(1,1,tα ) model are all special forms of the

44

GPM(1,1,N ) model with deterministic polynomial order, and this model is applicable to including but not

45

limited to the quasi homogeneous and non-homogeneous exponential sequences, or even some fluctuating

46

sequences. However, there are still some drawbacks that are worth to be further investigated. In this paper,

47

the major points are summarized as follows: (1) the polynomial order is usually unknown in practice, and

48

the empirical criterion [26] is given under the assumption that the original sequence is composed of the

49

exponential and polynomial functions without noise and thus is likely to be inaccurate; (2) the background

51 52 53 54 55

56 57

CE

AC

50

PT

42

coefficient is fixed at 0.5 in the existing model, which may result in poor predictions sometimes; (3) the least squared error, least squared percentage error and least absolute percentage error are three common criteria in model parameter estimation, nevertheless the relationships among them and their corresponding properties have not yet been researched. In order to address these drawbacks in the existing approaches, we focus on the polynomial order selection and parameter estimation. In detail, the main work lies in the following aspects: • The grey polynomial model with tuned background coefficient is derived from the fundamental form,

and the algorithmic framework for nonlinear parameter (polynomial order, background coefficient and 3

ACCEPTED MANUSCRIPT

58

59

model parameters) estimation is presented. • The background coefficient and model parameters are estimated by line-search based weighted least

60

square method in which the weighted objective function is turned out to be the unified representation

61

of the three above optimization criteria in grey prediction models.

63

• The properties including the comparison between the proposed model and other grey prediction models as well as the influence of affine transformation on modeling performance are deeply analyzed.

CR IP T

62

64

The remainder of this paper has the following structure: the proposed model and its parameter estimation

65

algorithm are elaborated in Section 2; the properties are researched in Section 3; the numerical example and

66

application are respectively presented in Sections 4 and 5; the conclusions are drawn in Section 6.

67

2. Construction of GPMB(1,1,N ) model

68

2.1. Representation of GPMB(1,1,N ) model

70

71

where

k X i=1

(1)

x(0) (k), k = 1, 2, · · · , n.

(2)

M

The grey polynomial model is suggested to be the following difference equation x(1) (k − 1) + x(1) (k) 2k − 1 k N +1 − (k − 1) + β0 + β1 + · · · + 2 2 N +1

N +1

βN ,

ED

x(0) (k) = α 73

x(0) (1), x(0) (2), · · · , x(0) (n) , n ≥ 4, its first order

accumulating generation sequence is defined as n o X (1) = x(1) (1), x(1) (2), · · · , x(1) (n) , x(1) (k) =

72



Assuming that the original sequence is X (0) =

AN US

69

and its whitening equation is expressed as

d (1) x (t) = αx(1) (t) + β0 + β1 t + · · · + βN tN , dt

(4)

PT

74

(3)

where the order of polynomial satisfies N ≤ n − 4 and α, β0 , β1 , · · · , βN are model parameters.

The difference expression Eq. (3) is a roughly discrete approximation of Eq. (4), and the accuracy of

76

this approximate formula can be further improved. Without loss of generality, consider the integration on the both side of Eq. (4) in the interval [k − 1, k], then Z k Z k Z k Z (1) (1) x (t)dt = α x (t)dt + β0 dt + β1

AC

77

CE

75

k−1

78 79

k−1

k−1

k

k−1

tdt + · · · + βN

Z

k

tN dt.

(5)

k−1

According to the Newton-Leibniz formula and the definition of accumulating generation in Eq. (2), the

left side of Eq. (5) can be expressed as Z k x(1) (t)dt = x(1) (k) − x(1) (k − 1) = x(0) (k)

(6)

k−1

80

and the right side of Eq. (5) is equal to Z k N +1 k N +1 − (k − 1) 2k − 1 α x(1) (t)dt + β0 + β1 + · · · + βN . 2 N +1 k−1 4

(7)

ACCEPTED MANUSCRIPT

81 82

Rk

x(1) (t)dt denotes the area between t-axis and the curve x(1) (t) in the interval [k − 1, k]. It was calculated using the trapezoid formula in the traditional GPM(1,1,N ) model It is notable that the integration term

k−1

83

[26] and many other recent researches [14–17, 19–25, 29–31]. But the trapezoid formula is not accurate

84

enough regardless of that x(1) (t) is convex, concave or non-convex in each subinterval as shown in Figure 1.

85

y

y  x (1) (t )

x (1 ) (k )

x (1 ) (k )

x (1) (k  1)

x (1) (k  1) k 1 (a) convex

k

y

y  x (1) (t )

x (1) (k  1)

k 1 (b) concave

t

k

t

Figure 1: Geometric schematic representation of the integration term

87

Rk

k−1

theorem for definite integrals, there always exists a real number λ ∈ [0, 1] such that Z k x(1) (t)dt = λx(1) (k − 1) + (1 − λ)x(1) (k).

(8)

Substituting Eqs. (6), (7) and (8) into Eq. (5) gives the more accurate approximate formula, termed as GPMB(1,1,N ) model, expressed as

N +1 h i 2k − 1 k N +1 − (k − 1) x(0) (k) = α λx(1) (k − 1) + (1 − λ)x(1) (k) + β0 + β1 + · · · + βN , 2 N +1

where the unknown parameter λ ∈ [0, 1] is called background coefficient.

91

2.2. Framework for polynomial order selection and parameter estimation

(9)

ED

90

According to Eq. (9), it can be seen that the model parameters α, β0 , · · · , βN , background coefficient

PT

92

t

x(1) (t)dt.

M

89

k

(c) nonconvex

Assume that x(1) (t) is a monotone function in each subinterval, then according to the first mean value

k−1

88

k 1

AN US

86

y  x (1) (t )

x (1 ) (k )

CR IP T

y

93

λ and polynomial order N have nonlinear relationship, which makes the ordinary least square method

94

inapplicable. In practice, the alternative polynomial order is always suggested to be N ∈ N = {0, 1, 2, 3} to

avoid overfitting and ill-condition problem [26, 33, 34], and therefore the enumeration strategy is acceptable.

96

On the basis of λ ∈ [0, 1], the line-search based weighted least square method is employed to search the

98 99

100 101 102 103

background coefficient and estimate the model parameters.

Line search algorithm is a numerical optimization method that can be applied to discontinuous or even

AC

97

CE

95

indifferentiable functions [34]. Each iteration is given by λi+1 = λi + ηpi

(10)

where η > 0 is the step length and pj = ±1 is the descent direction to guarantee the reduction of objective function.

Under each background coefficient, the model parameters α, β0 , β1 , · · · , βN are estimated from Eq. (9)

using the weighted least square method n o  T
−1 T 2 2 κ ˆ= α D W Y, ˆ βˆ0 βˆ1 · · · βˆN = arg min φ(κ) = kW (Y − Dκ)k2 = D T W 2 D κ

5

(11)

ACCEPTED MANUSCRIPT

105

106 107

where W = diag (w2 , w3 , · · · , wn ) is a positive definite matrix in which wk > 0 (k = 2, 3, · · · , n) is the weight of the k-th sample, and  λi x(1) (1) + (1 − λi )x(1) (2) 1  (1) (1) 1  λi x (2) + (1 − λi )x (3) D= .. ..   . . (1) (1) λi x (n − 1) + (1 − λi )x (n) 1

3 2 5 2

.. .

2n−1 2

··· ··· .. . ···



 (0)  x (2)   x(0) (3)    , Y =   ..  . ..   .   . nN +1 −(n−1)N +1 x(0) (n) 2N +1 −1 N +1 3N +1 −2N +1 N +1

Substituting the estimates α ˆ , βˆ0 , βˆ1 , · · · , βˆN into the general solution to the whitening equation (4), the

time response function (the details can be seen in reference [26]) is obtained as ˆ x ˆ(1) (t) = ceαt + γ0 + γ1 t + · · · + γN t N ,

110

(14)

constant c is optimized by minimizing the sum of the squared errors,

c

n  X

k=1

x

(1)

(k) − x ˆ

(1)

(k)

2

n N X e2αˆ − 1 X (1) = 2αn x (k) − γi k i e ˆ −1 i=0 k=1

!

ˆ eα(k−2) .

(15)

Substituting the polynomial coefficients in Eq. (14) and the estimated integration constant in Eq. (15) into the time response function Eq. (13), the time response sequence is obtained as ! # " N n X e2αˆ − 1 X i α(k−2) ˆ ˆ (1) (1) γi k e eαk + γ0 + γ1 k + · · · + γN k N , x (k) − x ˆ (k) = 2αn e ˆ −1 i=0

ED

112

−1   0 βˆ0 ..   ˆ  β1  .    . .    ..  N βˆN −ˆ α

(13)

Different form the traditional method taking x ˆ(1) (1) = x(1) (1) as the initial condition, the integration

cˆ = arg min 111

··· .. . .. . ···

AN US

109

where the polynomial coefficients γ0 , γ1 , · · · , γN are   −ˆ α 1 γ0   γ1   0 −ˆ α   γ= . = .. ..  ..    . . γN 0 0

M

108

(12)

N +1

CR IP T

104

(16)

k=1

116

PT

114 115

where k = 1, 2, · · · , n + p and p ∈ N+ is the number of samples to be predicted.

By using the first order inverse accumulating generation operator, the fitted and predicted values corre-

sponding to the original sequence can be obtained as n o ˆ (0) = x X ˆ(0) (1), x ˆ(0) (2), · · · , x ˆ(0) (n), x ˆ(0) (n + 1), · · · , x ˆ(0) (n + p) , where

CE

113

117 118

AC

x ˆ

(0)

(k) =

(

x ˆ(1) (1), x ˆ(1) (k) − x ˆ(1) (k − 1),

k=1 . k = 2, 3, · · · , n + p

(17)

(18)

The mean absolute percentage error (MAPE) is used to evaluate the fitting and predicting performance,

which is defined as MAPE =

T2 (0) X x (k) − x 1 ˆ(0) (k) × 100%, T2 − T1 + 1 x(0) (k)

(19)

k=T1

119

where MAPE is the fitting performance measure with T1 = 1 and T2 = n and the predicting performance

120

measure with T1 = n + 1 and T2 = n + p.

121 122

To sum up, the algorithmic framework for polynomial order selection, background coefficient search and parameter estimation is presented in Algorithm 1. 6

ACCEPTED MANUSCRIPT

Algorithm 1 Input: Original sequence X (0) , weighted matrix W . 1: Divide the original sequence into two parts: n o n o (0) (0) Xtrain = x(0) (1), x(0) (2), · · · , x(0) (n1 ) and Xtest = x(0) (n1 + 1), x(0) (n1 + 2), · · · , x(0) (n) ,

124

2.3. Computational steps of GPMB(1,1,N ) model

According to the above modeling processes, the computational steps of GPMB(1,1,N ) model can be

ED

123

M

AN US

CR IP T

where n denotes the number of samples in X (0) and traditionally n1 is the largest integer less than or equal to 54 n. 2: for all N ∈ N = {0, 1, 2, 3} do 3: Initialize the count i = 0 and the background coefficient λi = 0. 4: while λi ≤ 1 do 5: Calculate the estimates of model parameters by Eq. (11), the polynomial coefficients by Eq. (14) and the initial condition by Eq. (15); 6: Calculate the time response sequence by Eq. (16) and the fitted and predicted values corresponding to X (0) by Eq. (17); 7: Update the background coefficient by Eq. (10) in which the step length and descent direction can be set as η = 0.01 and pi = 1. 8: end while 9: Pick the best background coefficient and denote it as λopt . Here best is defined as having the smallest (0) MAPE on Xtrain , abbreviated as MAPEtrain . 10: end for 11: Select the best polynomial order and denote it as Nopt . Here best is defined as having the smallest (0) MAPE on Xtest , abbreviated as MAPEtest . Output: Nopt , λopt .

125

summarized as follows.

126

Step 1 Given the optimization criterion (i.e., the weighted matrix W ), compute the best polynomial order

130 131

132

133 134 135

PT

129

Step 2 Based on the whole sequence X (0) including all samples, calculate the time response sequence by substituting Nopt and λopt into Eq. (16); Step 3 Compute the fitted and predicted values corresponding to the original sequence by Eqs. (17) and

CE

128

and background coefficient according to Algorithm 1;

(18), as well as the performance measures by Eq. (19). 3. Affine property of GPMB(1,1,N ) model

AC

127

The weighted least square method is a unified representation of the common optimization criteria includ-

ing but not limited to the least squared error, least squared percentage error and least absolute percentage error. When wk = 1, the model parameters in Eq. (11) are estimated under the least squared error criterion.

136

ˆ = arg min κ κ

(

φ(κ) = kW (Y −

2 Dκ)k2

=

n  X

k=2

7

x

(0)

(k) − x ˆ

(0)

2 (k)

)

.

(20)

ACCEPTED MANUSCRIPT

137

When wk = 1/x(0) (k), the model parameters in Eq. (11) are estimated under the least squared percentage

138

error criterion. ˆ = arg min κ κ

140

φ(κ) = kW (Y −

2 Dκ)k2

=

2 n  (0) X x (k) − x ˆ(0) (k) x(0) (k)

k=2

)

.

(21)

  When wk = 1/ x(0) (k) x(0) (k) − x ˆ(0) (k) , the model parameters in Eq. (11) are estimated under the least

absolute percentage error criterion. ( ˆ = arg min κ κ

φ(κ) = kW (Y −

2 Dκ)k2

) n (0) (0) X x (k) − x ˆ (k) . = x(0) (k)

CR IP T

139

(

k=2

(22)

Due to that the sample weights is not independent of the fitted values, the objective function in Eq. (22) is non-differentiable and the analytic solution in Eq. (11) is not applicable any more. The objective function

143

in Eq. (21) is the approximation of that in Eq. (22) which is always solved by the intelligent algorithms

144

such as particle swarm optimization and genetic algorithm [4].

145

Remark 1. Under the least squared error criterion, i.e., the sample weights wk = 1, let the background

146

coefficient be 0.5 and el be the (n − 1) × 1 vector having 1 in the l-th element and 0’s in others, then in

150 151

M

149

(1) When N = 0, GPMB(1,1,N ) model yields the GM(1,1) model [4] defined as i α h (1) x(0) (k) = x (k − 1) + x(1) (k) + β0 , 2 and the estimates of model parameters are

 where D1 = D e1

 ˆ= α κ ˆ

ED

148

GPMB(1,1,N ) model

 e2 .

x(0) (k) =

153 154

155

156

T

= D1T D1


−1

D1T Y

i 2k − 1 α h (1) x (k − 1) + x(1) (k) + β0 + β1 , 2 2

CE

and the estimates of model parameters are

 where D2 = D e1

AC

152

βˆ0

(2) When N = 1, GPMB(1,1,N ) model yields the NGM(1,1) model [23] defined as

PT

147

AN US

141 142

e2

 e3 .

 ˆ= α κ ˆ

βˆ0

βˆ1

T

= D2T D2


−1

D2T Y

(3) When β1 = β2 = · · · = βN −1 = 0, GPMB(1,1,N ) model yields the GM(1,1,tα ) model [24] defined as x(0) (k) =

i α h (1) k N +1 − (k − 1)N +1 x (k − 1) + x(1) (k) + β0 + βN , 2 N +1

and the estimates of model parameters are

 where D3 = D e1

e2

 ˆ= α κ ˆ

 eN +2 .

βˆ0

βˆN

T

8

= D3T D3


−1

D3T Y

ACCEPTED MANUSCRIPT

157

Remark 1 indicates that the existing GM(1,1), NGM(1,1) and GM(1,1,tα ) model are all special cases

158

of GPMB(1,1,N ) model with fixed background coefficient. It is convenient to obtain GM(1,1), NGM(1,1),

159

160

GM(1,1,tα ) and other novel models by combining polynomial order and optimization criterion. n o (1) (1) (1) (1) (1) Theorem 1. Assume XAF = xAF (1), xAF (2), · · · , xAF (n) , where xAF (k) = ρx(1) (k) + ξ, ρ 6= 0, to be

161

the affine transformation of the accumulating generation sequence X (1) , then the fitted and predicted values

162

of GPMB(1,1,N ) models respectively constructed by X (1) and XAF satisfy x ˆAF (k) = ρˆ x(1) (k) + ξ, k =

163

1, 2, · · · , n + p.

Proof. From the definitions of accumulating generation and affine transformation, it follows that  (0)      (1) (1) xAF (2) − xAF (1) xAF (2) ρx(0) (2)   (0)   (1) (1) (0)   xAF (3)   xAF (3) − xAF (2)    ρx (3)     = = YAF =   = ρY .. ..    ..    .    .   . (0) (1) (1) (0) ρx (n) xAF (n) − xAF (n − 1) xAF (n) and

DAF



  ρ λx(1) (1) + (1 − λ)x(1) (2) + ξ     ρ λx(1) (2) + (1 − λ)x(1) (3) + ξ  = ..  .  (1)  ρ λx (n − 1) + (1 − λ)x(1) (n) + ξ

1 1 .. . 1

3 2 5 2

.. .

2n−1 2

··· ···

2N +1 −1 N +1 3N +1 −2N +1 N +1

··· ···

.. .

nN +1 −(n−1)N +1 N +1

where P and Q are two (N + 2)-order nonsingular matrices respectively defined as     ρ 0 ··· 0 1 0 ··· 0 0 1 · · · 0  ξ 1 · · · 0     . P = . . . and Q =  . . .  . . . ...  . . ..   .. ..  .. ..  0 0 ··· 1 0 0 ··· 1



   = DP Q  

(23)

(24)

169 170 171

172 173

PT

 0 0  ..  . .

(25)

1

Under the same optimization criterion, the weighted matrices always satisfy WAF = $W , where $ is a

CE

168

Because the coefficient ρ 6= 0, P is invertible and the inverses of P and Q are    1 0 ··· 0 1 0 ···    −ξ 1 ··· 0 ρ · · · 0 1   and Q−1 =  . P −1 =  . . .  .. . . . . . ..  ρ  .. ..  .. . . 0 0 ··· 0 0 ··· ρ constant only depending on ρ, for instance, $ equals

1 ρ

under the three above optimization criteria. Denoting  0  (1) 0 T , it ˆ AF = α the model parameters of GPMB(1,1,N ) model constructed by XAF as κ ˆ βˆ00 βˆ10 · · · βˆN follows that h
−1 T
−1 T 2 i T 2 2 ˆ AF = DAF ˆ . (26) κ WAF DAF DAF WAF YAF = ρQ−1 P −1 D T W 2 D D W Y = ρQ−1 P −1 κ

AC

167

ED

M

166

CR IP T

165

(1)

AN US

164

(1)

Substituting Eq. (25) into Eq. (26) gives the model parameters satisfying α ˆ0 = α ˆ , βˆ00 = ρβˆ0 − ξ α ˆ , βˆj0 =

ρβˆj , j = 2, · · · , N , and the polynomial coefficients satisfying −1   0  −ˆ    α 1 ··· 0 −ξ α ˆ + ρβˆ0 γ0 ργ0 + ξ   .  γ10      ργ1  ρβˆ1 0 −ˆ α .. 0         γAF =  .  =    =  .. , .   . .. .. ..  ..      .  ..  . N . 0 ˆ γN ργ N ρβN 0 0 · · · −ˆ α 9

(27)

ACCEPTED MANUSCRIPT

174

as well as the integration constant satisfying # ! " n N X e2αˆ − 1 X (1) α(k−2) ˆ i x (k) − = ρˆ c. cˆAF = ρ 2αn γi k e e ˆ −1 i=0

(28)

k=1

175

Substituting Eqs. (26), (27) and (28) into the time response sequence gives that " ! # n N N X X e2αˆ − 1 X (1) (1) i α(k−2) ˆ αk ˆ x ˆAF (k) = ρ 2αn x (k) − γ k e e + ρ γi k i + ξ = ρˆ x(1) (k) + ξ. i e ˆ −1 i=0 i=0

176

By the inverse affine transformation,

1 ρ



(1)

x ˆAF − ξ



(29)

CR IP T

k=1

= x ˆ(1) (k), that is to say, the fitted and predicted

177

values are invariant whether the affine transformation is applied to the accumulating generation sequence

178

or not.

179

Furthermore, if the translation coefficient ξ = 0, the affine transformation of the accumulating generation sequence is equivalent to the multiple transformation of the original sequence. From Theorem 1, the multiple transformation of the original sequence has no influence on modeling performance. Therefore, a suitable ρ

182

can be selected to reduce the condition number in Eq. (11) without change of accuracy [35].

183

4. Numerical example

AN US

180 181

In this section, the GPMB(1,1,N ) model is build on the original sequence shown in Table 1, and the

185

conventional GPM(1,1,N ) model [26] is also constructed to illustrate the improvements. The modeling

186

results are also compared with the commonly used prediction models: NDGM(1,1) [28], Verhulst [36] and

187

CPR (Cubic Polynomial Regression) [33].

ED

M

184

Table 1: The actual values of original sequence.

189

2 5.39 13 4.35

3 4.61 14 5.41

4 4.36 15 7.05

5 3.92 16 9.10

6 3.36 17 11.64

7 2.87 18 14.94

8 2.84 19 19.10

9 2.97 20 24.29

10 2.95 21 30.71

The original sequence is divided into two parts: the first part accounting for almost eighty percent of the original sequence (the first 16 elements) is used to build model and denoted as X , and the rest is used

to validate extrapolation ability and denoted as Y.

191

4.1. Polynomial order selection, background coefficient search and parameter estimation

193 194

AC

190

192

11 3.03

CE

188

1 6.03 12 3.72

PT

Index Actual value Index Actual value

(0)

According to Algorithm 1, the modeling sequence X is divided into the training sequence Xtrain = (0) {6.03, 5.39, 4.61, 4.36, 3.92, 3.36, 2.87, 2.84, 2.97, 2.95, 3.03, 3.72}, and the testing sequence Xtest = {4.35, 5.41, 7.05, 9.10}.

195

Under the least squared error criterion (i.e., wk = 1) and least square percentage error criterion (i.e., wk =

196

1/x(0) (k)), the fitting error MAPEtrain and predicting error MAPEtest versus the increasing of background

197

coefficient for GPMB(1,1,N ) model with N = 0, 1, 2, 3 are displayed in Figure 2.

198 199

It can be seen in Figure 2 that the GPMB(1,1,N ) models under the least squared error criterion and the least squared percentage error criterion has the similar trend. 10

ACCEPTED MANUSCRIPT

9

(b)

0

0.6 0.2 0.4 0.8 background coefficient N=0

13

51.5

8

51

6

50.5

60 1

4

62

14

0

0.2 0.4 0.6 0.8 background coefficient N=1

MAPE

train

(%)

12 61

12

60

11

0

0.2 0.4 0.6 0.8 background coefficient

59 1

6

8

4

4

2

52

10

0

6

50.5

4

4

0

0.2 0.4 0.6 0.8 background coefficient

50 1

2

0 1

16 12

8

51

0.2 0.4 0.6 0.8 background coefficient N=2

8

8 6

10

12

50 1

51.5

10

8

4

0

0.4 0.2 0.6 0.8 background coefficient

0 1

150

10

100

5

50

(%)

10

N=3

15

test

52

16

0

0

0.2 0.4 0.6 0.8 background coefficient N=3

12

0 1

200

10

160

8

120

6

80

4

40

2

0

0.8 0.2 0.4 0.6 background coefficient

MAPE

61

12

N=2

10

(%)

MAPE

10

52.5

test

62

N=1

14

MAPE

11

train

(%)

63

CR IP T

N=0

12

(a)

0 1

AN US

Figure 2: Trend of the MAPEtrain and MAPEtest with the increasing of background coefficient for GPMB(1,1,N ) model under (a) the least squared error criterion and (b) the least squared percentage error criterion.

Under the least squared error criterion, the MAPEtrain of GPMB(1,1,0) model reaches the minimum

201

9.63% with λopt = 0.45 while its MAPEtest keeps decreasing; the MAPEtrain of GPMB(1,1,1) model reaches

202

the minimum 5.78% with λopt = 0.51 while its MAPEtest keeps increasing; the MAPEtrain and MAPEtest of

203

the GPMB(1,1,2) and GPMB(1,1,3) model both increase after their decreasing, and their MAPEtrain s reach

204

the minimums 2.79% and 2.70% respective with λopt = 0.51 and λopt = 0.49, meanwhile their corresponding

205

MAPEtest s are 1.25% and 4.89%.

M

200

Under the least squared percentage error criterion, both MAPEtrain and MAPEtest of GPMB(1,1,0) model

207

keep decreasing with the increasing of background coefficient; the MAPEtrain of GPMB(1,1,1) model reaches

208

the minimum 5.71% with λopt = 0.52 while its MAPEtest keeps increasing; the MAPEtrain and MAPEtest

209

of the GPMB(1,1,2) and GPMB(1,1,3) model both increase after their decreasing, and their MAPEtrain s

210

reach the minimums 2.70% and 2.64% with λopt = 0.50 and λopt = 0.53, respectively, meanwhile their

211

corresponding MAPEtest s are 3.35% and 7.58%.

PT

ED

206

To sum up, whichever optimization criterion is chosen, the GPMB(1,1,0) and GPMB(1,1,1) model are

213

underfitted while the GPMB(1,1,3) model is overfitted, and only the GPMB(1,1,2) model having minimal

215

216 217 218 219

predicting error MAPEtest is the best. 4.2. Numerical results and comparisons with other models According to the computational steps in Section 2.3, GPMB(1,1,2) models under the two above criteria

AC

214

CE

212

are both built on the sequence X , and the results are summarized in Table 2. In addition, several affine

transformations with different coefficients are applied to the accumulating sequence to further verify the conclusions in Theorem 1.

220

It can be seen in Table 2 that the GPMB(1,1,2) model under the least squared percentage error criterion

221

outperforms that under the least squared error criterion due to that the evaluation measure MAPE is in

222 223 224

step with the optimization criterion. For the sake of comparisons, the models below are all build on the sequence X under the least squared error criterion. The difference stepwise ratio criterion in reference [26] indicates that GPM(1,1,0) is the best 11

ACCEPTED MANUSCRIPT

Table 2: Parameter values and errors of GPMB(1,1,2) model under different optimization criteria.

Least squared error Least squared percentage error

ρ 1.0 0.1 0.1 1.0 0.1 0.1

ξ 0.0 0.0 -0.6 0.0 0.0 -0.6

λopt 0.51 0.51 0.51 0.50 0.50 0.50

α ˆ 0.2020 0.2020 0.2020 0.1896 0.1896 0.1896

βˆ0 6.2129 0.62129 0.74248 6.2795 0.62795 0.74173

βˆ2 0.0769 0.00769 0.00769 0.0772 0.00772 0.00772

βˆ1 -1.8868 -0.18868 -0.18868 -1.8458 -0.18458 -0.18458

γ0 -3.1708 -0.31708 -0.91708 -4.4153 -0.44153 -1.04153

γ1 5.5725 0.55725 0.55725 5.4422 0.54422 0.54422

γ2 -0.3807 -0.03807 -0.03807 -0.4069 -0.04069 -0.04069

cˆ 3.2902 0.32902 0.32902 4.4918 0.44918 0.44918

Errors (%) MAPEX =2.45 MAPEY =1.75 MAPEX =2.46 MAPEY =1.19

CR IP T

Criteria

225

model, and the fitted and predicted values are listed in Table 3. The results of NDGM(1,1), Verhulst and

226

CPR are also obtained in Table 3.

Table 3: Fitted and predicted values of GPMB(1,1,2) and other models. GPM(1,1,0)

6.55 6.92 7.30 7.71 8.14

NDGM(1,1) Values 6.03 5.11 4.63 4.29 4.04 3.86 3.73 3.63 3.56 3.51 3.48 3.45 3.43 3.42 3.41 3.40

APE(%) 0.00 5.24 0.49 1.64 3.03 14.82 29.85 27.87 19.96 19.07 14.73 7.25 21.12 36.83 51.67 62.64 19.76 70.84 77.30 82.26 86.06 88.98 81.09

Verhulst

Values 6.03 0.77 0.87 0.97 1.10 1.23 1.39 1.56 1.75 1.96 2.20 2.46 2.75 3.07 3.43 3.82

APE(%) 0.00 85.73 81.21 77.64 72.02 63.30 51.71 45.18 41.16 33.55 27.49 33.87 36.74 43.17 51.34 58.00 50.13 63.47 68.40 72.60 76.18 79.21 71.97

AN US

APE(%) 12.22 46.17 33.56 25.83 12.90 7.29 32.62 41.50 42.85 51.84 56.08 34.23 21.19 2.88 16.64 31.82 29.35 43.72 53.71 61.77 68.26 73.49 60.19

M

Values 6.77 2.90 3.06 3.23 3.41 3.60 3.81 4.02 4.24 4.48 4.73 4.99 5.27 5.57 5.88 6.20

ED

11.66 15.07 19.42 24.91 31.80

APE(%) 0.29 1.08 3.52 2.35 3.10 1.40 7.95 1.78 5.49 2.54 3.26 2.97 0.43 1.07 0.94 1.00 2.45 0.15 0.87 1.66 2.54 3.53 1.75

3.39 3.39 3.39 3.39 3.38

4.25 4.72 5.23 5.79 6.38

CPR

Values 5.89 5.40 4.87 4.34 3.84 3.40 3.05 2.81 2.74 2.84 3.16 3.73 4.57 5.72 7.21 9.08 11.34 14.04 17.20 20.86 25.05

APE(%) 2.36 0.17 5.40 0.39 2.05 1.12 5.76 0.93 8.59 3.87 4.12 0.17 4.81 5.45 2.27 0.25 2.98 2.62 6.40 11.02 16.43 22.60 11.81

It can be seen in Table 3 that on the modeling sequence X , the MAPE of GPMB(1,1,2) model is close

CE

227

Values 6.05 5.33 4.77 4.26 3.80 3.41 3.10 2.89 2.81 2.88 3.13 3.61 4.37 5.47 6.98 9.01

PT

Actual values 1 6.03 2 5.39 3 4.61 4 4.36 5 3.92 6 3.36 7 2.87 8 2.84 9 2.97 10 2.95 11 3.03 12 3.72 13 4.35 14 5.41 15 7.05 16 9.10 MAPE (%) 17 11.64 18 14.94 19 19.10 20 24.29 21 30.71 MAPE (%)

GPMB(1,1,2)

Index

228

to that of CPR model, and much smaller than that of GPM(1,1,0), NDGM(1,1) and Verhulst model; on

229

the validating sequence Y, the MAPE of CPR model is more than 6 times larger than that of GPMB(1,1,2)

231

232

model which is also the smallest among these models. So GPMB(1,1,2) model can catch the overall trend

AC

230

accurately and has the best performance in this example. 5. Application

233

The annual per capita electricity consumption (ApCEC) of China is collected from the official website

234

National Bureau of Statistics of China (http://data.stats.gov.cn/english/easyquery.htm?cn=C01), as shown

235

in Table 4. The sequence from 2000 to 2010 is used to build the models, and that from 2011 to 2015 is used

236

to validate accuracy.

12

ACCEPTED MANUSCRIPT

Table 4: Annual per capita electricity consumption (kilowatt hour) of China form 2000 to 2015. Year ApCEC Year ApCEC

2000 1066.9 2008 2607.6

2001 1157.6 2009 2781.7

2002 1286.0 2010 3134.8

2003 1477.1 2011 3497.0

2004 1695.2 2012 3684.2

2005 1913.0 2013 3993.0

2006 2180.6 2014 4132.9

2007 2482.2 2015 4231.0

Under the least squared percentage error criterion, the polynomial order and background coefficient are

238

respectively selected as Nopt = 2 and λopt = 0.49 according to Algorithm 1. In particular, there only exists

239

slight difference between GPMB(1,1,2) model (MAPEtrain = 0.39%, MAPEtest = 9.13%) and GPMB(1,1,3)

240

model (MAPEtrain = 0.35%, MAPEtest = 9.83%). But GPMB(1,1,3) is an overfitted model with good

241

fitting but bad predicting accuracy. In addition, the results of GPMB(1,1,2) model are also compared to the

242

GPM(1,1,N ), NDGM(1,1), Verhulst and CPR model. The fitted and predicted values are listed in Table 5

243

and plotted in Figure 3.

CR IP T

237

GPMB(1,1,2)

GPM(1,1,0)

NDGM(1,1)

Verhulst

CPR

ED

M

Values APE(%) Values APE(%) Values APE(%) Values APE(%) Values APE(%) 1068.3 0.13 1040.2 2.50 1066.9 0.00 1066.9 0.00 1055.6 1.06 1149.7 0.69 1228.3 6.11 1147.2 0.90 1195.2 3.25 1158.7 0.10 1292.6 0.51 1365.7 6.20 1315.7 2.31 1357.8 5.58 1304.8 1.46 1484.7 0.51 1518.4 2.80 1496.4 1.31 1525.2 3.25 1485.8 0.59 1700.4 0.31 1688.2 0.41 1690.1 0.30 1703.9 0.51 1694.2 0.06 1927.6 0.76 1877.1 1.88 1897.8 0.79 1897.3 0.82 1922.0 0.47 2160.2 0.94 2087.0 4.29 2120.5 2.75 2108.2 3.32 2161.5 0.88 2395.4 3.50 2320.4 6.52 2359.4 4.95 2339.0 5.77 2404.8 3.12 2631.9 0.93 2580.0 1.06 2615.4 0.30 2592.1 0.59 2644.2 1.40 2869.0 3.14 2868.5 3.12 2890.0 3.89 2870.1 3.18 2871.8 3.24 3106.4 0.91 3189.3 1.74 3184.4 1.58 3175.7 1.31 3079.9 1.75 1.12 3.33 1.73 2.51 1.28 3343.9 4.38 3546.0 1.40 3500.0 0.09 3512.0 0.43 3260.6 6.76 3581.5 2.79 3942.7 7.02 3838.5 4.19 3882.0 5.37 3406.1 7.55 3819.1 4.36 4383.6 9.78 4201.4 5.22 4289.4 7.42 3508.7 12.13 4056.7 1.84 4873.9 17.93 4590.5 11.07 4738.1 14.64 3560.5 13.85 4294.4 1.50 5419.0 28.08 5007.7 18.36 5232.2 23.66 3553.7 16.01 2.97 12.84 7.78 10.31 11.26

PT

Year Actual values 2000 1066.9 2001 1157.6 2002 1286.0 2003 1477.1 2004 1695.2 2005 1913.0 2006 2180.6 2007 2482.2 2008 2607.6 2009 2781.7 2010 3134.8 MAPE(%) 2011 3497.0 2012 3684.2 2013 3993.0 2014 4132.9 2015 4231.0 MAPE(%)

AN US

Table 5: Fitted and predicted values of annual per capita electricity consumption of China form 2000 to 2015.

It can be seen in Table 5 that the MAPEs on the fitting sequence (from 2000 to 2010) of these models are all less than 5%, which indicates that they have almost the same fitting performance. But these models

246

differ form each other in the MAPEs on the predicting sequence (from 2011 to 2015). It is shown in Figure

247

3 that the annual per capita electricity consumption from 2011 to 2015 are overestimated by GPM(1,1,0),

248

NDGM(1,1) and Verhulst model, and on the contrary they are underestimated by CPR model. Only

250

251

AC

249

CE

244 245

GPMB(1,1,2) model presents an accurate estimation of the annual per capita electricity consumption and thus outperforms the other four models in this case study. 6. Conclusion

252

Different form GM(1,1) model (or NGM(1,1) model) based on the hypothesis that the original sequence

253

is in accord with homogeneous (or non-homogeneous) exponential trend, GPMB(1,1,N ) model is proposed

254

for a more general sequence with including but not limited to homogeneous and non-homogeneous trend.

255

In this paper, the main contributions are briefly summarized as follows: 13

6000 5500 5000 4500 4000

Actual values GMPB(1,1,2) GMP(1,1,0) NDGM(1,1) Verhulst CPR

Sequence for buliding models ←

3500 3000 2500 2000 1500 1000 2000 2001 2002

2003 2004

2005 2006

2007 2008 Year

2009 2010

CR IP T

Fitted and predicted values of ApCEC (kilowatt hour)

ACCEPTED MANUSCRIPT

2011 2012

2013 2014

2015

AN US

Figure 3: Fitted and predicted values from various models.

256

1. The most commonly used three optimization criteria, i.e., the least squared error criterion, the least

257

squared percentage error criterion and the least absolute percentage error criterion in grey prediction

258

model, are uniformly represented by weighted least square method.

2. The algorithm for polynomial order selection, background search and nonlinear parameter estimation of GPMB(1,1,N ) model is proposed, and the numerical example and application indicate that it is

261

more robust than the conventional difference stepwise ratio criterion in reference [26].

M

259 260

3. The mathematical analysis indicates that GM(1,1), NGM(1,1) and GM(1,1,tα ) model are all special

263

forms of GPMB(1,1,N ) model, and that the multiple transformation for original sequence and the affine

264

transformation for accumulating sequence are proved to be independent of modeling performance.

ED

262

Although the numerical example and application have verified the effectiveness of GPMB(1,1,N ) model,

266

large-scale simulation experiments under different signal-to-noise ratios still need to be designed to quantify

267

the influence of noise level and sample size on the modeling performance in the future work.

268

other hand, it should be also noticed that the model parameters in Eq. (11) are estimated based on the

269

difference equation (9) while the fitted and predicted values are obtained according to the analytic solution

270

to whitening equation (4). The skip between these two equations may introduce new errors [27]. Therefore,

271

the feature work can be focused on constructing a model based on Eq. (3) directly

273

x(0) (k) = ax(1) (k − 1) + b0 + b1 k + · · · + bN k N ,

where the polynomial feature and the model parameters can be together obtained by using the `1 -regularized least square method [37], arg min a,b0 ,··· ,bN

274

On the

CE

AC

272

PT

265

n  X

k=2

x

(0)

(k) − x ˆ

(0)

! N 2 X (k) + ϑ |a| + |bi | i=0

and ϑ > 0 is the regularized parameter. 14

(30)

ACCEPTED MANUSCRIPT

275

Acknowledgement Thank the editors and the anonymous referees for their insightful comments to improve the paper. This

277

work is supported by National Natural Science Foundation of China under grant 71671090, Aeronautical

278

Science Foundation of China under grant 2016ZG52068 and Qinglan Project for excellent youth or middle-

279

aged academic leaders in Jiangsu Province, China.

280

References

281

References

289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331

AN US

288

M

287

ED

285 286

PT

284

[1] J. L. Deng, Control problems of grey systems, Systems & Control Letters 1 (1982) 288–294. [2] M. S. Yin, Fifteen years of grey system theory research: A historical review and bibliometric analysis, Expert Systems with Applications 40 (2013) 2767–2775. [3] S. F. Liu, Y. J. Yang, N. M. Xie, J. Forrest, New progress of grey system theory in the new millennium, Grey Systems: Theory & Application 6 (2016) 2–31. [4] N. M. Xie, R. Z. Wang, A historic review of grey forecasting models, Journal of Grey System 29 (2017) 1–29. [5] X. Li, Y. Tao, Y. Zheng, Model GM(1,1,β) and its applicable region, Grey Systems: Theory & Application 3 (2013) 266–275. [6] D. C. Li, C. Yeh, C. J. Chang, An improved grey-based approach for early manufacturing data forecasting, Computers & Industrial Engineering 57 (2009) 1161–1167. [7] D. C. Li, C. J. Chang, C. C. Chen, W. C. Chen, Forecasting short-term electricity consumption using the adaptive grey-based approachan asian case, Omega 40 (2012) 767–773. [8] Z. Zhao, J. Wang, J. Zhao, Z. Su, Using a grey model optimized by differential evolution algorithm to forecast the per capita annual net income of rural households in china, Omega 40 (2012) 525–532. [9] Y. G. Dang, S. F. Liu, K. J. Chen, The GM models that x(n) be taken as initial value, Kybernetes 33 (2004) 247–254. [10] T. L. Tien, A new grey prediction model FGM(1, 1), Mathematical & Computer Modelling 49 (2009) 1416–1426. [11] Y. H. Wang, Y. G. Dang, Y. Q. Li, S. F. Liu, An approach to increase prediction precision of GM(1,1) model based on optimization of the initial condition, Expert Systems with Applications 37 (2010) 5640–5644. [12] J. Xu, T. Tan, M. Tu, L. Qi, Improvement of grey models by least squares, Expert Systems with Applications 38 (2011) 13961–13966. [13] Y. H. Wang, Q. Liu, J. R. Tang, C. W. Bin, L. X. Zhong, Optimization approach of background value and initial item for improving prediction precision of GM(1,1) model, Journal of Systems Engineering and Electronics 25 (2014) 77–82. [14] L. F. Wu, S. F. Liu, L. G. Yao, S. L. Yan, D. L. Liu, Grey system model with the fractional order accumulation, Communications in Nonlinear Science & Numerical Simulation 18 (2013) 1775–1785. [15] L. F. Wu, S. F. Liu, Y. J. Yang, A gray model with a time varying weighted generating operator, IEEE Transactions on Systems, Man & Cybernetics: Systems 46 (2016) 427–433. [16] C. I. Chen, S. J. Huang, The necessary and sufficient condition for GM(1,1) grey prediction model, Applied Mathematics & Computation 219 (2013) 6152–6162. [17] J. Liu, X. P. Xiao, J. H. Guo, S. H. Mao, Error and its upper bound estimation between the solutions of GM(1,1) grey forecasting models, Applied Mathematics & Computation 246 (2014) 648–660. [18] X. P. Xiao, H. Guo, S. H. Mao, The modeling mechanism, extension and optimization of grey GM(1,1) model, Applied Mathematical Modelling 38 (2014) 1896–1910. [19] T. X. Yao, S. F. Liu, N. M. Xie, On the properties of small sample of GM(1,1) model, Applied Mathematical Modelling 33 (2009) 1894–1903. [20] L. F. Wu, S. F. Liu, L. G. Yao, S. L. Yan, The effect of sample size on the grey system model, Applied Mathematical Modelling 37 (2013) 6577–6583. [21] M. Tabaszewski, C. Cempel, Using a set of GM(1,1) models to predict values of diagnostic symptoms, Mechanical Systems & Signal Processing s5253 (2015) 416–425. [22] R. Guo, D. Guo, Random fuzzy variable foundation for grey differential equation modeling, Soft Computing 13 (2009) 185–201. [23] J. Cui, S. F. Liu, B. Zeng, N. M. Xie, A novel grey forecasting model and its optimization, Applied Mathematical Modelling 37 (2013) 4399–4406. [24] W. Y. Qian, Y. G. Dang, S. F. Liu, Grey GM(1,1,tα ) model with time power and its application, Systems Engineering — Theory & Practice 32 (2012) 2247–2252. [25] X. Ma, Y. S. Hu, Z. B. Liu, A novel kernel regularized nonhomogeneous grey model and its applications, Communications in Nonlinear Science & Numerical Simulation 48 (2016) 51–62. [26] D. Luo, B. L. Wei, Grey forecasting model with polynomial term and its optimization, Journal of Grey System 29 (2017) 58–69. [27] N. M. Xie, S. F. Liu, Discrete grey forecasting model and its optimization, Applied Mathematical Modelling 33 (2009) 1173–1186.

CE

283

AC

282

CR IP T

276

15

ACCEPTED MANUSCRIPT

339 340 341 342 343 344 345 346 347 348 349 350

CR IP T

338

AN US

337

M

336

ED

335

PT

334

[28] N. M. Xie, S. F. Liu, Y. J. Yang, C. Q. Yuan, On novel grey forecasting model based on non-homogeneous index sequence, Applied Mathematical Modelling 37 (2013) 5059–5068. [29] B. Zeng, C. M. Luo, S. F. Liu, Y. Bai, C. Li, Development of an optimization method for the GM(1,N ) model, Engineering Applications of Artificial Intelligence 55 (2016) 353–362. [30] B. Zeng, C. M. Luo, S. F. Liu, C. Li, A novel multi-variable grey forecasting model and its application in forecasting the amount of motor vehicles in Beijing, Computers & Industrial Engineering 101 (2016) 479–489. [31] X. Ma, Z. B. Liu, The kernel-based nonlinear multivariate grey model, Applied Mathematical Modelling 56 (2018) 217–238. [32] B. Zeng, C. Li, Improved multi-variable grey forecasting model with a dynamic background-value coefficient and its application, Computers & Industrial Engineering 118 (2018) 278–290. [33] C. M. Bishop, Pattern Recognition and Machine Learning (Information Science and Statistics), Springer-Verlag New York, Inc., 2006. [34] N. Jorge, J. W. Stephen, Numerical Optimization, Springer-Verlag New York, Inc., 2006. [35] K. Jing, Y. Z. Liu, Morbidity problem of GM(1,1) model, Control & Decision 31 (2016) 869–874. [36] C. I. Chen, H. L. Chen, S. P. Chen, Forecasting of foreign exchange rates of taiwan’s major trading partners by novel nonlinear Grey Bernoulli model NGBM(1,1), Communications in Nonlinear Science & Numerical Simulation 13 (2008) 1194–1204. [37] K. Fountoulakis, J. Gondzio, Performance of first- and second-order methods for `1 -regularized least squares problems, Computational Optimization & Applications 65 (2016) 605–635.

CE

333

AC

332

16